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