[301] | 1 | // ----------------------------------------------------------------------- |
---|
| 2 | // Statistics.cc: Member functions for the templated StatsContainer |
---|
| 3 | // class. |
---|
| 4 | // ----------------------------------------------------------------------- |
---|
| 5 | // Copyright (C) 2006, Matthew Whiting, ATNF |
---|
| 6 | // |
---|
| 7 | // This program is free software; you can redistribute it and/or modify it |
---|
| 8 | // under the terms of the GNU General Public License as published by the |
---|
| 9 | // Free Software Foundation; either version 2 of the License, or (at your |
---|
| 10 | // option) any later version. |
---|
| 11 | // |
---|
| 12 | // Duchamp is distributed in the hope that it will be useful, but WITHOUT |
---|
| 13 | // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or |
---|
| 14 | // FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
---|
| 15 | // for more details. |
---|
| 16 | // |
---|
| 17 | // You should have received a copy of the GNU General Public License |
---|
| 18 | // along with Duchamp; if not, write to the Free Software Foundation, |
---|
| 19 | // Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA |
---|
| 20 | // |
---|
| 21 | // Correspondence concerning Duchamp may be directed to: |
---|
| 22 | // Internet email: Matthew.Whiting [at] atnf.csiro.au |
---|
| 23 | // Postal address: Dr. Matthew Whiting |
---|
| 24 | // Australia Telescope National Facility, CSIRO |
---|
| 25 | // PO Box 76 |
---|
| 26 | // Epping NSW 1710 |
---|
| 27 | // AUSTRALIA |
---|
| 28 | // ----------------------------------------------------------------------- |
---|
[201] | 29 | #include <iostream> |
---|
[393] | 30 | #include <duchamp/Utils/Statistics.hh> |
---|
| 31 | #include <duchamp/Utils/utils.hh> |
---|
[190] | 32 | |
---|
[201] | 33 | namespace Statistics |
---|
[190] | 34 | { |
---|
[266] | 35 | |
---|
[271] | 36 | template <class T> |
---|
| 37 | float madfmToSigma(T madfm){ |
---|
[266] | 38 | return float(madfm)/correctionFactor; |
---|
[282] | 39 | } |
---|
[266] | 40 | template float madfmToSigma<int>(int madfm); |
---|
| 41 | template float madfmToSigma<long>(long madfm); |
---|
| 42 | template float madfmToSigma<float>(float madfm); |
---|
| 43 | template float madfmToSigma<double>(double madfm); |
---|
| 44 | |
---|
[222] | 45 | //-------------------------------------------------------------------- |
---|
[282] | 46 | float sigmaToMADFM(float sigma){ |
---|
[275] | 47 | return float(sigma)*correctionFactor; |
---|
[282] | 48 | } |
---|
[222] | 49 | //-------------------------------------------------------------------- |
---|
[275] | 50 | //-------------------------------------------------------------------- |
---|
[222] | 51 | |
---|
| 52 | template <class Type> |
---|
| 53 | StatsContainer<Type>::StatsContainer(const StatsContainer<Type>& s) |
---|
[201] | 54 | { |
---|
[271] | 55 | /** |
---|
[365] | 56 | * The copy constructor for the StatsContainer class. Just uses |
---|
| 57 | * the assignment operator. |
---|
[271] | 58 | */ |
---|
[365] | 59 | operator=(s); |
---|
[201] | 60 | } |
---|
| 61 | template StatsContainer<int>::StatsContainer(const StatsContainer<int>& s); |
---|
| 62 | template StatsContainer<long>::StatsContainer(const StatsContainer<long>& s); |
---|
| 63 | template StatsContainer<float>::StatsContainer(const StatsContainer<float>& s); |
---|
| 64 | template StatsContainer<double>::StatsContainer(const StatsContainer<double>& s); |
---|
| 65 | //-------------------------------------------------------------------- |
---|
[192] | 66 | |
---|
[222] | 67 | template <class Type> |
---|
[271] | 68 | StatsContainer<Type>& |
---|
| 69 | StatsContainer<Type>::operator= (const StatsContainer<Type>& s) |
---|
[201] | 70 | { |
---|
[271] | 71 | /** |
---|
| 72 | * The assignment operator for the StatsContainer class. |
---|
| 73 | */ |
---|
| 74 | |
---|
[270] | 75 | if(this == &s) return *this; |
---|
[201] | 76 | this->defined = s.defined; |
---|
| 77 | this->mean = s.mean; |
---|
| 78 | this->stddev = s.stddev; |
---|
| 79 | this->median = s.median; |
---|
| 80 | this->madfm = s.madfm; |
---|
| 81 | this->threshold = s.threshold; |
---|
| 82 | this->pThreshold = s.pThreshold; |
---|
| 83 | this->useRobust = s.useRobust; |
---|
| 84 | this->useFDR = s.useFDR; |
---|
[270] | 85 | return *this; |
---|
[201] | 86 | } |
---|
| 87 | template StatsContainer<int>& StatsContainer<int>::operator= (const StatsContainer<int>& s); |
---|
| 88 | template StatsContainer<long>& StatsContainer<long>::operator= (const StatsContainer<long>& s); |
---|
| 89 | template StatsContainer<float>& StatsContainer<float>::operator= (const StatsContainer<float>& s); |
---|
| 90 | template StatsContainer<double>& StatsContainer<double>::operator= (const StatsContainer<double>& s); |
---|
| 91 | //-------------------------------------------------------------------- |
---|
[192] | 92 | |
---|
[222] | 93 | template <class Type> |
---|
[271] | 94 | float StatsContainer<Type>::getThresholdSNR() |
---|
| 95 | { |
---|
| 96 | /** |
---|
| 97 | * The SNR is defined in terms of excess over the middle estimator |
---|
| 98 | * in units of the spread estimator. |
---|
| 99 | */ |
---|
| 100 | return (threshold - this->getMiddle())/this->getSpread(); |
---|
| 101 | } |
---|
| 102 | template float StatsContainer<int>::getThresholdSNR(); |
---|
| 103 | template float StatsContainer<long>::getThresholdSNR(); |
---|
| 104 | template float StatsContainer<float>::getThresholdSNR(); |
---|
| 105 | template float StatsContainer<double>::getThresholdSNR(); |
---|
| 106 | //-------------------------------------------------------------------- |
---|
| 107 | |
---|
| 108 | template <class Type> |
---|
| 109 | void StatsContainer<Type>::setThresholdSNR(float snr) |
---|
| 110 | { |
---|
| 111 | /** |
---|
| 112 | * The SNR is defined in terms of excess over the middle estimator |
---|
| 113 | * in units of the spread estimator. |
---|
| 114 | */ |
---|
| 115 | threshold=this->getMiddle() + snr*this->getSpread(); |
---|
| 116 | } |
---|
| 117 | template void StatsContainer<int>::setThresholdSNR(float snr); |
---|
| 118 | template void StatsContainer<long>::setThresholdSNR(float snr); |
---|
| 119 | template void StatsContainer<float>::setThresholdSNR(float snr); |
---|
| 120 | template void StatsContainer<double>::setThresholdSNR(float snr); |
---|
| 121 | //-------------------------------------------------------------------- |
---|
| 122 | |
---|
| 123 | template <class Type> |
---|
[481] | 124 | float StatsContainer<Type>::valueToSNR(float value) |
---|
[271] | 125 | { |
---|
| 126 | /** |
---|
| 127 | * The SNR is defined in terms of excess over the middle estimator |
---|
| 128 | * in units of the spread estimator. |
---|
| 129 | */ |
---|
| 130 | return (value - this->getMiddle())/this->getSpread(); |
---|
| 131 | } |
---|
[481] | 132 | template float StatsContainer<int>::valueToSNR(float value); |
---|
| 133 | template float StatsContainer<long>::valueToSNR(float value); |
---|
| 134 | template float StatsContainer<float>::valueToSNR(float value); |
---|
| 135 | template float StatsContainer<double>::valueToSNR(float value); |
---|
[271] | 136 | //-------------------------------------------------------------------- |
---|
[481] | 137 | |
---|
| 138 | template <class Type> |
---|
| 139 | float StatsContainer<Type>::snrToValue(float snr) |
---|
| 140 | { |
---|
| 141 | /** |
---|
| 142 | * The SNR is defined in terms of excess over the middle estimator |
---|
| 143 | * in units of the spread estimator. |
---|
| 144 | */ |
---|
| 145 | return snr * this->getSpread() + this->getMiddle(); |
---|
| 146 | } |
---|
| 147 | template float StatsContainer<int>::snrToValue(float value); |
---|
| 148 | template float StatsContainer<long>::snrToValue(float value); |
---|
| 149 | template float StatsContainer<float>::snrToValue(float value); |
---|
| 150 | template float StatsContainer<double>::snrToValue(float value); |
---|
| 151 | //-------------------------------------------------------------------- |
---|
[271] | 152 | |
---|
| 153 | template <class Type> |
---|
[475] | 154 | void StatsContainer<Type>::setMiddle(float middle) |
---|
| 155 | { |
---|
| 156 | /** |
---|
| 157 | * The middle value is determined by the StatsContainer::useRobust |
---|
| 158 | * flag -- it will be either the median (if true), or the mean (if |
---|
| 159 | * false). |
---|
| 160 | */ |
---|
| 161 | if(useRobust) this->median = Type(middle); |
---|
| 162 | else this->mean = middle; |
---|
| 163 | } |
---|
| 164 | template void StatsContainer<int>::setMiddle(float middle); |
---|
| 165 | template void StatsContainer<long>::setMiddle(float middle); |
---|
| 166 | template void StatsContainer<float>::setMiddle(float middle); |
---|
| 167 | template void StatsContainer<double>::setMiddle(float middle); |
---|
| 168 | //-------------------------------------------------------------------- |
---|
| 169 | |
---|
| 170 | template <class Type> |
---|
[271] | 171 | float StatsContainer<Type>::getMiddle() |
---|
| 172 | { |
---|
| 173 | /** |
---|
| 174 | * The middle value is determined by the StatsContainer::useRobust |
---|
| 175 | * flag -- it will be either the median (if true), or the mean (if |
---|
| 176 | * false). |
---|
| 177 | */ |
---|
[385] | 178 | if(useRobust) return float(this->median); |
---|
| 179 | else return this->mean; |
---|
[271] | 180 | } |
---|
| 181 | template float StatsContainer<int>::getMiddle(); |
---|
| 182 | template float StatsContainer<long>::getMiddle(); |
---|
| 183 | template float StatsContainer<float>::getMiddle(); |
---|
| 184 | template float StatsContainer<double>::getMiddle(); |
---|
| 185 | //-------------------------------------------------------------------- |
---|
| 186 | |
---|
| 187 | template <class Type> |
---|
[475] | 188 | void StatsContainer<Type>::setSpread(float spread){ |
---|
| 189 | /** |
---|
| 190 | * The spread value is set according to the |
---|
| 191 | * StatsContainer::useRobust flag -- it will be either the madfm |
---|
| 192 | * (if true), or the rms (if false). If robust, the spread value will be |
---|
| 193 | * converted to a madfm from an equivalent rms under the assumption of |
---|
| 194 | * Gaussianity, using the Statistics::sigmaToMADFM function. |
---|
| 195 | */ |
---|
| 196 | if(useRobust) this->madfm = Type(sigmaToMADFM(spread)); |
---|
| 197 | else this->stddev = spread; |
---|
| 198 | } |
---|
| 199 | template void StatsContainer<int>::setSpread(float spread); |
---|
| 200 | template void StatsContainer<long>::setSpread(float spread); |
---|
| 201 | template void StatsContainer<float>::setSpread(float spread); |
---|
| 202 | template void StatsContainer<double>::setSpread(float spread); |
---|
| 203 | //-------------------------------------------------------------------- |
---|
| 204 | |
---|
| 205 | template <class Type> |
---|
[271] | 206 | float StatsContainer<Type>::getSpread(){ |
---|
| 207 | /** |
---|
| 208 | * The spread value returned is determined by the |
---|
| 209 | * StatsContainer::useRobust flag -- it will be either the madfm |
---|
| 210 | * (if true), or the rms (if false). If robust, the madfm will be |
---|
| 211 | * converted to an equivalent rms under the assumption of |
---|
| 212 | * Gaussianity, using the Statistics::madfmToSigma function. |
---|
| 213 | */ |
---|
[385] | 214 | if(useRobust) return madfmToSigma(this->madfm); |
---|
| 215 | else return this->stddev; |
---|
[271] | 216 | } |
---|
| 217 | template float StatsContainer<int>::getSpread(); |
---|
| 218 | template float StatsContainer<long>::getSpread(); |
---|
| 219 | template float StatsContainer<float>::getSpread(); |
---|
| 220 | template float StatsContainer<double>::getSpread(); |
---|
| 221 | //-------------------------------------------------------------------- |
---|
| 222 | |
---|
| 223 | template <class Type> |
---|
[275] | 224 | void StatsContainer<Type>::scaleNoise(float scale) |
---|
| 225 | { |
---|
| 226 | /** |
---|
| 227 | * Multiply the noise parameters (stddev & madfm) by a given |
---|
| 228 | * factor, and adjust the threshold. |
---|
| 229 | */ |
---|
| 230 | float snr = (threshold - this->getMiddle())/this->getSpread(); |
---|
| 231 | this->madfm = Type(this->madfm*scale); |
---|
| 232 | this->stddev *= scale; |
---|
| 233 | this->threshold = this->getMiddle() + snr*this->getSpread(); |
---|
| 234 | } |
---|
| 235 | template void StatsContainer<int>::scaleNoise(float scale); |
---|
| 236 | template void StatsContainer<long>::scaleNoise(float scale); |
---|
| 237 | template void StatsContainer<float>::scaleNoise(float scale); |
---|
| 238 | template void StatsContainer<double>::scaleNoise(float scale); |
---|
| 239 | //-------------------------------------------------------------------- |
---|
| 240 | |
---|
| 241 | template <class Type> |
---|
[271] | 242 | float StatsContainer<Type>::getPValue(float value) |
---|
| 243 | { |
---|
| 244 | /** |
---|
| 245 | * Get the "probability", under the assumption of normality, of a |
---|
| 246 | * value occuring. |
---|
| 247 | * |
---|
| 248 | * It is defined by \f$0.5 \operatorname{erfc}(z/\sqrt{2})\f$, where |
---|
| 249 | * \f$z=(x-\mu)/\sigma\f$. We need the factor of 0.5 here, as we are |
---|
| 250 | * only considering the positive tail of the distribution -- we |
---|
| 251 | * don't care about negative detections. |
---|
| 252 | */ |
---|
| 253 | |
---|
| 254 | float zStat = (value - this->getMiddle()) / this->getSpread(); |
---|
| 255 | return 0.5 * erfc( zStat / M_SQRT2 ); |
---|
| 256 | } |
---|
| 257 | template float StatsContainer<int>::getPValue(float value); |
---|
| 258 | template float StatsContainer<long>::getPValue(float value); |
---|
| 259 | template float StatsContainer<float>::getPValue(float value); |
---|
| 260 | template float StatsContainer<double>::getPValue(float value); |
---|
| 261 | //-------------------------------------------------------------------- |
---|
| 262 | |
---|
| 263 | template <class Type> |
---|
| 264 | bool StatsContainer<Type>::isDetection(float value){ |
---|
| 265 | /** |
---|
| 266 | * Compares the value given to the correct threshold, depending on |
---|
| 267 | * the value of the StatsContainer::useFDR flag. |
---|
| 268 | */ |
---|
| 269 | if(useFDR) return (this->getPValue(value) < this->pThreshold); |
---|
| 270 | else return (value > this->threshold); |
---|
| 271 | } |
---|
| 272 | template bool StatsContainer<int>::isDetection(float value); |
---|
| 273 | template bool StatsContainer<long>::isDetection(float value); |
---|
| 274 | template bool StatsContainer<float>::isDetection(float value); |
---|
| 275 | template bool StatsContainer<double>::isDetection(float value); |
---|
| 276 | //-------------------------------------------------------------------- |
---|
| 277 | |
---|
| 278 | template <class Type> |
---|
[222] | 279 | void StatsContainer<Type>::calculate(Type *array, long size) |
---|
[201] | 280 | { |
---|
[222] | 281 | /** |
---|
| 282 | * Calculate all four statistics for all elements of a given |
---|
| 283 | * array. |
---|
| 284 | * |
---|
| 285 | * \param array The input data array. |
---|
| 286 | * \param size The length of the input array |
---|
| 287 | */ |
---|
[275] | 288 | // findNormalStats(array, size, this->mean, this->stddev); |
---|
| 289 | // findMedianStats(array, size, this->median, this->madfm); |
---|
| 290 | findAllStats(array,size,this->mean,this->stddev,this->median,this->madfm); |
---|
[201] | 291 | this->defined = true; |
---|
| 292 | } |
---|
[222] | 293 | template void StatsContainer<int>::calculate(int *array, long size); |
---|
| 294 | template void StatsContainer<long>::calculate(long *array, long size); |
---|
| 295 | template void StatsContainer<float>::calculate(float *array, long size); |
---|
| 296 | template void StatsContainer<double>::calculate(double *array, long size); |
---|
[201] | 297 | //-------------------------------------------------------------------- |
---|
[190] | 298 | |
---|
[222] | 299 | template <class Type> |
---|
[275] | 300 | void StatsContainer<Type>::calculate(Type *array, long size, bool *mask) |
---|
[201] | 301 | { |
---|
[222] | 302 | /** |
---|
| 303 | * Calculate all four statistics for a subset of a given |
---|
| 304 | * array. The subset is defined by an array of bool |
---|
| 305 | * variables. |
---|
| 306 | * |
---|
| 307 | * \param array The input data array. |
---|
| 308 | * \param size The length of the input array |
---|
[275] | 309 | * \param mask An array of the same length that says whether to |
---|
| 310 | * include each member of the array in the calculations. Use a |
---|
| 311 | * value if mask=true. |
---|
[222] | 312 | */ |
---|
[275] | 313 | // findNormalStats(array, size, mask, this->mean, this->stddev); |
---|
| 314 | // findMedianStats(array, size, mask, this->median, this->madfm); |
---|
| 315 | findAllStats(array, size, mask, |
---|
| 316 | this->mean, this->stddev, this->median, this->madfm); |
---|
[201] | 317 | this->defined = true; |
---|
| 318 | } |
---|
[275] | 319 | template void StatsContainer<int>::calculate(int *array, long size, bool *mask); |
---|
| 320 | template void StatsContainer<long>::calculate(long *array, long size, bool *mask); |
---|
| 321 | template void StatsContainer<float>::calculate(float *array, long size, bool *mask); |
---|
| 322 | template void StatsContainer<double>::calculate(double *array, long size, bool *mask); |
---|
[201] | 323 | //-------------------------------------------------------------------- |
---|
[190] | 324 | |
---|
[222] | 325 | template <class Type> |
---|
| 326 | std::ostream& operator<< (std::ostream& theStream, StatsContainer<Type> &s) |
---|
[201] | 327 | { |
---|
[221] | 328 | /** |
---|
| 329 | * Prints out the four key statistics to the requested stream. |
---|
| 330 | */ |
---|
[201] | 331 | theStream << "Mean = " << s.mean << "\t" |
---|
| 332 | << "Std.Dev. = " << s.stddev << "\n" |
---|
| 333 | << "Median = " << s.median << "\t" |
---|
| 334 | << "MADFM = " << s.madfm << "\n"; |
---|
[270] | 335 | return theStream; |
---|
[201] | 336 | } |
---|
| 337 | template std::ostream& operator<<<int> (std::ostream& theStream, StatsContainer<int> &s); |
---|
| 338 | template std::ostream& operator<<<long> (std::ostream& theStream, StatsContainer<long> &s); |
---|
| 339 | template std::ostream& operator<<<float> (std::ostream& theStream, StatsContainer<float> &s); |
---|
| 340 | template std::ostream& operator<<<double> (std::ostream& theStream, StatsContainer<double> &s); |
---|
| 341 | |
---|
| 342 | |
---|
| 343 | |
---|
| 344 | } // matches: namespace Statistics { |
---|