[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> |
---|
[190] | 30 | #include <Utils/Statistics.hh> |
---|
| 31 | #include <Utils/utils.hh> |
---|
| 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 | //-------------------------------------------------------------------- |
---|
[275] | 46 | template <class T> |
---|
[282] | 47 | float sigmaToMADFM(float sigma){ |
---|
[275] | 48 | return float(sigma)*correctionFactor; |
---|
[282] | 49 | } |
---|
[222] | 50 | //-------------------------------------------------------------------- |
---|
[275] | 51 | //-------------------------------------------------------------------- |
---|
[222] | 52 | |
---|
| 53 | template <class Type> |
---|
| 54 | StatsContainer<Type>::StatsContainer(const StatsContainer<Type>& s) |
---|
[201] | 55 | { |
---|
[271] | 56 | /** |
---|
[365] | 57 | * The copy constructor for the StatsContainer class. Just uses |
---|
| 58 | * the assignment operator. |
---|
[271] | 59 | */ |
---|
[365] | 60 | operator=(s); |
---|
[201] | 61 | } |
---|
| 62 | template StatsContainer<int>::StatsContainer(const StatsContainer<int>& s); |
---|
| 63 | template StatsContainer<long>::StatsContainer(const StatsContainer<long>& s); |
---|
| 64 | template StatsContainer<float>::StatsContainer(const StatsContainer<float>& s); |
---|
| 65 | template StatsContainer<double>::StatsContainer(const StatsContainer<double>& s); |
---|
| 66 | //-------------------------------------------------------------------- |
---|
[192] | 67 | |
---|
[222] | 68 | template <class Type> |
---|
[271] | 69 | StatsContainer<Type>& |
---|
| 70 | StatsContainer<Type>::operator= (const StatsContainer<Type>& s) |
---|
[201] | 71 | { |
---|
[271] | 72 | /** |
---|
| 73 | * The assignment operator for the StatsContainer class. |
---|
| 74 | */ |
---|
| 75 | |
---|
[270] | 76 | if(this == &s) return *this; |
---|
[201] | 77 | this->defined = s.defined; |
---|
| 78 | this->mean = s.mean; |
---|
| 79 | this->stddev = s.stddev; |
---|
| 80 | this->median = s.median; |
---|
| 81 | this->madfm = s.madfm; |
---|
| 82 | this->threshold = s.threshold; |
---|
| 83 | this->pThreshold = s.pThreshold; |
---|
| 84 | this->useRobust = s.useRobust; |
---|
| 85 | this->useFDR = s.useFDR; |
---|
[270] | 86 | return *this; |
---|
[201] | 87 | } |
---|
| 88 | template StatsContainer<int>& StatsContainer<int>::operator= (const StatsContainer<int>& s); |
---|
| 89 | template StatsContainer<long>& StatsContainer<long>::operator= (const StatsContainer<long>& s); |
---|
| 90 | template StatsContainer<float>& StatsContainer<float>::operator= (const StatsContainer<float>& s); |
---|
| 91 | template StatsContainer<double>& StatsContainer<double>::operator= (const StatsContainer<double>& s); |
---|
| 92 | //-------------------------------------------------------------------- |
---|
[192] | 93 | |
---|
[222] | 94 | template <class Type> |
---|
[271] | 95 | float StatsContainer<Type>::getThresholdSNR() |
---|
| 96 | { |
---|
| 97 | /** |
---|
| 98 | * The SNR is defined in terms of excess over the middle estimator |
---|
| 99 | * in units of the spread estimator. |
---|
| 100 | */ |
---|
| 101 | return (threshold - this->getMiddle())/this->getSpread(); |
---|
| 102 | } |
---|
| 103 | template float StatsContainer<int>::getThresholdSNR(); |
---|
| 104 | template float StatsContainer<long>::getThresholdSNR(); |
---|
| 105 | template float StatsContainer<float>::getThresholdSNR(); |
---|
| 106 | template float StatsContainer<double>::getThresholdSNR(); |
---|
| 107 | //-------------------------------------------------------------------- |
---|
| 108 | |
---|
| 109 | template <class Type> |
---|
| 110 | void StatsContainer<Type>::setThresholdSNR(float snr) |
---|
| 111 | { |
---|
| 112 | /** |
---|
| 113 | * The SNR is defined in terms of excess over the middle estimator |
---|
| 114 | * in units of the spread estimator. |
---|
| 115 | */ |
---|
| 116 | threshold=this->getMiddle() + snr*this->getSpread(); |
---|
| 117 | } |
---|
| 118 | template void StatsContainer<int>::setThresholdSNR(float snr); |
---|
| 119 | template void StatsContainer<long>::setThresholdSNR(float snr); |
---|
| 120 | template void StatsContainer<float>::setThresholdSNR(float snr); |
---|
| 121 | template void StatsContainer<double>::setThresholdSNR(float snr); |
---|
| 122 | //-------------------------------------------------------------------- |
---|
| 123 | |
---|
| 124 | template <class Type> |
---|
| 125 | float StatsContainer<Type>::getSNR(float value) |
---|
| 126 | { |
---|
| 127 | /** |
---|
| 128 | * The SNR is defined in terms of excess over the middle estimator |
---|
| 129 | * in units of the spread estimator. |
---|
| 130 | */ |
---|
| 131 | return (value - this->getMiddle())/this->getSpread(); |
---|
| 132 | } |
---|
| 133 | template float StatsContainer<int>::getSNR(float value); |
---|
| 134 | template float StatsContainer<long>::getSNR(float value); |
---|
| 135 | template float StatsContainer<float>::getSNR(float value); |
---|
| 136 | template float StatsContainer<double>::getSNR(float value); |
---|
| 137 | //-------------------------------------------------------------------- |
---|
| 138 | |
---|
| 139 | template <class Type> |
---|
| 140 | float StatsContainer<Type>::getMiddle() |
---|
| 141 | { |
---|
| 142 | /** |
---|
| 143 | * The middle value is determined by the StatsContainer::useRobust |
---|
| 144 | * flag -- it will be either the median (if true), or the mean (if |
---|
| 145 | * false). |
---|
| 146 | */ |
---|
[385] | 147 | if(useRobust) return float(this->median); |
---|
| 148 | else return this->mean; |
---|
[271] | 149 | } |
---|
| 150 | template float StatsContainer<int>::getMiddle(); |
---|
| 151 | template float StatsContainer<long>::getMiddle(); |
---|
| 152 | template float StatsContainer<float>::getMiddle(); |
---|
| 153 | template float StatsContainer<double>::getMiddle(); |
---|
| 154 | //-------------------------------------------------------------------- |
---|
| 155 | |
---|
| 156 | template <class Type> |
---|
| 157 | float StatsContainer<Type>::getSpread(){ |
---|
| 158 | /** |
---|
| 159 | * The spread value returned is determined by the |
---|
| 160 | * StatsContainer::useRobust flag -- it will be either the madfm |
---|
| 161 | * (if true), or the rms (if false). If robust, the madfm will be |
---|
| 162 | * converted to an equivalent rms under the assumption of |
---|
| 163 | * Gaussianity, using the Statistics::madfmToSigma function. |
---|
| 164 | */ |
---|
[385] | 165 | if(useRobust) return madfmToSigma(this->madfm); |
---|
| 166 | else return this->stddev; |
---|
[271] | 167 | } |
---|
| 168 | template float StatsContainer<int>::getSpread(); |
---|
| 169 | template float StatsContainer<long>::getSpread(); |
---|
| 170 | template float StatsContainer<float>::getSpread(); |
---|
| 171 | template float StatsContainer<double>::getSpread(); |
---|
| 172 | //-------------------------------------------------------------------- |
---|
| 173 | |
---|
| 174 | template <class Type> |
---|
[275] | 175 | void StatsContainer<Type>::scaleNoise(float scale) |
---|
| 176 | { |
---|
| 177 | /** |
---|
| 178 | * Multiply the noise parameters (stddev & madfm) by a given |
---|
| 179 | * factor, and adjust the threshold. |
---|
| 180 | */ |
---|
| 181 | float snr = (threshold - this->getMiddle())/this->getSpread(); |
---|
| 182 | this->madfm = Type(this->madfm*scale); |
---|
| 183 | this->stddev *= scale; |
---|
| 184 | this->threshold = this->getMiddle() + snr*this->getSpread(); |
---|
| 185 | } |
---|
| 186 | template void StatsContainer<int>::scaleNoise(float scale); |
---|
| 187 | template void StatsContainer<long>::scaleNoise(float scale); |
---|
| 188 | template void StatsContainer<float>::scaleNoise(float scale); |
---|
| 189 | template void StatsContainer<double>::scaleNoise(float scale); |
---|
| 190 | //-------------------------------------------------------------------- |
---|
| 191 | |
---|
| 192 | template <class Type> |
---|
[271] | 193 | float StatsContainer<Type>::getPValue(float value) |
---|
| 194 | { |
---|
| 195 | /** |
---|
| 196 | * Get the "probability", under the assumption of normality, of a |
---|
| 197 | * value occuring. |
---|
| 198 | * |
---|
| 199 | * It is defined by \f$0.5 \operatorname{erfc}(z/\sqrt{2})\f$, where |
---|
| 200 | * \f$z=(x-\mu)/\sigma\f$. We need the factor of 0.5 here, as we are |
---|
| 201 | * only considering the positive tail of the distribution -- we |
---|
| 202 | * don't care about negative detections. |
---|
| 203 | */ |
---|
| 204 | |
---|
| 205 | float zStat = (value - this->getMiddle()) / this->getSpread(); |
---|
| 206 | return 0.5 * erfc( zStat / M_SQRT2 ); |
---|
| 207 | } |
---|
| 208 | template float StatsContainer<int>::getPValue(float value); |
---|
| 209 | template float StatsContainer<long>::getPValue(float value); |
---|
| 210 | template float StatsContainer<float>::getPValue(float value); |
---|
| 211 | template float StatsContainer<double>::getPValue(float value); |
---|
| 212 | //-------------------------------------------------------------------- |
---|
| 213 | |
---|
| 214 | template <class Type> |
---|
| 215 | bool StatsContainer<Type>::isDetection(float value){ |
---|
| 216 | /** |
---|
| 217 | * Compares the value given to the correct threshold, depending on |
---|
| 218 | * the value of the StatsContainer::useFDR flag. |
---|
| 219 | */ |
---|
| 220 | if(useFDR) return (this->getPValue(value) < this->pThreshold); |
---|
| 221 | else return (value > this->threshold); |
---|
| 222 | } |
---|
| 223 | template bool StatsContainer<int>::isDetection(float value); |
---|
| 224 | template bool StatsContainer<long>::isDetection(float value); |
---|
| 225 | template bool StatsContainer<float>::isDetection(float value); |
---|
| 226 | template bool StatsContainer<double>::isDetection(float value); |
---|
| 227 | //-------------------------------------------------------------------- |
---|
| 228 | |
---|
| 229 | template <class Type> |
---|
[222] | 230 | void StatsContainer<Type>::calculate(Type *array, long size) |
---|
[201] | 231 | { |
---|
[222] | 232 | /** |
---|
| 233 | * Calculate all four statistics for all elements of a given |
---|
| 234 | * array. |
---|
| 235 | * |
---|
| 236 | * \param array The input data array. |
---|
| 237 | * \param size The length of the input array |
---|
| 238 | */ |
---|
[275] | 239 | // findNormalStats(array, size, this->mean, this->stddev); |
---|
| 240 | // findMedianStats(array, size, this->median, this->madfm); |
---|
| 241 | findAllStats(array,size,this->mean,this->stddev,this->median,this->madfm); |
---|
[201] | 242 | this->defined = true; |
---|
| 243 | } |
---|
[222] | 244 | template void StatsContainer<int>::calculate(int *array, long size); |
---|
| 245 | template void StatsContainer<long>::calculate(long *array, long size); |
---|
| 246 | template void StatsContainer<float>::calculate(float *array, long size); |
---|
| 247 | template void StatsContainer<double>::calculate(double *array, long size); |
---|
[201] | 248 | //-------------------------------------------------------------------- |
---|
[190] | 249 | |
---|
[222] | 250 | template <class Type> |
---|
[275] | 251 | void StatsContainer<Type>::calculate(Type *array, long size, bool *mask) |
---|
[201] | 252 | { |
---|
[222] | 253 | /** |
---|
| 254 | * Calculate all four statistics for a subset of a given |
---|
| 255 | * array. The subset is defined by an array of bool |
---|
| 256 | * variables. |
---|
| 257 | * |
---|
| 258 | * \param array The input data array. |
---|
| 259 | * \param size The length of the input array |
---|
[275] | 260 | * \param mask An array of the same length that says whether to |
---|
| 261 | * include each member of the array in the calculations. Use a |
---|
| 262 | * value if mask=true. |
---|
[222] | 263 | */ |
---|
[275] | 264 | // findNormalStats(array, size, mask, this->mean, this->stddev); |
---|
| 265 | // findMedianStats(array, size, mask, this->median, this->madfm); |
---|
| 266 | findAllStats(array, size, mask, |
---|
| 267 | this->mean, this->stddev, this->median, this->madfm); |
---|
[201] | 268 | this->defined = true; |
---|
| 269 | } |
---|
[275] | 270 | template void StatsContainer<int>::calculate(int *array, long size, bool *mask); |
---|
| 271 | template void StatsContainer<long>::calculate(long *array, long size, bool *mask); |
---|
| 272 | template void StatsContainer<float>::calculate(float *array, long size, bool *mask); |
---|
| 273 | template void StatsContainer<double>::calculate(double *array, long size, bool *mask); |
---|
[201] | 274 | //-------------------------------------------------------------------- |
---|
[190] | 275 | |
---|
[222] | 276 | template <class Type> |
---|
| 277 | std::ostream& operator<< (std::ostream& theStream, StatsContainer<Type> &s) |
---|
[201] | 278 | { |
---|
[221] | 279 | /** |
---|
| 280 | * Prints out the four key statistics to the requested stream. |
---|
| 281 | */ |
---|
[201] | 282 | theStream << "Mean = " << s.mean << "\t" |
---|
| 283 | << "Std.Dev. = " << s.stddev << "\n" |
---|
| 284 | << "Median = " << s.median << "\t" |
---|
| 285 | << "MADFM = " << s.madfm << "\n"; |
---|
[270] | 286 | return theStream; |
---|
[201] | 287 | } |
---|
| 288 | template std::ostream& operator<<<int> (std::ostream& theStream, StatsContainer<int> &s); |
---|
| 289 | template std::ostream& operator<<<long> (std::ostream& theStream, StatsContainer<long> &s); |
---|
| 290 | template std::ostream& operator<<<float> (std::ostream& theStream, StatsContainer<float> &s); |
---|
| 291 | template std::ostream& operator<<<double> (std::ostream& theStream, StatsContainer<double> &s); |
---|
| 292 | |
---|
| 293 | |
---|
| 294 | |
---|
| 295 | } // matches: namespace Statistics { |
---|