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