| Diploma Thesis Percolation Simulation C++ Sourcecode Documentation |
00001 // statistics.h 00002 // (C) 2000 Andreas Krueger 00003 // v2.5 - last change: 27.9.2001 00004 // 00005 // average mu 00006 // variance sigma^2 00007 // skewness alpha_3 00008 // kurtosis alpha_4 00009 // 00010 // of all elements of a list 00011 // AFTER these elements are each preprocessed by a given function 00012 // e.g. linear or bitone_notzero 00013 00014 #include <list> 00015 #include <float.h> // for DBL_MAX 00016 00017 #ifndef LOADED_MEASURE 00018 #include "measure.h" 00019 #endif 00020 00021 namespace statistics { 00022 00023 template <class ELEM> class manipdata; 00024 template <class RES,class ELEM> class statistical; 00025 template <class RES> class muvarskewkurt; 00026 class minmax; 00027 00028 /* // predefine them here is not possible using Digital-cxx compiler 00029 template <class RES,class ELEM> 00030 muvarskewkurt<RES> calculate_moments(ELEM (manipdata<ELEM>::*function)(ELEM,ELEM), ELEM fn_arg, 00031 std::list<ELEM> &elements); 00032 template <class RES> void sdev_of_mean(long n); 00033 */ 00034 template <class RES> ostream& operator<<(ostream& os, muvarskewkurt<RES> A); 00035 00036 00037 00038 // here it starts: 00039 00040 // the data preprocessing functions 00041 template <class ELEM> class manipdata { 00042 public: 00043 // this is the default: DO NOTHING 00044 ELEM linear(ELEM x, ELEM dummy) {return (x);}; 00045 00046 // special manipulation functions: 00047 ELEM square(ELEM x, ELEM dummy) {return (x*x);}; 00048 00049 ELEM notzero(ELEM x, ELEM dummy) {return (x!=0);}; 00050 ELEM bitone_notzero(ELEM x,ELEM dummy){return (x&1!=0);}; 00051 00052 ELEM bit_notzero(ELEM x, ELEM bit) {return ((x&(ELEM)pow((int)2,(int)bit-1))!=0);}; 00053 ELEM allbits_notzero(ELEM x, ELEM dim) {return (x==(ELEM)pow((int)2,(int)dim)-1);}; 00054 00055 ELEM shift(ELEM x, ELEM shift) {return (x+shift);}; 00056 }; 00057 00058 00059 00060 template <class RES, class ELEM> class statistical { 00061 public: 00062 // four moments -> four "outer" functions 00063 RES lineardis(ELEM x, RES mu){return (RES)x-mu;}; 00064 RES squaredis(ELEM x, RES mu){return pow((RES)x-mu,2);}; 00065 RES cubicdis (ELEM x, RES mu){return pow((RES)x-mu,3);}; 00066 RES pow4dis (ELEM x, RES mu){return pow((RES)x-mu,4);}; 00067 00068 // the core of the game: sum up the list and divide by size... 00069 // mu is the first moment (=linear expectation value) or zero 00070 // inner_fn(x,fn_arg) manipulates x before it is summed 00071 RES average_of_list(RES (statistical<RES,ELEM>::*outer_fn)(ELEM,RES), RES mu, 00072 ELEM (manipdata<ELEM>::*inner_fn)(ELEM,ELEM), ELEM fn_arg, 00073 std::list<ELEM> &L); 00074 // now do it again for the square distances from average 00075 RES standarddeviation_of_list(RES av, 00076 RES (statistical<RES,ELEM>::*outer_fn)(ELEM,RES), RES mu, 00077 ELEM (manipdata<ELEM>::*inner_fn)(ELEM,ELEM), ELEM fn_arg, 00078 std::list<ELEM> &L); 00079 00080 // this is used by average, variance, s3, s4 00081 measure<RES> general_moment(RES (statistical<RES,ELEM>::*distfn)(ELEM,RES), RES mu, 00082 ELEM (manipdata<ELEM>::*inner_fn)(ELEM,ELEM),ELEM fn_arg, 00083 std::list<ELEM> & elements); 00084 // the four moments 00085 measure<RES> average (ELEM (manipdata<ELEM>::*inner_fn)(ELEM,ELEM), ELEM fn_arg, 00086 std::list<ELEM>& elements); 00087 measure<RES> variance(ELEM (manipdata<ELEM>::*inner_fn)(ELEM,ELEM), ELEM fn_arg, 00088 std::list<ELEM>& elements, RES mu); 00089 measure<RES> s3 (ELEM (manipdata<ELEM>::*inner_fn)(ELEM,ELEM), ELEM fn_arg, 00090 std::list<ELEM>& elements, RES mu); 00091 measure<RES> s4 (ELEM (manipdata<ELEM>::*inner_fn)(ELEM,ELEM), ELEM fn_arg, 00092 std::list<ELEM>& elements, RES mu); 00093 // these two are "z-standardized" by 1/sigma^k 00094 measure<RES> skewness(ELEM (manipdata<ELEM>::*inner_fn)(ELEM,ELEM), ELEM fn_arg, 00095 std::list<ELEM>& elements, RES mu, RES sigma); 00096 measure<RES> kurtosis(ELEM (manipdata<ELEM>::*inner_fn)(ELEM,ELEM), ELEM fn_arg, 00097 std::list<ELEM>& elements, RES mu, RES sigma); 00098 00099 // skewness ( >0 long tail to the right = steep rise at the left side) 00100 // kurtosis ( >3 leptokurtic = peak higher than normal-distribution, not platykurtic) 00101 }; 00102 00103 00104 // --> AVERAGE <-- 00105 // The average & standard-deviation of a whole list are calculated. 00106 // The two preprocessing functions outer_fn and inner_fn 00107 // process each element before it is summed: 00108 // sum_of [ outer_fn( inner_fn(ELEMENT) ) ] 00109 // 00110 // Why? 00111 // The inner-function can be used to manipulate the stored element 00112 // eg. to test if the element is zero - and add the resulting bool 00113 // The outer-function is used to square, cubic,... the manipulated-element - mu 00114 // eg. to calculate the average of square(manip(elem) - mu) or cubic(manip(elem) - mu) 00115 // to get the higher moments variance, skewness, kurtosis 00116 00117 template <class RES, class ELEM> 00118 RES statistical<RES,ELEM>::average_of_list(RES (statistical<RES,ELEM>::*outer_fn)(ELEM,RES), RES mu, 00119 ELEM (manipdata<ELEM>::*inner_fn)(ELEM,ELEM), ELEM fn_arg, 00120 std::list<ELEM> &L){ 00121 manipdata<ELEM> *processdata=NULL; 00122 statistical<RES,ELEM> *distfn=NULL; 00123 std::list<ELEM>::iterator one; 00124 00125 RES sum=(RES)0; 00126 for (one=L.begin();one!=L.end();one++){ 00127 sum += (distfn->*outer_fn)((processdata->*inner_fn)(*one, fn_arg), mu); 00128 } 00129 return sum/(RES)L.size(); 00130 } 00131 00132 template <class RES, class ELEM> 00133 RES statistical<RES,ELEM>::standarddeviation_of_list(RES av, 00134 RES (statistical<RES,ELEM>::*outer_fn)(ELEM,RES), RES mu, 00135 ELEM (manipdata<ELEM>::*inner_fn)(ELEM,ELEM), ELEM fn_arg, 00136 std::list<ELEM> &L){ 00137 manipdata<ELEM> *processdata=NULL; 00138 statistical<RES,ELEM> *distfn=NULL; 00139 std::list<ELEM>::iterator one; 00140 00141 RES sdev=(RES)0; 00142 RES diff; 00143 for (one=L.begin();one!=L.end();one++){ 00144 diff=(distfn->*outer_fn)((processdata->*inner_fn)(*one, fn_arg), mu) - av ; 00145 sdev +=diff*diff; 00146 } 00147 if (L.size()<2) return sqrt(sdev); // changed in version2.2: 00148 return sqrt(sdev/(RES)(L.size()-1)); // 1/N -> 1/(N-1) 00149 } 00150 00151 // --> MOMENTS <-- 00152 // the first 4 moments are calculated 00153 // of a list of elements of type ELEM 00154 // returned is a tupel (average, stand.deviation) 00155 // of class measure<RES> 00156 00157 template <class RES, class ELEM> 00158 measure<RES> 00159 statistical<RES,ELEM>::general_moment( 00160 RES (statistical<RES,ELEM>::*distfn)(ELEM,RES), RES mu, 00161 ELEM (manipdata<ELEM>::*inner_fn)(ELEM,ELEM),ELEM fn_arg, 00162 std::list<ELEM> & elements) 00163 { 00164 statistical<RES,ELEM> *s=NULL; 00165 measure<RES> temp; 00166 temp.set_av (s->average_of_list(distfn, mu, inner_fn, fn_arg, elements)); 00167 temp.set_sdev(s->standarddeviation_of_list(temp.av(), distfn, mu, inner_fn, fn_arg, elements) ); 00168 return temp; 00169 } 00170 00171 template <class RES, class ELEM> 00172 measure<RES> 00173 statistical<RES,ELEM>::average (ELEM (manipdata<ELEM>::*inner_fn)(ELEM,ELEM),ELEM fn_arg, 00174 std::list<ELEM> & elements){ 00175 RES (statistical<RES,ELEM>::*distfn)(ELEM,RES)= &statistical<RES,ELEM>::lineardis; 00176 return general_moment(distfn, 0, inner_fn, fn_arg, elements); 00177 } 00178 00179 template <class RES, class ELEM> 00180 measure<RES> statistical<RES,ELEM>::variance(ELEM (manipdata<ELEM>::*inner_fn)(ELEM,ELEM),ELEM fn_arg, 00181 std::list<ELEM> & elements, RES mu){ 00182 RES (statistical<RES,ELEM>::*distfn)(ELEM,RES)= &statistical<RES,ELEM>::squaredis; 00183 return general_moment(distfn, mu, inner_fn, fn_arg, elements); 00184 } 00185 00186 // s3 and s4 are the third and fourth moments, 00187 // but not z-standardized yet: 00188 template <class RES, class ELEM> 00189 measure<RES> statistical<RES,ELEM>::s3(ELEM (manipdata<ELEM>::*inner_fn)(ELEM,ELEM),ELEM fn_arg, 00190 std::list<ELEM> & elements, RES mu){ 00191 RES (statistical<RES,ELEM>::*distfn)(ELEM,RES)= &statistical<RES,ELEM>::cubicdis; 00192 return general_moment(distfn, mu, inner_fn, fn_arg, elements); 00193 } 00194 00195 template <class RES, class ELEM> 00196 measure<RES> statistical<RES,ELEM>::s4(ELEM (manipdata<ELEM>::*inner_fn)(ELEM,ELEM), ELEM fn_arg, 00197 std::list<ELEM> & elements, RES mu){ 00198 RES (statistical<RES,ELEM>::*distfn)(ELEM,RES)= &statistical<RES,ELEM>::pow4dis; 00199 return general_moment(distfn, mu, inner_fn, fn_arg, elements); 00200 } 00201 00202 // skewness and kurtosis are z-standardized: 00203 template <class RES, class ELEM> 00204 measure<RES> statistical<RES,ELEM>::skewness(ELEM (manipdata<ELEM>::*inner_fn)(ELEM,ELEM), ELEM fn_arg, 00205 std::list<ELEM> & elements, RES mu, 00206 RES sigma){ // changed in v2.3: sigma not measure<RES> anymore 00207 if (sigma==0) return 0; 00208 else return s3(inner_fn, fn_arg, elements, mu)/pow(sigma,3); 00209 } 00210 00211 template <class RES, class ELEM> 00212 measure<RES> statistical<RES,ELEM>::kurtosis(ELEM (manipdata<ELEM>::*inner_fn)(ELEM,ELEM), ELEM fn_arg, 00213 std::list<ELEM> & elements, RES mu, 00214 RES sigma){ // changed in v2.3: sigma not measure<RES> anymore 00215 if (sigma==0) return 0; 00216 return s4(inner_fn, fn_arg, elements, mu)/pow(sigma,4); 00217 } 00218 00219 // end of class "statistical" 00220 00221 template<class RES, class ELEM> 00222 void calculate_all_moments(ELEM (manipdata<ELEM>::*function)(ELEM,ELEM), ELEM fn_arg, 00223 std::list<ELEM> &elements, 00224 measure<RES> &average, measure<RES> &variance, 00225 measure<RES> &skewness,measure<RES> &kurtosis) 00226 { 00227 statistical<RES, ELEM> *moments=NULL; 00228 average =moments->average (function, fn_arg, elements); 00229 variance=moments->variance(function, fn_arg, elements, average.av()); 00230 skewness=moments->skewness(function, fn_arg, elements, average.av(), sqrt(variance)); 00231 kurtosis=moments->kurtosis(function, fn_arg, elements, average.av(), sqrt(variance)); 00232 } 00233 00234 00235 template <class RES> class muvarskewkurt { 00236 public: 00237 measure<RES> mu; 00238 measure<RES> variance; 00239 measure<RES> skewness; 00240 measure<RES> kurtosis; 00241 void zero(); 00242 void sdev_of_mean(long N); 00243 }; 00244 template <class RES> 00245 void muvarskewkurt<RES>::zero (){ 00246 mu=(RES)0; 00247 variance=(RES)0; 00248 skewness=(RES)0; 00249 kurtosis=(RES)0; 00250 } 00251 template <class RES> ostream& operator<<(ostream& os, muvarskewkurt<RES> A){ 00252 os <<"mu="<<A.mu<<" var="<<A.variance<<" skw="<<A.skewness<<" krt="<<A.kurtosis; 00253 return os; 00254 } 00255 00256 template <class RES> // new in version 2.2 00257 void muvarskewkurt<RES>::sdev_of_mean(long n){ 00258 mu.set_sdev (mu.sdev() /sqrt((RES)n)); // manipulation! sdev of mean ~ 1/sqrt(nn) 00259 variance.set_sdev(variance.sdev()/sqrt((RES)n)); // manipulation! sdev of mean ~ 1/sqrt(nn) 00260 skewness.set_sdev(skewness.sdev()/sqrt((RES)n)); // manipulation! sdev of mean ~ 1/sqrt(nn) 00261 kurtosis.set_sdev(kurtosis.sdev()/sqrt((RES)n)); // manipulation! sdev of mean ~ 1/sqrt(nn) 00262 } 00263 00264 template <class RES, class ELEM> 00265 muvarskewkurt<RES> calculate_moments(ELEM (manipdata<ELEM>::*function)(ELEM,ELEM), 00266 ELEM fn_arg, 00267 std::list<ELEM> &elements, 00268 RES dummy) 00269 { 00270 statistical<RES, ELEM> *moments=NULL; 00271 muvarskewkurt<RES> temp; 00272 temp.mu =moments->average (function, fn_arg, elements); 00273 temp.variance=moments->variance(function, fn_arg, elements, temp.mu.av()); 00274 temp.skewness=moments->skewness(function, fn_arg, elements, temp.mu.av(), sqrt(temp.variance.av())); 00275 temp.kurtosis=moments->kurtosis(function, fn_arg, elements, temp.mu.av(), sqrt(temp.variance.av())); 00276 return temp; 00277 } 00278 00279 00280 void test_statistics(){ 00281 cout <<"\nHi!\nTest statistics with a list of random numbers."; 00282 cout <<"\nHow many:"<<flush; 00283 int N=0; 00284 while (N<1) cin>>N; 00285 std::list<float> L; 00286 for (int i=0;i<N;i++) L.push_back(rand()/(float)RAND_MAX); 00287 muvarskewkurt<double> A; 00288 int dummy=42; 00289 double dummy2=1.0; // for calculate_moments result type 00290 A=calculate_moments<double,float>(&manipdata<float>::linear, dummy, L, dummy2); 00291 cout <<"The distribution has these statistical properties:\n"; 00292 cout <<"(average mu, variance, skewness, kurtosis)\n" <<A<<endl; 00293 } 00294 00295 00296 class minmax { 00297 public: 00298 double min; 00299 double max; 00300 minmax(); 00301 minmax(double min_,double max_); 00302 void reset(); // has to be set to [min;max]=[+inf;-inf] before being used 00303 void set_if_more_extreme(double comparer); 00304 }; 00305 00306 minmax::minmax(){ 00307 max=-DBL_MAX; min=DBL_MAX; 00308 } 00309 void minmax::reset(){ 00310 max=-DBL_MAX; min=DBL_MAX; 00311 } 00312 00313 minmax::minmax(double min_,double max_){ min=min_; max=max_; } 00314 00315 void minmax::set_if_more_extreme(double comparer){ 00316 if (comparer>max) max=comparer; 00317 if (comparer<min) min=comparer; 00318 } 00319 00320 template <class T> 00321 T minof2(T a, T b){ return (a<b)?a:b; } 00322 template <class T> 00323 T maxof2(T a, T b){ return (a<b)?b:a; } 00324 00325 00326 } // end of namespace statistics
| Diploma Thesis Sourcecode
Documentation check out the text and the executable binaries |