CLHEP VERSION Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // -*- C++ -*- 00002 // $Id: testRandDists.cc,v 1.11 2011/07/11 15:55:45 garren Exp $ 00003 // ---------------------------------------------------------------------- 00004 00005 // ---------------------------------------------------------------------- 00006 // 00007 // testRandDists -- tests of the correctness of random distributions 00008 // 00009 // Usage: 00010 // testRandDists < testRandDists.dat > testRandDists.log 00011 // 00012 // Currently tested: 00013 // RandGauss 00014 // RandGeneral 00015 // 00016 // M. Fischler 5/17/99 Reconfigured to be suitable for use with 00017 // an automated validation script - will return 00018 // 0 if validation is OK, or a mask indicating 00019 // where problems were found. 00020 // M. Fischler 5/18/99 Added test for RandGeneral. 00021 // Evgenyi T. 5/20/99 Vetted for compilation on various CLHEP/CERN 00022 // platforms. 00023 // M. Fischler 5/26/99 Extended distribution test to intervals of .5 00024 // sigma and to moments up to the sixth. 00025 // M. Fischler 10/29/99 Added validation for RandPoisson. 00026 // M. Fischler 11/09/99 Made gammln static to avoid (harmless) 00027 // confusion with the gammln in RandPoisson. 00028 // M. Fischler 2/04/99 Added validation for the Q and T versions of 00029 // Poisson and Gauss 00030 // M. Fischler 11/04/04 Add kludge to gaussianTest to deal with 00031 // different behaviour under optimization on 00032 // some compilers (gcc 2.95.2) 00033 // This behaviour was only seen with stepwise 00034 // RandGeneral and appears to be solely a 00035 // function of the test program. 00036 // 00037 // ---------------------------------------------------------------------- 00038 00039 #include "CLHEP/Units/GlobalPhysicalConstants.h" // used to provoke shadowing warnings 00040 #include "CLHEP/Random/Randomize.h" 00041 #include "CLHEP/Random/RandGaussQ.h" 00042 #include "CLHEP/Random/RandGaussT.h" 00043 #include "CLHEP/Random/RandPoissonQ.h" 00044 #include "CLHEP/Random/RandPoissonT.h" 00045 #include "CLHEP/Random/RandSkewNormal.h" 00046 #include "CLHEP/Random/defs.h" 00047 #include <iostream> 00048 #include <iomanip> 00049 #include <cmath> // double abs() 00050 #include <stdlib.h> // int abs() 00051 #include <cstdlib> // for exit() 00052 00053 using std::cin; 00054 using std::cout; 00055 using std::cerr; 00056 using std::endl; 00057 using std::setprecision; 00058 using namespace CLHEP; 00059 //#ifndef _WIN32 00060 //using std::exp; 00061 //#endif 00062 00063 // Tolerance of deviation from expected results 00064 static const double REJECT = 4.0; 00065 00066 // Mask bits to form a word indicating which if any dists were "bad" 00067 static const int GaussBAD = 1 << 0; 00068 static const int GeneralBAD = 1 << 1; 00069 static const int PoissonBAD = 1 << 2; 00070 static const int GaussQBAD = 1 << 3; 00071 static const int GaussTBAD = 1 << 4; 00072 static const int PoissonQBAD = 1 << 5; 00073 static const int PoissonTBAD = 1 << 6; 00074 static const int SkewNormalBAD = 1 << 7; 00075 00076 00077 // ********************** 00078 // 00079 // SECTION I - General tools for the various tests 00080 // 00081 // ********************** 00082 00083 static 00084 double gammln(double x) { 00085 // Note: This uses the gammln algorith in Numerical Recipes. 00086 // In the "old" RandPoisson there is a slightly different algorithm, 00087 // which mathematically is identical to this one. The advantage of 00088 // the modified method is one fewer division by x (in exchange for 00089 // doing one subtraction of 1 from x). The advantage of the method 00090 // here comes when .00001 < x < .65: In this range, the alternate 00091 // method produces results which have errors 10-100 times those 00092 // of this method (though still less than 1.0E-10). If we package 00093 // either method, we should use the reflection formula (6.1.4) so 00094 // that the user can never get inaccurate results, even for x very 00095 // small. The test for x < 1 is as costly as a divide, but so be it. 00096 00097 double y, tmp, ser; 00098 static double c[6] = { 00099 76.18009172947146, 00100 -86.50532032941677, 00101 24.01409824083091, 00102 -1.231739572450155, 00103 0.001208650973866179, 00104 -0.000005395239384953 }; 00105 y = x; 00106 tmp = x + 5.5; 00107 tmp -= (x+.5)*std::log(tmp); 00108 ser = 1.000000000190015; 00109 for (int i = 0; i < 6; i++) { 00110 ser += c[i]/(++y); 00111 } 00112 double ans = (-tmp + std::log (std::sqrt(CLHEP::twopi)*ser/x)); 00113 return ans; 00114 } 00115 00116 static 00117 double gser(double a, double x) { 00118 const int ITMAX = 100; 00119 const double EPS = 1.0E-8; 00120 double ap = a; 00121 double sum = 1/a; 00122 double del = sum; 00123 for (int n=0; n < ITMAX; n++) { 00124 ap++; 00125 del *= x/ap; 00126 sum += del; 00127 if (std::fabs(del) < std::fabs(sum)*EPS) { 00128 return sum*std::exp(-x+a*std::log(x)-gammln(a)); 00129 } 00130 } 00131 cout << "Problem - inaccurate gser " << a << ", " << x << "\n"; 00132 return sum*std::exp(-x+a*std::log(x)-gammln(a)); 00133 } 00134 00135 static 00136 double gcf(double a, double x) { 00137 const int ITMAX = 100; 00138 const double EPS = 1.0E-8; 00139 const double VERYSMALL = 1.0E-100; 00140 double b = x+1-a; 00141 double c = 1/VERYSMALL; 00142 double d = 1/b; 00143 double h = d; 00144 for (int i = 1; i <= ITMAX; i++) { 00145 double an = -i*(i-a); 00146 b += 2; 00147 d = an*d + b; 00148 if (std::fabs(d) < VERYSMALL) d = VERYSMALL; 00149 c = b + an/c; 00150 if (std::fabs(c) < VERYSMALL) c = VERYSMALL; 00151 d = 1/d; 00152 double del = d*c; 00153 h *= del; 00154 if (std::fabs(del-1.0) < EPS) { 00155 return std::exp(-x+a*std::log(x)-gammln(a))*h; 00156 } 00157 } 00158 cout << "Problem - inaccurate gcf " << a << ", " << x << "\n"; 00159 return std::exp(-x+a*std::log(x)-gammln(a))*h; 00160 } 00161 00162 static 00163 double gammp (double a, double x) { 00164 if (x < a+1) { 00165 return gser(a,x); 00166 } else { 00167 return 1-gcf(a,x); 00168 } 00169 } 00170 00171 00172 // ********************** 00173 // 00174 // SECTION II - Validation of specific distributions 00175 // 00176 // ********************** 00177 00178 // ------------ 00179 // gaussianTest 00180 // ------------ 00181 00182 bool gaussianTest ( HepRandom & dist, double mu, 00183 double sigma, int nNumbers ) { 00184 00185 bool good = true; 00186 double worstSigma = 0; 00187 00188 // We will accumulate mean and moments up to the sixth, 00189 // The second moment should be sigma**2, the fourth 3 sigma**4, 00190 // the sixth 15 sigma**6. The expected variance in these is 00191 // (for the m-th moment with m even) (2m-1)!! (m-1)!!**2 / n 00192 // (for the m-th moment with m odd) (2m-1)!! m!!**2 / n 00193 // We also do a histogram with bins every half sigma. 00194 00195 double sumx = 0; 00196 double sumx2 = 0; 00197 double sumx3 = 0; 00198 double sumx4 = 0; 00199 double sumx5 = 0; 00200 double sumx6 = 0; 00201 int counts[11]; 00202 int ncounts[11]; 00203 int ciu; 00204 for (ciu = 0; ciu < 11; ciu++) { 00205 counts[ciu] = 0; 00206 ncounts[ciu] = 0; 00207 } 00208 00209 int oldprecision = cout.precision(); 00210 cout.precision(5); 00211 // hack so that gcc 4.3 puts x and u into memory instead of a register 00212 volatile double x; 00213 volatile double u; 00214 int ipr = nNumbers / 10 + 1; 00215 for (int ifire = 0; ifire < nNumbers; ifire++) { 00216 x = dist(); // We avoid fire() because that is not virtual 00217 // in HepRandom. 00218 if( x < mu - 12.0*sigma ) { 00219 cout << "x = " << x << "\n"; 00220 } 00221 if ( (ifire % ipr) == 0 ) { 00222 cout << ifire << endl; 00223 } 00224 sumx += x; 00225 sumx2 += x*x; 00226 sumx3 += x*x*x; 00227 sumx4 += x*x*x*x; 00228 sumx5 += x*x*x*x*x; 00229 sumx6 += x*x*x*x*x*x; 00230 u = (x - mu) / sigma; 00231 if ( u >= 0 ) { 00232 ciu = (int)(2*u); 00233 if (ciu>10) ciu = 10; 00234 counts[ciu]++; 00235 } else { 00236 ciu = (int)(2*(-u)); 00237 if (ciu>10) ciu = 10; 00238 ncounts[ciu]++; 00239 } 00240 } 00241 00242 double mean = sumx / nNumbers; 00243 double u2 = sumx2/nNumbers - mean*mean; 00244 double u3 = sumx3/nNumbers - 3*sumx2*mean/nNumbers + 2*mean*mean*mean; 00245 double u4 = sumx4/nNumbers - 4*sumx3*mean/nNumbers 00246 + 6*sumx2*mean*mean/nNumbers - 3*mean*mean*mean*mean; 00247 double u5 = sumx5/nNumbers - 5*sumx4*mean/nNumbers 00248 + 10*sumx3*mean*mean/nNumbers 00249 - 10*sumx2*mean*mean*mean/nNumbers 00250 + 4*mean*mean*mean*mean*mean; 00251 double u6 = sumx6/nNumbers - 6*sumx5*mean/nNumbers 00252 + 15*sumx4*mean*mean/nNumbers 00253 - 20*sumx3*mean*mean*mean/nNumbers 00254 + 15*sumx2*mean*mean*mean*mean/nNumbers 00255 - 5*mean*mean*mean*mean*mean*mean; 00256 00257 cout << "Mean (should be close to " << mu << "): " << mean << endl; 00258 cout << "Second moment (should be close to " << sigma*sigma << 00259 "): " << u2 << endl; 00260 cout << "Third moment (should be close to zero): " << u3 << endl; 00261 cout << "Fourth moment (should be close to " << 3*sigma*sigma*sigma*sigma << 00262 "): " << u4 << endl; 00263 cout << "Fifth moment (should be close to zero): " << u5 << endl; 00264 cout << "Sixth moment (should be close to " 00265 << 15*sigma*sigma*sigma*sigma*sigma*sigma 00266 << "): " << u6 << endl; 00267 00268 // For large N, the variance squared in the scaled 2nd, 3rd 4th 5th and 00269 // 6th moments are roughly 2/N, 6/N, 96/N, 720/N and 10170/N respectively. 00270 // Based on this, we can judge how many sigma a result represents: 00271 00272 double del1 = std::sqrt ( (double) nNumbers ) * std::abs(mean - mu) / sigma; 00273 double del2 = std::sqrt ( nNumbers/2.0 ) * std::abs(u2 - sigma*sigma) / (sigma*sigma); 00274 double del3 = std::sqrt ( nNumbers/6.0 ) * std::abs(u3) / (sigma*sigma*sigma); 00275 double sigma4 = sigma*sigma*sigma*sigma; 00276 double del4 = std::sqrt ( nNumbers/96.0 ) * std::abs(u4 - 3 * sigma4) / sigma4; 00277 double del5 = std::sqrt ( nNumbers/720.0 ) * std::abs(u5) / (sigma*sigma4); 00278 double del6 = std::sqrt ( nNumbers/10170.0 ) * std::abs(u6 - 15*sigma4*sigma*sigma) 00279 / (sigma4*sigma*sigma); 00280 00281 cout << " These represent " << 00282 del1 << ", " << del2 << ", " << del3 << ", \n" 00283 <<" " << del4 << ", " << del5 << ", " << del6 00284 <<"\n standard deviations from expectations\n"; 00285 if ( del1 > worstSigma ) worstSigma = del1; 00286 if ( del2 > worstSigma ) worstSigma = del2; 00287 if ( del3 > worstSigma ) worstSigma = del3; 00288 if ( del4 > worstSigma ) worstSigma = del4; 00289 if ( del5 > worstSigma ) worstSigma = del5; 00290 if ( del6 > worstSigma ) worstSigma = del6; 00291 00292 if ( del1 > REJECT || del2 > REJECT || del3 > REJECT 00293 || del4 > REJECT || del5 > REJECT || del6 > REJECT ) { 00294 cout << "REJECT hypothesis that this distribution is correct!!\n"; 00295 good = false; 00296 } 00297 00298 // The variance of the bin counts is given by a Poisson estimate (std::sqrt(npq)). 00299 00300 double table[11] = { // Table of integrated density in each range: 00301 .191462, // 0.0 - 0.5 sigma 00302 .149882, // 0.5 - 1.0 sigma 00303 .091848, // 1.0 - 1.5 sigma 00304 .044057, // 1.5 - 2.0 sigma 00305 .016540, // 2.0 - 2.5 sigma 00306 .004860, // 2.5 - 3.0 sigma 00307 .001117, // 3.0 - 3.5 sigma 00308 .000201, // 3.5 - 4.0 sigma 00309 2.83E-5, // 4.0 - 4.5 sigma 00310 3.11E-6, // 4.5 - 5.0 sigma 00311 3.87E-7 // 5.0 sigma and up 00312 }; 00313 00314 for (int m1 = 0; m1 < 11; m1++) { 00315 double expect = table[m1]*nNumbers; 00316 double sig = std::sqrt ( table[m1] * (1.0-table[m1]) * nNumbers ); 00317 cout.precision(oldprecision); 00318 cout << "Between " << m1/2.0 << " sigma and " 00319 << m1/2.0+.5 << " sigma (should be about " << expect << "):\n " 00320 << " " 00321 << ncounts[m1] << " negative and " << counts[m1] << " positive " << "\n"; 00322 cout.precision(5); 00323 double negSigs = std::abs ( ncounts[m1] - expect ) / sig; 00324 double posSigs = std::abs ( counts[m1] - expect ) / sig; 00325 cout << " These represent " << 00326 negSigs << " and " << posSigs << " sigma from expectations\n"; 00327 if ( negSigs > REJECT || posSigs > REJECT ) { 00328 cout << "REJECT hypothesis that this distribution is correct!!\n"; 00329 good = false; 00330 } 00331 if ( negSigs > worstSigma ) worstSigma = negSigs; 00332 if ( posSigs > worstSigma ) worstSigma = posSigs; 00333 } 00334 00335 cout << "\n The worst deviation encountered (out of about 25) was " 00336 << worstSigma << " sigma \n\n"; 00337 00338 cout.precision(oldprecision); 00339 00340 return good; 00341 00342 } // gaussianTest() 00343 00344 00345 00346 // ------------ 00347 // skewNormalTest 00348 // ------------ 00349 00350 bool skewNormalTest ( HepRandom & dist, double k, int nNumbers ) { 00351 00352 bool good = true; 00353 double worstSigma = 0; 00354 00355 // We will accumulate mean and moments up to the sixth, 00356 // The second moment should be sigma**2, the fourth 3 sigma**4. 00357 // The expected variance in these is 00358 // (for the m-th moment with m even) (2m-1)!! (m-1)!!**2 / n 00359 // (for the m-th moment with m odd) (2m-1)!! m!!**2 / n 00360 00361 double sumx = 0; 00362 double sumx2 = 0; 00363 double sumx3 = 0; 00364 double sumx4 = 0; 00365 double sumx5 = 0; 00366 double sumx6 = 0; 00367 00368 int oldprecision = cout.precision(); 00369 cout.precision(5); 00370 // hack so that gcc 4.3 puts x into memory instead of a register 00371 volatile double x; 00372 // calculate mean and sigma 00373 double delta = k / std::sqrt( 1 + k*k ); 00374 double mu = delta/std::sqrt(CLHEP::halfpi); 00375 double mom2 = 1.; 00376 double mom3 = 3*delta*(1-(delta*delta)/3.)/std::sqrt(CLHEP::halfpi); 00377 double mom4 = 3.; 00378 double mom5 = 15*delta*(1-2.*(delta*delta)/3.+(delta*delta*delta*delta)/5.)/std::sqrt(CLHEP::halfpi); 00379 double mom6 = 15.; 00380 00381 int ipr = nNumbers / 10 + 1; 00382 for (int ifire = 0; ifire < nNumbers; ifire++) { 00383 x = dist(); // We avoid fire() because that is not virtual 00384 // in HepRandom. 00385 if( x < mu - 12.0 ) { 00386 cout << "x = " << x << "\n"; 00387 } 00388 if ( (ifire % ipr) == 0 ) { 00389 cout << ifire << endl; 00390 } 00391 sumx += x; 00392 sumx2 += x*x; 00393 sumx3 += x*x*x; 00394 sumx4 += x*x*x*x; 00395 sumx5 += x*x*x*x*x; 00396 sumx6 += x*x*x*x*x*x; 00397 } 00398 00399 double mean = sumx / nNumbers; 00400 double u2 = sumx2/nNumbers; 00401 double u3 = sumx3/nNumbers; 00402 double u4 = sumx4/nNumbers; 00403 double u5 = sumx5/nNumbers; 00404 double u6 = sumx6/nNumbers; 00405 00406 cout << "Mean (should be close to " << mu << "): " << mean << endl; 00407 cout << "Second moment (should be close to " << mom2 << "): " << u2 << endl; 00408 cout << "Third moment (should be close to " << mom3 << "): " << u3 << endl; 00409 cout << "Fourth moment (should be close to " << mom4 << "): " << u4 << endl; 00410 cout << "Fifth moment (should be close to " << mom5 << "): " << u5 << endl; 00411 cout << "Sixth moment (should be close to " << mom6 << "): " << u6 << endl; 00412 00413 double del1 = std::sqrt ( (double) nNumbers ) * std::abs(mean - mu); 00414 double del2 = std::sqrt ( nNumbers/2.0 ) * std::abs(u2 - mom2); 00415 double del3 = std::sqrt ( nNumbers/(15.-mom3*mom3) ) * std::abs(u3 - mom3 ); 00416 double del4 = std::sqrt ( nNumbers/96.0 ) * std::abs(u4 - mom4); 00417 double del5 = std::sqrt ( nNumbers/(945.-mom5*mom5) ) * std::abs(u5 - mom5 ); 00418 double del6 = std::sqrt ( nNumbers/10170.0 ) * std::abs(u6 - mom6); 00419 00420 cout << " These represent " << 00421 del1 << ", " << del2 << ", " << del3 << ", \n" 00422 <<" " << del4 << ", " << del5 << ", " << del6 00423 <<"\n standard deviations from expectations\n"; 00424 if ( del1 > worstSigma ) worstSigma = del1; 00425 if ( del2 > worstSigma ) worstSigma = del2; 00426 if ( del3 > worstSigma ) worstSigma = del3; 00427 if ( del4 > worstSigma ) worstSigma = del4; 00428 if ( del5 > worstSigma ) worstSigma = del5; 00429 if ( del6 > worstSigma ) worstSigma = del6; 00430 00431 if ( del1 > REJECT || del2 > REJECT || del3 > REJECT || 00432 del4 > REJECT || del5 > REJECT || del6 > REJECT ) { 00433 cout << "REJECT hypothesis that this distribution is correct!!\n"; 00434 good = false; 00435 } 00436 00437 cout << "\n The worst deviation encountered (out of about 25) was " 00438 << worstSigma << " sigma \n\n"; 00439 00440 cout.precision(oldprecision); 00441 00442 return good; 00443 00444 } // skewNormalTest() 00445 00446 00447 // ------------ 00448 // poissonTest 00449 // ------------ 00450 00451 class poisson { 00452 double mu_; 00453 public: 00454 poisson(double mu) : mu_(mu) {} 00455 double operator()(int r) { 00456 double logAnswer = -mu_ + r*std::log(mu_) - gammln(r+1); 00457 return std::exp(logAnswer); 00458 } 00459 }; 00460 00461 double* createRefDist ( poisson pdist, int N, 00462 int MINBIN, int MAXBINS, int clumping, 00463 int& firstBin, int& lastBin ) { 00464 00465 // Create the reference distribution -- making sure there are more than 00466 // 20 points at each value. The entire tail will be rolled up into one 00467 // value (at each end). We shall end up with some range of bins starting 00468 // at 0 or more, and ending at MAXBINS-1 or less. 00469 00470 double * refdist = new double [MAXBINS]; 00471 00472 int c = 0; // c is the number of the clump, that is, the member number 00473 // of the refdist array. 00474 int ic = 0; // ic is the number within the clump; mod clumping 00475 int r = 0; // r is the value of the variate. 00476 00477 // Determine the first bin: at least 20 entries must be at the level 00478 // of that bin (so that we won't immediately dip belpw 20) but the number 00479 // to enter is cumulative up to that bin. 00480 00481 double start = 0; 00482 double binc; 00483 while ( c < MAXBINS ) { 00484 for ( ic=0, binc=0; ic < clumping; ic++, r++ ) { 00485 binc += pdist(r) * N; 00486 } 00487 start += binc; 00488 if (binc >= MINBIN) break; 00489 c++; 00490 if ( c > MAXBINS/3 ) { 00491 cout << "The number of samples supplied " << N << 00492 " is too small to set up a chi^2 to test this distribution.\n"; 00493 exit(-1); 00494 } 00495 } 00496 firstBin = c; 00497 refdist[firstBin] = start; 00498 c++; 00499 00500 // Fill all the other bins until one has less than 20 items. 00501 double next = 0; 00502 while ( c < MAXBINS ) { 00503 for ( ic=0, binc=0; ic < clumping; ic++, r++ ) { 00504 binc += pdist(r) * N; 00505 } 00506 next = binc; 00507 if (next < MINBIN) break; 00508 refdist[c] = next; 00509 c++; 00510 } 00511 00512 // Shove all the remaining items into last bin. 00513 lastBin = c-1; 00514 next += refdist[lastBin]; 00515 while ( c < MAXBINS ) { 00516 for ( ic=0, binc=0; ic < clumping; ic++, r++ ) { 00517 binc += pdist(r) * N; 00518 } 00519 next += binc; 00520 c++; 00521 } 00522 refdist[lastBin] = next; 00523 00524 return refdist; 00525 00526 } // createRefDist() 00527 00528 00529 bool poissonTest ( RandPoisson & dist, double mu, int N ) { 00530 00531 // Three tests will be done: 00532 // 00533 // A chi-squared test will be used to test the hypothesis that the 00534 // generated distribution of N numbers matches the proper Poisson distribution. 00535 // 00536 // The same test will be applied to the distribution of numbers "clumping" 00537 // together std::sqrt(mu) bins. This will detect small deviations over several 00538 // touching bins, when mu is not small. 00539 // 00540 // The mean and second moment are checked against their theoretical values. 00541 00542 bool good = true; 00543 00544 int clumping = int(std::sqrt(mu)); 00545 if (clumping <= 1) clumping = 2; 00546 const int MINBIN = 20; 00547 const int MAXBINS = 1000; 00548 int firstBin; 00549 int lastBin; 00550 int firstBin2; 00551 int lastBin2; 00552 00553 poisson pdist(mu); 00554 00555 double* refdist = createRefDist( pdist, N, 00556 MINBIN, MAXBINS, 1, firstBin, lastBin); 00557 double* refdist2 = createRefDist( pdist, N, 00558 MINBIN, MAXBINS, clumping, firstBin2, lastBin2); 00559 00560 // Now roll the random dists, treating the tails in the same way as we go. 00561 00562 double sum = 0; 00563 double moment = 0; 00564 00565 double* samples = new double [MAXBINS]; 00566 double* samples2 = new double [MAXBINS]; 00567 int r; 00568 for (r = 0; r < MAXBINS; r++) { 00569 samples[r] = 0; 00570 samples2[r] = 0; 00571 } 00572 00573 int r1; 00574 int r2; 00575 for (int i = 0; i < N; i++) { 00576 r = dist.fire(); 00577 sum += r; 00578 moment += (r - mu)*(r - mu); 00579 r1 = r; 00580 if (r1 < firstBin) r1 = firstBin; 00581 if (r1 > lastBin) r1 = lastBin; 00582 samples[r1] += 1; 00583 r2 = r/clumping; 00584 if (r2 < firstBin2) r2 = firstBin2; 00585 if (r2 > lastBin2) r2 = lastBin2; 00586 samples2[r2] += 1; 00587 } 00588 00589 // #ifdef DIAGNOSTIC 00590 int k; 00591 for (k = firstBin; k <= lastBin; k++) { 00592 cout << k << " " << samples[k] << " " << refdist[k] << " " << 00593 (samples[k]-refdist[k])*(samples[k]-refdist[k])/refdist[k] << "\n"; 00594 } 00595 cout << "----\n"; 00596 for (k = firstBin2; k <= lastBin2; k++) { 00597 cout << k << " " << samples2[k] << " " << refdist2[k] << "\n"; 00598 } 00599 // #endif // DIAGNOSTIC 00600 00601 // Now find chi^2 for samples[] to apply the first test 00602 00603 double chi2 = 0; 00604 for ( r = firstBin; r <= lastBin; r++ ) { 00605 double delta = (samples[r] - refdist[r]); 00606 chi2 += delta*delta/refdist[r]; 00607 } 00608 int degFreedom = (lastBin - firstBin + 1) - 1; 00609 00610 // and finally, p. Since we only care about it for small values, 00611 // and never care about it past the 10% level, we can use the approximations 00612 // CL(chi^2,n) = 1/std::sqrt(CLHEP::twopi) * ErrIntC ( y ) with 00613 // y = std::sqrt(2*chi2) - std::sqrt(2*n-1) 00614 // errIntC (y) = std::exp((-y^2)/2)/(y*std::sqrt(CLHEP::twopi)) 00615 00616 double pval; 00617 pval = 1.0 - gammp ( .5*degFreedom , .5*chi2 ); 00618 00619 cout << "Chi^2 is " << chi2 << " on " << degFreedom << " degrees of freedom." 00620 << " p = " << pval << "\n"; 00621 00622 delete[] refdist; 00623 delete[] samples; 00624 00625 // Repeat the chi^2 and p for the clumped sample, to apply the second test 00626 00627 chi2 = 0; 00628 for ( r = firstBin2; r <= lastBin2; r++ ) { 00629 double delta = (samples2[r] - refdist2[r]); 00630 chi2 += delta*delta/refdist2[r]; 00631 } 00632 degFreedom = (lastBin2 - firstBin2 + 1) - 1; 00633 double pval2; 00634 pval2 = 1.0 - gammp ( .5*degFreedom , .5*chi2 ); 00635 00636 cout << "Clumps: Chi^2 is " << chi2 << " on " << degFreedom << 00637 " degrees of freedom." << " p = " << pval2 << "\n"; 00638 00639 delete[] refdist2; 00640 delete[] samples2; 00641 00642 // Check out the mean and sigma to apply the third test 00643 00644 double mean = sum / N; 00645 double sigma = std::sqrt( moment / (N-1) ); 00646 00647 double deviationMean = std::fabs(mean - mu)/(std::sqrt(mu/N)); 00648 double expectedSigma2Variance = (2*N*mu*mu/(N-1) + mu) / N; 00649 double deviationSigma = std::fabs(sigma*sigma-mu)/std::sqrt(expectedSigma2Variance); 00650 00651 cout << "Mean (should be " << mu << ") is " << mean << "\n"; 00652 cout << "Sigma (should be " << std::sqrt(mu) << ") is " << sigma << "\n"; 00653 00654 cout << "These are " << deviationMean << " and " << deviationSigma << 00655 " standard deviations from expected values\n\n"; 00656 00657 // If either p-value for the chi-squared tests is less that .0001, or 00658 // either the mean or sigma are more than 3.5 standard deviations off, 00659 // then reject the validation. This would happen by chance one time 00660 // in 2000. Since we will be validating for several values of mu, the 00661 // net chance of false rejection remains acceptable. 00662 00663 if ( (pval < .0001) || (pval2 < .0001) || 00664 (deviationMean > 3.5) || (deviationSigma > 3.5) ) { 00665 good = false; 00666 cout << "REJECT this distributon!!!\n"; 00667 } 00668 00669 return good; 00670 00671 } // poissonTest() 00672 00673 00674 // ********************** 00675 // 00676 // SECTION III - Tests of each distribution class 00677 // 00678 // ********************** 00679 00680 // --------- 00681 // RandGauss 00682 // --------- 00683 00684 int testRandGauss() { 00685 00686 cout << "\n--------------------------------------------\n"; 00687 cout << "Test of RandGauss distribution \n\n"; 00688 00689 long seed; 00690 cout << "Please enter an integer seed: "; 00691 cin >> seed; cout << seed << "\n"; 00692 if (seed == 0) { 00693 cout << "Moving on to next test...\n"; 00694 return 0; 00695 } 00696 00697 int nNumbers; 00698 cout << "How many numbers should we generate: "; 00699 cin >> nNumbers; cout << nNumbers << "\n"; 00700 00701 double mu; 00702 double sigma; 00703 cout << "Enter mu: "; 00704 cin >> mu; cout << mu << "\n"; 00705 00706 cout << "Enter sigma: "; 00707 cin >> sigma; cout << sigma << "\n"; 00708 00709 cout << "\nInstantiating distribution utilizing TripleRand engine...\n"; 00710 TripleRand eng (seed); 00711 RandGauss dist (eng, mu, sigma); 00712 00713 cout << "\n Sample fire(): \n"; 00714 double x; 00715 00716 x = dist.fire(); 00717 cout << x; 00718 00719 cout << "\n Testing operator() ... \n"; 00720 00721 bool good = gaussianTest ( dist, mu, sigma, nNumbers ); 00722 00723 if (good) { 00724 return 0; 00725 } else { 00726 return GaussBAD; 00727 } 00728 00729 } // testRandGauss() 00730 00731 00732 00733 // --------- 00734 // SkewNormal 00735 // --------- 00736 00737 int testSkewNormal() { 00738 00739 cout << "\n--------------------------------------------\n"; 00740 cout << "Test of SkewNormal distribution \n\n"; 00741 00742 long seed; 00743 cout << "Please enter an integer seed: "; 00744 cin >> seed; cout << seed << "\n"; 00745 if (seed == 0) { 00746 cout << "Moving on to next test...\n"; 00747 return 0; 00748 } 00749 00750 int nNumbers; 00751 cout << "How many numbers should we generate: "; 00752 cin >> nNumbers; cout << nNumbers << "\n"; 00753 00754 double k; 00755 cout << "Enter k: "; 00756 cin >> k; cout << k << "\n"; 00757 00758 cout << "\nInstantiating distribution utilizing TripleRand engine...\n"; 00759 TripleRand eng (seed); 00760 RandSkewNormal dist (eng, k); 00761 00762 cout << "\n Sample fire(): \n"; 00763 double x; 00764 00765 x = dist.fire(); 00766 cout << x; 00767 00768 cout << "\n Testing operator() ... \n"; 00769 00770 bool good = skewNormalTest ( dist, k, nNumbers ); 00771 00772 if (good) { 00773 return 0; 00774 } else { 00775 return SkewNormalBAD; 00776 } 00777 00778 } // testSkewNormal() 00779 00780 00781 00782 // --------- 00783 // RandGaussT 00784 // --------- 00785 00786 int testRandGaussT() { 00787 00788 cout << "\n--------------------------------------------\n"; 00789 cout << "Test of RandGaussT distribution \n\n"; 00790 00791 long seed; 00792 cout << "Please enter an integer seed: "; 00793 cin >> seed; cout << seed << "\n"; 00794 if (seed == 0) { 00795 cout << "Moving on to next test...\n"; 00796 return 0; 00797 } 00798 00799 int nNumbers; 00800 cout << "How many numbers should we generate: "; 00801 cin >> nNumbers; cout << nNumbers << "\n"; 00802 00803 double mu; 00804 double sigma; 00805 cout << "Enter mu: "; 00806 cin >> mu; cout << mu << "\n"; 00807 00808 cout << "Enter sigma: "; 00809 cin >> sigma; cout << sigma << "\n"; 00810 00811 cout << "\nInstantiating distribution utilizing TripleRand engine...\n"; 00812 TripleRand eng (seed); 00813 RandGaussT dist (eng, mu, sigma); 00814 00815 cout << "\n Sample fire(): \n"; 00816 double x; 00817 00818 x = dist.fire(); 00819 cout << x; 00820 00821 cout << "\n Testing operator() ... \n"; 00822 00823 bool good = gaussianTest ( dist, mu, sigma, nNumbers ); 00824 00825 if (good) { 00826 return 0; 00827 } else { 00828 return GaussTBAD; 00829 } 00830 00831 } // testRandGaussT() 00832 00833 00834 00835 // --------- 00836 // RandGaussQ 00837 // --------- 00838 00839 int testRandGaussQ() { 00840 00841 cout << "\n--------------------------------------------\n"; 00842 cout << "Test of RandGaussQ distribution \n\n"; 00843 00844 long seed; 00845 cout << "Please enter an integer seed: "; 00846 cin >> seed; cout << seed << "\n"; 00847 if (seed == 0) { 00848 cout << "Moving on to next test...\n"; 00849 return 0; 00850 } 00851 00852 int nNumbers; 00853 cout << "How many numbers should we generate: "; 00854 cin >> nNumbers; cout << nNumbers << "\n"; 00855 00856 if (nNumbers >= 20000000) { 00857 cout << "With that many samples RandGaussQ need not pass validation...\n"; 00858 } 00859 00860 double mu; 00861 double sigma; 00862 cout << "Enter mu: "; 00863 cin >> mu; cout << mu << "\n"; 00864 00865 cout << "Enter sigma: "; 00866 cin >> sigma; cout << sigma << "\n"; 00867 00868 cout << "\nInstantiating distribution utilizing DualRand engine...\n"; 00869 DualRand eng (seed); 00870 RandGaussQ dist (eng, mu, sigma); 00871 00872 cout << "\n Sample fire(): \n"; 00873 double x; 00874 00875 x = dist.fire(); 00876 cout << x; 00877 00878 cout << "\n Testing operator() ... \n"; 00879 00880 bool good = gaussianTest ( dist, mu, sigma, nNumbers ); 00881 00882 if (good) { 00883 return 0; 00884 } else { 00885 return GaussQBAD; 00886 } 00887 00888 } // testRandGaussQ() 00889 00890 00891 // --------- 00892 // RandPoisson 00893 // --------- 00894 00895 int testRandPoisson() { 00896 00897 cout << "\n--------------------------------------------\n"; 00898 cout << "Test of RandPoisson distribution \n\n"; 00899 00900 long seed; 00901 cout << "Please enter an integer seed: "; 00902 cin >> seed; cout << seed << "\n"; 00903 if (seed == 0) { 00904 cout << "Moving on to next test...\n"; 00905 return 0; 00906 } 00907 00908 cout << "\nInstantiating distribution utilizing TripleRand engine...\n"; 00909 TripleRand eng (seed); 00910 00911 int nNumbers; 00912 cout << "How many numbers should we generate for each mu: "; 00913 cin >> nNumbers; cout << nNumbers << "\n"; 00914 00915 bool good = true; 00916 00917 while (true) { 00918 double mu; 00919 cout << "Enter a value for mu: "; 00920 cin >> mu; cout << mu << "\n"; 00921 if (mu == 0) break; 00922 00923 RandPoisson dist (eng, mu); 00924 00925 cout << "\n Sample fire(): \n"; 00926 double x; 00927 00928 x = dist.fire(); 00929 cout << x; 00930 00931 cout << "\n Testing operator() ... \n"; 00932 00933 bool this_good = poissonTest ( dist, mu, nNumbers ); 00934 if (!this_good) { 00935 cout << "\n Poisson distribution for mu = " << mu << " is incorrect!!!\n"; 00936 } 00937 good &= this_good; 00938 } // end of the while(true) 00939 00940 if (good) { 00941 return 0; 00942 } else { 00943 return PoissonBAD; 00944 } 00945 00946 } // testRandPoisson() 00947 00948 00949 // --------- 00950 // RandPoissonQ 00951 // --------- 00952 00953 int testRandPoissonQ() { 00954 00955 cout << "\n--------------------------------------------\n"; 00956 cout << "Test of RandPoissonQ distribution \n\n"; 00957 00958 long seed; 00959 cout << "Please enter an integer seed: "; 00960 cin >> seed; cout << seed << "\n"; 00961 if (seed == 0) { 00962 cout << "Moving on to next test...\n"; 00963 return 0; 00964 } 00965 00966 cout << "\nInstantiating distribution utilizing TripleRand engine...\n"; 00967 TripleRand eng (seed); 00968 00969 int nNumbers; 00970 cout << "How many numbers should we generate for each mu: "; 00971 cin >> nNumbers; cout << nNumbers << "\n"; 00972 00973 bool good = true; 00974 00975 while (true) { 00976 double mu; 00977 cout << "Enter a value for mu: "; 00978 cin >> mu; cout << mu << "\n"; 00979 if (mu == 0) break; 00980 00981 RandPoissonQ dist (eng, mu); 00982 00983 cout << "\n Sample fire(): \n"; 00984 double x; 00985 00986 x = dist.fire(); 00987 cout << x; 00988 00989 cout << "\n Testing operator() ... \n"; 00990 00991 bool this_good = poissonTest ( dist, mu, nNumbers ); 00992 if (!this_good) { 00993 cout << "\n Poisson distribution for mu = " << mu << " is incorrect!!!\n"; 00994 } 00995 good &= this_good; 00996 } // end of the while(true) 00997 00998 if (good) { 00999 return 0; 01000 } else { 01001 return PoissonQBAD; 01002 } 01003 01004 } // testRandPoissonQ() 01005 01006 01007 // --------- 01008 // RandPoissonT 01009 // --------- 01010 01011 int testRandPoissonT() { 01012 01013 cout << "\n--------------------------------------------\n"; 01014 cout << "Test of RandPoissonT distribution \n\n"; 01015 01016 long seed; 01017 cout << "Please enter an integer seed: "; 01018 cin >> seed; cout << seed << "\n"; 01019 if (seed == 0) { 01020 cout << "Moving on to next test...\n"; 01021 return 0; 01022 } 01023 01024 cout << "\nInstantiating distribution utilizing TripleRand engine...\n"; 01025 TripleRand eng (seed); 01026 01027 int nNumbers; 01028 cout << "How many numbers should we generate for each mu: "; 01029 cin >> nNumbers; cout << nNumbers << "\n"; 01030 01031 bool good = true; 01032 01033 while (true) { 01034 double mu; 01035 cout << "Enter a value for mu: "; 01036 cin >> mu; cout << mu << "\n"; 01037 if (mu == 0) break; 01038 01039 RandPoissonT dist (eng, mu); 01040 01041 cout << "\n Sample fire(): \n"; 01042 double x; 01043 01044 x = dist.fire(); 01045 cout << x; 01046 01047 cout << "\n Testing operator() ... \n"; 01048 01049 bool this_good = poissonTest ( dist, mu, nNumbers ); 01050 if (!this_good) { 01051 cout << "\n Poisson distribution for mu = " << mu << " is incorrect!!!\n"; 01052 } 01053 good &= this_good; 01054 } // end of the while(true) 01055 01056 if (good) { 01057 return 0; 01058 } else { 01059 return PoissonTBAD; 01060 } 01061 01062 } // testRandPoissonT() 01063 01064 01065 // ----------- 01066 // RandGeneral 01067 // ----------- 01068 01069 int testRandGeneral() { 01070 01071 cout << "\n--------------------------------------------\n"; 01072 cout << "Test of RandGeneral distribution (using a Gaussian shape)\n\n"; 01073 01074 bool good; 01075 01076 long seed; 01077 cout << "Please enter an integer seed: "; 01078 cin >> seed; cout << seed << "\n"; 01079 if (seed == 0) { 01080 cout << "Moving on to next test...\n"; 01081 return 0; 01082 } 01083 01084 int nNumbers; 01085 cout << "How many numbers should we generate: "; 01086 cin >> nNumbers; cout << nNumbers << "\n"; 01087 01088 double mu; 01089 double sigma; 01090 mu = .5; // Since randGeneral always ranges from 0 to 1 01091 sigma = .06; 01092 01093 cout << "Enter sigma: "; 01094 cin >> sigma; cout << sigma << "\n"; 01095 // We suggest sigma be .06. This leaves room for 8 sigma 01096 // in the distribution. If it is much smaller, the number 01097 // of bins necessary to expect a good match will increase. 01098 // If sigma is much larger, the cutoff before 5 sigma can 01099 // cause the Gaussian hypothesis to be rejected. At .14, for 01100 // example, the 4th moment is 7 sigma away from expectation. 01101 01102 int nBins; 01103 cout << "Enter nBins for stepwise pdf test: "; 01104 cin >> nBins; cout << nBins << "\n"; 01105 // We suggest at least 10000 bins; fewer would risk 01106 // false rejection because the step-function curve 01107 // does not match an actual Gaussian. At 10000 bins, 01108 // a million-hit test does not have the resolving power 01109 // to tell the boxy pdf from the true Gaussian. At 5000 01110 // bins, it does. 01111 01112 double xBins = nBins; 01113 double* aProbFunc = new double [nBins]; 01114 double x; 01115 for ( int iBin = 0; iBin < nBins; iBin++ ) { 01116 x = iBin / (xBins-1); 01117 aProbFunc [iBin] = std::exp ( - (x-mu)*(x-mu) / (2*sigma*sigma) ); 01118 } 01119 // Note that this pdf is not normalized; RandGeneral does that 01120 01121 cout << "\nInstantiating distribution utilizing Ranlux64 engine...\n"; 01122 Ranlux64Engine eng (seed, 3); 01123 01124 { // Open block for testing type 1 - step function pdf 01125 01126 RandGeneral dist (eng, aProbFunc, nBins, 1); 01127 delete[] aProbFunc; 01128 01129 double* garbage = new double[nBins]; 01130 // We wish to verify that deleting the pdf 01131 // after instantiating the engine is fine. 01132 for ( int gBin = 0; gBin < nBins; gBin++ ) { 01133 garbage [gBin] = 1; 01134 } 01135 01136 cout << "\n Sample fire(): \n"; 01137 01138 x = dist.fire(); 01139 cout << x; 01140 01141 cout << "\n Testing operator() ... \n"; 01142 01143 good = gaussianTest ( dist, mu, sigma, nNumbers ); 01144 01145 delete[] garbage; 01146 01147 } // Close block for testing type 1 - step function pdf 01148 // dist goes out of scope but eng is supposed to stick around; 01149 // by closing this block we shall verify that! 01150 01151 cout << "Enter nBins for linearized pdf test: "; 01152 cin >> nBins; cout << nBins << "\n"; 01153 // We suggest at least 1000 bins; fewer would risk 01154 // false rejection because the non-smooth curve 01155 // does not match an actual Gaussian. At 1000 bins, 01156 // a million-hit test does not resolve the non-smoothness; 01157 // at 300 bins it does. 01158 01159 xBins = nBins; 01160 aProbFunc = new double [nBins]; 01161 for ( int jBin = 0; jBin < nBins; jBin++ ) { 01162 x = jBin / (xBins-1); 01163 aProbFunc [jBin] = std::exp ( - (x-mu)*(x-mu) / (2*sigma*sigma) ); 01164 } 01165 // Note that this pdf is not normalized; RandGeneral does that 01166 01167 RandGeneral dist (eng, aProbFunc, nBins, 0); 01168 01169 cout << "\n Sample operator(): \n"; 01170 01171 x = dist(); 01172 cout << x; 01173 01174 cout << "\n Testing operator() ... \n"; 01175 01176 bool good2 = gaussianTest ( dist, mu, sigma, nNumbers ); 01177 good = good && good2; 01178 01179 if (good) { 01180 return 0; 01181 } else { 01182 return GeneralBAD; 01183 } 01184 01185 } // testRandGeneral() 01186 01187 01188 01189 01190 // ********************** 01191 // 01192 // SECTION IV - Main 01193 // 01194 // ********************** 01195 01196 int main() { 01197 01198 int mask = 0; 01199 01200 mask |= testRandGauss(); 01201 mask |= testRandGaussQ(); 01202 mask |= testRandGaussT(); 01203 01204 mask |= testRandGeneral(); 01205 01206 mask |= testRandPoisson(); 01207 mask |= testRandPoissonQ(); 01208 mask |= testRandPoissonT(); 01209 01210 mask |= testSkewNormal(); // k = 0 (gaussian) 01211 mask |= testSkewNormal(); // k = -2 01212 mask |= testSkewNormal(); // k = 1 01213 mask |= testSkewNormal(); // k = 5 01214 01215 return mask > 0 ? -mask : mask; 01216 } 01217