CLHEP VERSION Reference Documentation
   
CLHEP Home Page     CLHEP Documentation     CLHEP Bug Reports

testRandDists.cc
Go to the documentation of this file.
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