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

testRandDists.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 // $Id: testRandDists.cc,v 1.11 2011/07/11 15:55:45 garren Exp $
3 // ----------------------------------------------------------------------
4 
5 // ----------------------------------------------------------------------
6 //
7 // testRandDists -- tests of the correctness of random distributions
8 //
9 // Usage:
10 // testRandDists < testRandDists.dat > testRandDists.log
11 //
12 // Currently tested:
13 // RandGauss
14 // RandGeneral
15 //
16 // M. Fischler 5/17/99 Reconfigured to be suitable for use with
17 // an automated validation script - will return
18 // 0 if validation is OK, or a mask indicating
19 // where problems were found.
20 // M. Fischler 5/18/99 Added test for RandGeneral.
21 // Evgenyi T. 5/20/99 Vetted for compilation on various CLHEP/CERN
22 // platforms.
23 // M. Fischler 5/26/99 Extended distribution test to intervals of .5
24 // sigma and to moments up to the sixth.
25 // M. Fischler 10/29/99 Added validation for RandPoisson.
26 // M. Fischler 11/09/99 Made gammln static to avoid (harmless)
27 // confusion with the gammln in RandPoisson.
28 // M. Fischler 2/04/99 Added validation for the Q and T versions of
29 // Poisson and Gauss
30 // M. Fischler 11/04/04 Add kludge to gaussianTest to deal with
31 // different behaviour under optimization on
32 // some compilers (gcc 2.95.2)
33 // This behaviour was only seen with stepwise
34 // RandGeneral and appears to be solely a
35 // function of the test program.
36 //
37 // ----------------------------------------------------------------------
38 
39 #include "CLHEP/Units/GlobalPhysicalConstants.h" // used to provoke shadowing warnings
40 #include "CLHEP/Random/Randomize.h"
41 #include "CLHEP/Random/RandGaussQ.h"
42 #include "CLHEP/Random/RandGaussT.h"
43 #include "CLHEP/Random/RandPoissonQ.h"
44 #include "CLHEP/Random/RandPoissonT.h"
45 #include "CLHEP/Random/RandSkewNormal.h"
46 #include "CLHEP/Random/defs.h"
47 #include <iostream>
48 #include <iomanip>
49 #include <cmath> // double abs()
50 #include <stdlib.h> // int abs()
51 #include <cstdlib> // for exit()
52 
53 using std::cin;
54 using std::cout;
55 using std::cerr;
56 using std::endl;
57 using std::setprecision;
58 using namespace CLHEP;
59 //#ifndef _WIN32
60 //using std::exp;
61 //#endif
62 
63 // Tolerance of deviation from expected results
64 static const double REJECT = 4.0;
65 
66 // Mask bits to form a word indicating which if any dists were "bad"
67 static const int GaussBAD = 1 << 0;
68 static const int GeneralBAD = 1 << 1;
69 static const int PoissonBAD = 1 << 2;
70 static const int GaussQBAD = 1 << 3;
71 static const int GaussTBAD = 1 << 4;
72 static const int PoissonQBAD = 1 << 5;
73 static const int PoissonTBAD = 1 << 6;
74 static const int SkewNormalBAD = 1 << 7;
75 
76 
77 // **********************
78 //
79 // SECTION I - General tools for the various tests
80 //
81 // **********************
82 
83 static
84 double gammln(double x) {
85  // Note: This uses the gammln algorith in Numerical Recipes.
86  // In the "old" RandPoisson there is a slightly different algorithm,
87  // which mathematically is identical to this one. The advantage of
88  // the modified method is one fewer division by x (in exchange for
89  // doing one subtraction of 1 from x). The advantage of the method
90  // here comes when .00001 < x < .65: In this range, the alternate
91  // method produces results which have errors 10-100 times those
92  // of this method (though still less than 1.0E-10). If we package
93  // either method, we should use the reflection formula (6.1.4) so
94  // that the user can never get inaccurate results, even for x very
95  // small. The test for x < 1 is as costly as a divide, but so be it.
96 
97  double y, tmp, ser;
98  static double c[6] = {
99  76.18009172947146,
100  -86.50532032941677,
101  24.01409824083091,
102  -1.231739572450155,
103  0.001208650973866179,
104  -0.000005395239384953 };
105  y = x;
106  tmp = x + 5.5;
107  tmp -= (x+.5)*std::log(tmp);
108  ser = 1.000000000190015;
109  for (int i = 0; i < 6; i++) {
110  ser += c[i]/(++y);
111  }
112  double ans = (-tmp + std::log (std::sqrt(CLHEP::twopi)*ser/x));
113  return ans;
114 }
115 
116 static
117 double gser(double a, double x) {
118  const int ITMAX = 100;
119  const double EPS = 1.0E-8;
120  double ap = a;
121  double sum = 1/a;
122  double del = sum;
123  for (int n=0; n < ITMAX; n++) {
124  ap++;
125  del *= x/ap;
126  sum += del;
127  if (std::fabs(del) < std::fabs(sum)*EPS) {
128  return sum*std::exp(-x+a*std::log(x)-gammln(a));
129  }
130  }
131  cout << "Problem - inaccurate gser " << a << ", " << x << "\n";
132  return sum*std::exp(-x+a*std::log(x)-gammln(a));
133 }
134 
135 static
136 double gcf(double a, double x) {
137  const int ITMAX = 100;
138  const double EPS = 1.0E-8;
139  const double VERYSMALL = 1.0E-100;
140  double b = x+1-a;
141  double c = 1/VERYSMALL;
142  double d = 1/b;
143  double h = d;
144  for (int i = 1; i <= ITMAX; i++) {
145  double an = -i*(i-a);
146  b += 2;
147  d = an*d + b;
148  if (std::fabs(d) < VERYSMALL) d = VERYSMALL;
149  c = b + an/c;
150  if (std::fabs(c) < VERYSMALL) c = VERYSMALL;
151  d = 1/d;
152  double del = d*c;
153  h *= del;
154  if (std::fabs(del-1.0) < EPS) {
155  return std::exp(-x+a*std::log(x)-gammln(a))*h;
156  }
157  }
158  cout << "Problem - inaccurate gcf " << a << ", " << x << "\n";
159  return std::exp(-x+a*std::log(x)-gammln(a))*h;
160 }
161 
162 static
163 double gammp (double a, double x) {
164  if (x < a+1) {
165  return gser(a,x);
166  } else {
167  return 1-gcf(a,x);
168  }
169 }
170 
171 
172 // **********************
173 //
174 // SECTION II - Validation of specific distributions
175 //
176 // **********************
177 
178 // ------------
179 // gaussianTest
180 // ------------
181 
182 bool gaussianTest ( HepRandom & dist, double mu,
183  double sigma, int nNumbers ) {
184 
185  bool good = true;
186  double worstSigma = 0;
187 
188 // We will accumulate mean and moments up to the sixth,
189 // The second moment should be sigma**2, the fourth 3 sigma**4,
190 // the sixth 15 sigma**6. The expected variance in these is
191 // (for the m-th moment with m even) (2m-1)!! (m-1)!!**2 / n
192 // (for the m-th moment with m odd) (2m-1)!! m!!**2 / n
193 // We also do a histogram with bins every half sigma.
194 
195  double sumx = 0;
196  double sumx2 = 0;
197  double sumx3 = 0;
198  double sumx4 = 0;
199  double sumx5 = 0;
200  double sumx6 = 0;
201  int counts[11];
202  int ncounts[11];
203  int ciu;
204  for (ciu = 0; ciu < 11; ciu++) {
205  counts[ciu] = 0;
206  ncounts[ciu] = 0;
207  }
208 
209  int oldprecision = cout.precision();
210  cout.precision(5);
211  // hack so that gcc 4.3 puts x and u into memory instead of a register
212  volatile double x;
213  volatile double u;
214  int ipr = nNumbers / 10 + 1;
215  for (int ifire = 0; ifire < nNumbers; ifire++) {
216  x = dist(); // We avoid fire() because that is not virtual
217  // in HepRandom.
218  if( x < mu - 12.0*sigma ) {
219  cout << "x = " << x << "\n";
220  }
221  if ( (ifire % ipr) == 0 ) {
222  cout << ifire << endl;
223  }
224  sumx += x;
225  sumx2 += x*x;
226  sumx3 += x*x*x;
227  sumx4 += x*x*x*x;
228  sumx5 += x*x*x*x*x;
229  sumx6 += x*x*x*x*x*x;
230  u = (x - mu) / sigma;
231  if ( u >= 0 ) {
232  ciu = (int)(2*u);
233  if (ciu>10) ciu = 10;
234  counts[ciu]++;
235  } else {
236  ciu = (int)(2*(-u));
237  if (ciu>10) ciu = 10;
238  ncounts[ciu]++;
239  }
240  }
241 
242  double mean = sumx / nNumbers;
243  double u2 = sumx2/nNumbers - mean*mean;
244  double u3 = sumx3/nNumbers - 3*sumx2*mean/nNumbers + 2*mean*mean*mean;
245  double u4 = sumx4/nNumbers - 4*sumx3*mean/nNumbers
246  + 6*sumx2*mean*mean/nNumbers - 3*mean*mean*mean*mean;
247  double u5 = sumx5/nNumbers - 5*sumx4*mean/nNumbers
248  + 10*sumx3*mean*mean/nNumbers
249  - 10*sumx2*mean*mean*mean/nNumbers
250  + 4*mean*mean*mean*mean*mean;
251  double u6 = sumx6/nNumbers - 6*sumx5*mean/nNumbers
252  + 15*sumx4*mean*mean/nNumbers
253  - 20*sumx3*mean*mean*mean/nNumbers
254  + 15*sumx2*mean*mean*mean*mean/nNumbers
255  - 5*mean*mean*mean*mean*mean*mean;
256 
257  cout << "Mean (should be close to " << mu << "): " << mean << endl;
258  cout << "Second moment (should be close to " << sigma*sigma <<
259  "): " << u2 << endl;
260  cout << "Third moment (should be close to zero): " << u3 << endl;
261  cout << "Fourth moment (should be close to " << 3*sigma*sigma*sigma*sigma <<
262  "): " << u4 << endl;
263  cout << "Fifth moment (should be close to zero): " << u5 << endl;
264  cout << "Sixth moment (should be close to "
265  << 15*sigma*sigma*sigma*sigma*sigma*sigma
266  << "): " << u6 << endl;
267 
268  // For large N, the variance squared in the scaled 2nd, 3rd 4th 5th and
269  // 6th moments are roughly 2/N, 6/N, 96/N, 720/N and 10170/N respectively.
270  // Based on this, we can judge how many sigma a result represents:
271 
272  double del1 = std::sqrt ( (double) nNumbers ) * std::abs(mean - mu) / sigma;
273  double del2 = std::sqrt ( nNumbers/2.0 ) * std::abs(u2 - sigma*sigma) / (sigma*sigma);
274  double del3 = std::sqrt ( nNumbers/6.0 ) * std::abs(u3) / (sigma*sigma*sigma);
275  double sigma4 = sigma*sigma*sigma*sigma;
276  double del4 = std::sqrt ( nNumbers/96.0 ) * std::abs(u4 - 3 * sigma4) / sigma4;
277  double del5 = std::sqrt ( nNumbers/720.0 ) * std::abs(u5) / (sigma*sigma4);
278  double del6 = std::sqrt ( nNumbers/10170.0 ) * std::abs(u6 - 15*sigma4*sigma*sigma)
279  / (sigma4*sigma*sigma);
280 
281  cout << " These represent " <<
282  del1 << ", " << del2 << ", " << del3 << ", \n"
283  <<" " << del4 << ", " << del5 << ", " << del6
284  <<"\n standard deviations from expectations\n";
285  if ( del1 > worstSigma ) worstSigma = del1;
286  if ( del2 > worstSigma ) worstSigma = del2;
287  if ( del3 > worstSigma ) worstSigma = del3;
288  if ( del4 > worstSigma ) worstSigma = del4;
289  if ( del5 > worstSigma ) worstSigma = del5;
290  if ( del6 > worstSigma ) worstSigma = del6;
291 
292  if ( del1 > REJECT || del2 > REJECT || del3 > REJECT
293  || del4 > REJECT || del5 > REJECT || del6 > REJECT ) {
294  cout << "REJECT hypothesis that this distribution is correct!!\n";
295  good = false;
296  }
297 
298  // The variance of the bin counts is given by a Poisson estimate (std::sqrt(npq)).
299 
300  double table[11] = { // Table of integrated density in each range:
301  .191462, // 0.0 - 0.5 sigma
302  .149882, // 0.5 - 1.0 sigma
303  .091848, // 1.0 - 1.5 sigma
304  .044057, // 1.5 - 2.0 sigma
305  .016540, // 2.0 - 2.5 sigma
306  .004860, // 2.5 - 3.0 sigma
307  .001117, // 3.0 - 3.5 sigma
308  .000201, // 3.5 - 4.0 sigma
309  2.83E-5, // 4.0 - 4.5 sigma
310  3.11E-6, // 4.5 - 5.0 sigma
311  3.87E-7 // 5.0 sigma and up
312  };
313 
314  for (int m1 = 0; m1 < 11; m1++) {
315  double expect = table[m1]*nNumbers;
316  double sig = std::sqrt ( table[m1] * (1.0-table[m1]) * nNumbers );
317  cout.precision(oldprecision);
318  cout << "Between " << m1/2.0 << " sigma and "
319  << m1/2.0+.5 << " sigma (should be about " << expect << "):\n "
320  << " "
321  << ncounts[m1] << " negative and " << counts[m1] << " positive " << "\n";
322  cout.precision(5);
323  double negSigs = std::abs ( ncounts[m1] - expect ) / sig;
324  double posSigs = std::abs ( counts[m1] - expect ) / sig;
325  cout << " These represent " <<
326  negSigs << " and " << posSigs << " sigma from expectations\n";
327  if ( negSigs > REJECT || posSigs > REJECT ) {
328  cout << "REJECT hypothesis that this distribution is correct!!\n";
329  good = false;
330  }
331  if ( negSigs > worstSigma ) worstSigma = negSigs;
332  if ( posSigs > worstSigma ) worstSigma = posSigs;
333  }
334 
335  cout << "\n The worst deviation encountered (out of about 25) was "
336  << worstSigma << " sigma \n\n";
337 
338  cout.precision(oldprecision);
339 
340  return good;
341 
342 } // gaussianTest()
343 
344 
345 
346 // ------------
347 // skewNormalTest
348 // ------------
349 
350 bool skewNormalTest ( HepRandom & dist, double k, int nNumbers ) {
351 
352  bool good = true;
353  double worstSigma = 0;
354 
355 // We will accumulate mean and moments up to the sixth,
356 // The second moment should be sigma**2, the fourth 3 sigma**4.
357 // The expected variance in these is
358 // (for the m-th moment with m even) (2m-1)!! (m-1)!!**2 / n
359 // (for the m-th moment with m odd) (2m-1)!! m!!**2 / n
360 
361  double sumx = 0;
362  double sumx2 = 0;
363  double sumx3 = 0;
364  double sumx4 = 0;
365  double sumx5 = 0;
366  double sumx6 = 0;
367 
368  int oldprecision = cout.precision();
369  cout.precision(5);
370  // hack so that gcc 4.3 puts x into memory instead of a register
371  volatile double x;
372  // calculate mean and sigma
373  double delta = k / std::sqrt( 1 + k*k );
374  double mu = delta/std::sqrt(CLHEP::halfpi);
375  double mom2 = 1.;
376  double mom3 = 3*delta*(1-(delta*delta)/3.)/std::sqrt(CLHEP::halfpi);
377  double mom4 = 3.;
378  double mom5 = 15*delta*(1-2.*(delta*delta)/3.+(delta*delta*delta*delta)/5.)/std::sqrt(CLHEP::halfpi);
379  double mom6 = 15.;
380 
381  int ipr = nNumbers / 10 + 1;
382  for (int ifire = 0; ifire < nNumbers; ifire++) {
383  x = dist(); // We avoid fire() because that is not virtual
384  // in HepRandom.
385  if( x < mu - 12.0 ) {
386  cout << "x = " << x << "\n";
387  }
388  if ( (ifire % ipr) == 0 ) {
389  cout << ifire << endl;
390  }
391  sumx += x;
392  sumx2 += x*x;
393  sumx3 += x*x*x;
394  sumx4 += x*x*x*x;
395  sumx5 += x*x*x*x*x;
396  sumx6 += x*x*x*x*x*x;
397  }
398 
399  double mean = sumx / nNumbers;
400  double u2 = sumx2/nNumbers;
401  double u3 = sumx3/nNumbers;
402  double u4 = sumx4/nNumbers;
403  double u5 = sumx5/nNumbers;
404  double u6 = sumx6/nNumbers;
405 
406  cout << "Mean (should be close to " << mu << "): " << mean << endl;
407  cout << "Second moment (should be close to " << mom2 << "): " << u2 << endl;
408  cout << "Third moment (should be close to " << mom3 << "): " << u3 << endl;
409  cout << "Fourth moment (should be close to " << mom4 << "): " << u4 << endl;
410  cout << "Fifth moment (should be close to " << mom5 << "): " << u5 << endl;
411  cout << "Sixth moment (should be close to " << mom6 << "): " << u6 << endl;
412 
413  double del1 = std::sqrt ( (double) nNumbers ) * std::abs(mean - mu);
414  double del2 = std::sqrt ( nNumbers/2.0 ) * std::abs(u2 - mom2);
415  double del3 = std::sqrt ( nNumbers/(15.-mom3*mom3) ) * std::abs(u3 - mom3 );
416  double del4 = std::sqrt ( nNumbers/96.0 ) * std::abs(u4 - mom4);
417  double del5 = std::sqrt ( nNumbers/(945.-mom5*mom5) ) * std::abs(u5 - mom5 );
418  double del6 = std::sqrt ( nNumbers/10170.0 ) * std::abs(u6 - mom6);
419 
420  cout << " These represent " <<
421  del1 << ", " << del2 << ", " << del3 << ", \n"
422  <<" " << del4 << ", " << del5 << ", " << del6
423  <<"\n standard deviations from expectations\n";
424  if ( del1 > worstSigma ) worstSigma = del1;
425  if ( del2 > worstSigma ) worstSigma = del2;
426  if ( del3 > worstSigma ) worstSigma = del3;
427  if ( del4 > worstSigma ) worstSigma = del4;
428  if ( del5 > worstSigma ) worstSigma = del5;
429  if ( del6 > worstSigma ) worstSigma = del6;
430 
431  if ( del1 > REJECT || del2 > REJECT || del3 > REJECT ||
432  del4 > REJECT || del5 > REJECT || del6 > REJECT ) {
433  cout << "REJECT hypothesis that this distribution is correct!!\n";
434  good = false;
435  }
436 
437  cout << "\n The worst deviation encountered (out of about 25) was "
438  << worstSigma << " sigma \n\n";
439 
440  cout.precision(oldprecision);
441 
442  return good;
443 
444 } // skewNormalTest()
445 
446 
447 // ------------
448 // poissonTest
449 // ------------
450 
451 class poisson {
452  double mu_;
453  public:
454  poisson(double mu) : mu_(mu) {}
455  double operator()(int r) {
456  double logAnswer = -mu_ + r*std::log(mu_) - gammln(r+1);
457  return std::exp(logAnswer);
458  }
459 };
460 
461 double* createRefDist ( poisson pdist, int N,
462  int MINBIN, int MAXBINS, int clumping,
463  int& firstBin, int& lastBin ) {
464 
465  // Create the reference distribution -- making sure there are more than
466  // 20 points at each value. The entire tail will be rolled up into one
467  // value (at each end). We shall end up with some range of bins starting
468  // at 0 or more, and ending at MAXBINS-1 or less.
469 
470  double * refdist = new double [MAXBINS];
471 
472  int c = 0; // c is the number of the clump, that is, the member number
473  // of the refdist array.
474  int ic = 0; // ic is the number within the clump; mod clumping
475  int r = 0; // r is the value of the variate.
476 
477  // Determine the first bin: at least 20 entries must be at the level
478  // of that bin (so that we won't immediately dip belpw 20) but the number
479  // to enter is cumulative up to that bin.
480 
481  double start = 0;
482  double binc;
483  while ( c < MAXBINS ) {
484  for ( ic=0, binc=0; ic < clumping; ic++, r++ ) {
485  binc += pdist(r) * N;
486  }
487  start += binc;
488  if (binc >= MINBIN) break;
489  c++;
490  if ( c > MAXBINS/3 ) {
491  cout << "The number of samples supplied " << N <<
492  " is too small to set up a chi^2 to test this distribution.\n";
493  exit(-1);
494  }
495  }
496  firstBin = c;
497  refdist[firstBin] = start;
498  c++;
499 
500  // Fill all the other bins until one has less than 20 items.
501  double next = 0;
502  while ( c < MAXBINS ) {
503  for ( ic=0, binc=0; ic < clumping; ic++, r++ ) {
504  binc += pdist(r) * N;
505  }
506  next = binc;
507  if (next < MINBIN) break;
508  refdist[c] = next;
509  c++;
510  }
511 
512  // Shove all the remaining items into last bin.
513  lastBin = c-1;
514  next += refdist[lastBin];
515  while ( c < MAXBINS ) {
516  for ( ic=0, binc=0; ic < clumping; ic++, r++ ) {
517  binc += pdist(r) * N;
518  }
519  next += binc;
520  c++;
521  }
522  refdist[lastBin] = next;
523 
524  return refdist;
525 
526 } // createRefDist()
527 
528 
529 bool poissonTest ( RandPoisson & dist, double mu, int N ) {
530 
531 // Three tests will be done:
532 //
533 // A chi-squared test will be used to test the hypothesis that the
534 // generated distribution of N numbers matches the proper Poisson distribution.
535 //
536 // The same test will be applied to the distribution of numbers "clumping"
537 // together std::sqrt(mu) bins. This will detect small deviations over several
538 // touching bins, when mu is not small.
539 //
540 // The mean and second moment are checked against their theoretical values.
541 
542  bool good = true;
543 
544  int clumping = int(std::sqrt(mu));
545  if (clumping <= 1) clumping = 2;
546  const int MINBIN = 20;
547  const int MAXBINS = 1000;
548  int firstBin;
549  int lastBin;
550  int firstBin2;
551  int lastBin2;
552 
553  poisson pdist(mu);
554 
555  double* refdist = createRefDist( pdist, N,
556  MINBIN, MAXBINS, 1, firstBin, lastBin);
557  double* refdist2 = createRefDist( pdist, N,
558  MINBIN, MAXBINS, clumping, firstBin2, lastBin2);
559 
560  // Now roll the random dists, treating the tails in the same way as we go.
561 
562  double sum = 0;
563  double moment = 0;
564 
565  double* samples = new double [MAXBINS];
566  double* samples2 = new double [MAXBINS];
567  int r;
568  for (r = 0; r < MAXBINS; r++) {
569  samples[r] = 0;
570  samples2[r] = 0;
571  }
572 
573  int r1;
574  int r2;
575  for (int i = 0; i < N; i++) {
576  r = dist.fire();
577  sum += r;
578  moment += (r - mu)*(r - mu);
579  r1 = r;
580  if (r1 < firstBin) r1 = firstBin;
581  if (r1 > lastBin) r1 = lastBin;
582  samples[r1] += 1;
583  r2 = r/clumping;
584  if (r2 < firstBin2) r2 = firstBin2;
585  if (r2 > lastBin2) r2 = lastBin2;
586  samples2[r2] += 1;
587  }
588 
589 // #ifdef DIAGNOSTIC
590  int k;
591  for (k = firstBin; k <= lastBin; k++) {
592  cout << k << " " << samples[k] << " " << refdist[k] << " " <<
593  (samples[k]-refdist[k])*(samples[k]-refdist[k])/refdist[k] << "\n";
594  }
595  cout << "----\n";
596  for (k = firstBin2; k <= lastBin2; k++) {
597  cout << k << " " << samples2[k] << " " << refdist2[k] << "\n";
598  }
599 // #endif // DIAGNOSTIC
600 
601  // Now find chi^2 for samples[] to apply the first test
602 
603  double chi2 = 0;
604  for ( r = firstBin; r <= lastBin; r++ ) {
605  double delta = (samples[r] - refdist[r]);
606  chi2 += delta*delta/refdist[r];
607  }
608  int degFreedom = (lastBin - firstBin + 1) - 1;
609 
610  // and finally, p. Since we only care about it for small values,
611  // and never care about it past the 10% level, we can use the approximations
612  // CL(chi^2,n) = 1/std::sqrt(CLHEP::twopi) * ErrIntC ( y ) with
613  // y = std::sqrt(2*chi2) - std::sqrt(2*n-1)
614  // errIntC (y) = std::exp((-y^2)/2)/(y*std::sqrt(CLHEP::twopi))
615 
616  double pval;
617  pval = 1.0 - gammp ( .5*degFreedom , .5*chi2 );
618 
619  cout << "Chi^2 is " << chi2 << " on " << degFreedom << " degrees of freedom."
620  << " p = " << pval << "\n";
621 
622  delete[] refdist;
623  delete[] samples;
624 
625  // Repeat the chi^2 and p for the clumped sample, to apply the second test
626 
627  chi2 = 0;
628  for ( r = firstBin2; r <= lastBin2; r++ ) {
629  double delta = (samples2[r] - refdist2[r]);
630  chi2 += delta*delta/refdist2[r];
631  }
632  degFreedom = (lastBin2 - firstBin2 + 1) - 1;
633  double pval2;
634  pval2 = 1.0 - gammp ( .5*degFreedom , .5*chi2 );
635 
636  cout << "Clumps: Chi^2 is " << chi2 << " on " << degFreedom <<
637  " degrees of freedom." << " p = " << pval2 << "\n";
638 
639  delete[] refdist2;
640  delete[] samples2;
641 
642  // Check out the mean and sigma to apply the third test
643 
644  double mean = sum / N;
645  double sigma = std::sqrt( moment / (N-1) );
646 
647  double deviationMean = std::fabs(mean - mu)/(std::sqrt(mu/N));
648  double expectedSigma2Variance = (2*N*mu*mu/(N-1) + mu) / N;
649  double deviationSigma = std::fabs(sigma*sigma-mu)/std::sqrt(expectedSigma2Variance);
650 
651  cout << "Mean (should be " << mu << ") is " << mean << "\n";
652  cout << "Sigma (should be " << std::sqrt(mu) << ") is " << sigma << "\n";
653 
654  cout << "These are " << deviationMean << " and " << deviationSigma <<
655  " standard deviations from expected values\n\n";
656 
657  // If either p-value for the chi-squared tests is less that .0001, or
658  // either the mean or sigma are more than 3.5 standard deviations off,
659  // then reject the validation. This would happen by chance one time
660  // in 2000. Since we will be validating for several values of mu, the
661  // net chance of false rejection remains acceptable.
662 
663  if ( (pval < .0001) || (pval2 < .0001) ||
664  (deviationMean > 3.5) || (deviationSigma > 3.5) ) {
665  good = false;
666  cout << "REJECT this distributon!!!\n";
667  }
668 
669  return good;
670 
671 } // poissonTest()
672 
673 
674 // **********************
675 //
676 // SECTION III - Tests of each distribution class
677 //
678 // **********************
679 
680 // ---------
681 // RandGauss
682 // ---------
683 
685 
686  cout << "\n--------------------------------------------\n";
687  cout << "Test of RandGauss distribution \n\n";
688 
689  long seed;
690  cout << "Please enter an integer seed: ";
691  cin >> seed; cout << seed << "\n";
692  if (seed == 0) {
693  cout << "Moving on to next test...\n";
694  return 0;
695  }
696 
697  int nNumbers;
698  cout << "How many numbers should we generate: ";
699  cin >> nNumbers; cout << nNumbers << "\n";
700 
701  double mu;
702  double sigma;
703  cout << "Enter mu: ";
704  cin >> mu; cout << mu << "\n";
705 
706  cout << "Enter sigma: ";
707  cin >> sigma; cout << sigma << "\n";
708 
709  cout << "\nInstantiating distribution utilizing TripleRand engine...\n";
710  TripleRand eng (seed);
711  RandGauss dist (eng, mu, sigma);
712 
713  cout << "\n Sample fire(): \n";
714  double x;
715 
716  x = dist.fire();
717  cout << x;
718 
719  cout << "\n Testing operator() ... \n";
720 
721  bool good = gaussianTest ( dist, mu, sigma, nNumbers );
722 
723  if (good) {
724  return 0;
725  } else {
726  return GaussBAD;
727  }
728 
729 } // testRandGauss()
730 
731 
732 
733 // ---------
734 // SkewNormal
735 // ---------
736 
738 
739  cout << "\n--------------------------------------------\n";
740  cout << "Test of SkewNormal distribution \n\n";
741 
742  long seed;
743  cout << "Please enter an integer seed: ";
744  cin >> seed; cout << seed << "\n";
745  if (seed == 0) {
746  cout << "Moving on to next test...\n";
747  return 0;
748  }
749 
750  int nNumbers;
751  cout << "How many numbers should we generate: ";
752  cin >> nNumbers; cout << nNumbers << "\n";
753 
754  double k;
755  cout << "Enter k: ";
756  cin >> k; cout << k << "\n";
757 
758  cout << "\nInstantiating distribution utilizing TripleRand engine...\n";
759  TripleRand eng (seed);
760  RandSkewNormal dist (eng, k);
761 
762  cout << "\n Sample fire(): \n";
763  double x;
764 
765  x = dist.fire();
766  cout << x;
767 
768  cout << "\n Testing operator() ... \n";
769 
770  bool good = skewNormalTest ( dist, k, nNumbers );
771 
772  if (good) {
773  return 0;
774  } else {
775  return SkewNormalBAD;
776  }
777 
778 } // testSkewNormal()
779 
780 
781 
782 // ---------
783 // RandGaussT
784 // ---------
785 
787 
788  cout << "\n--------------------------------------------\n";
789  cout << "Test of RandGaussT distribution \n\n";
790 
791  long seed;
792  cout << "Please enter an integer seed: ";
793  cin >> seed; cout << seed << "\n";
794  if (seed == 0) {
795  cout << "Moving on to next test...\n";
796  return 0;
797  }
798 
799  int nNumbers;
800  cout << "How many numbers should we generate: ";
801  cin >> nNumbers; cout << nNumbers << "\n";
802 
803  double mu;
804  double sigma;
805  cout << "Enter mu: ";
806  cin >> mu; cout << mu << "\n";
807 
808  cout << "Enter sigma: ";
809  cin >> sigma; cout << sigma << "\n";
810 
811  cout << "\nInstantiating distribution utilizing TripleRand engine...\n";
812  TripleRand eng (seed);
813  RandGaussT dist (eng, mu, sigma);
814 
815  cout << "\n Sample fire(): \n";
816  double x;
817 
818  x = dist.fire();
819  cout << x;
820 
821  cout << "\n Testing operator() ... \n";
822 
823  bool good = gaussianTest ( dist, mu, sigma, nNumbers );
824 
825  if (good) {
826  return 0;
827  } else {
828  return GaussTBAD;
829  }
830 
831 } // testRandGaussT()
832 
833 
834 
835 // ---------
836 // RandGaussQ
837 // ---------
838 
840 
841  cout << "\n--------------------------------------------\n";
842  cout << "Test of RandGaussQ distribution \n\n";
843 
844  long seed;
845  cout << "Please enter an integer seed: ";
846  cin >> seed; cout << seed << "\n";
847  if (seed == 0) {
848  cout << "Moving on to next test...\n";
849  return 0;
850  }
851 
852  int nNumbers;
853  cout << "How many numbers should we generate: ";
854  cin >> nNumbers; cout << nNumbers << "\n";
855 
856  if (nNumbers >= 20000000) {
857  cout << "With that many samples RandGaussQ need not pass validation...\n";
858  }
859 
860  double mu;
861  double sigma;
862  cout << "Enter mu: ";
863  cin >> mu; cout << mu << "\n";
864 
865  cout << "Enter sigma: ";
866  cin >> sigma; cout << sigma << "\n";
867 
868  cout << "\nInstantiating distribution utilizing DualRand engine...\n";
869  DualRand eng (seed);
870  RandGaussQ dist (eng, mu, sigma);
871 
872  cout << "\n Sample fire(): \n";
873  double x;
874 
875  x = dist.fire();
876  cout << x;
877 
878  cout << "\n Testing operator() ... \n";
879 
880  bool good = gaussianTest ( dist, mu, sigma, nNumbers );
881 
882  if (good) {
883  return 0;
884  } else {
885  return GaussQBAD;
886  }
887 
888 } // testRandGaussQ()
889 
890 
891 // ---------
892 // RandPoisson
893 // ---------
894 
896 
897  cout << "\n--------------------------------------------\n";
898  cout << "Test of RandPoisson distribution \n\n";
899 
900  long seed;
901  cout << "Please enter an integer seed: ";
902  cin >> seed; cout << seed << "\n";
903  if (seed == 0) {
904  cout << "Moving on to next test...\n";
905  return 0;
906  }
907 
908  cout << "\nInstantiating distribution utilizing TripleRand engine...\n";
909  TripleRand eng (seed);
910 
911  int nNumbers;
912  cout << "How many numbers should we generate for each mu: ";
913  cin >> nNumbers; cout << nNumbers << "\n";
914 
915  bool good = true;
916 
917  while (true) {
918  double mu;
919  cout << "Enter a value for mu: ";
920  cin >> mu; cout << mu << "\n";
921  if (mu == 0) break;
922 
923  RandPoisson dist (eng, mu);
924 
925  cout << "\n Sample fire(): \n";
926  double x;
927 
928  x = dist.fire();
929  cout << x;
930 
931  cout << "\n Testing operator() ... \n";
932 
933  bool this_good = poissonTest ( dist, mu, nNumbers );
934  if (!this_good) {
935  cout << "\n Poisson distribution for mu = " << mu << " is incorrect!!!\n";
936  }
937  good &= this_good;
938  } // end of the while(true)
939 
940  if (good) {
941  return 0;
942  } else {
943  return PoissonBAD;
944  }
945 
946 } // testRandPoisson()
947 
948 
949 // ---------
950 // RandPoissonQ
951 // ---------
952 
954 
955  cout << "\n--------------------------------------------\n";
956  cout << "Test of RandPoissonQ distribution \n\n";
957 
958  long seed;
959  cout << "Please enter an integer seed: ";
960  cin >> seed; cout << seed << "\n";
961  if (seed == 0) {
962  cout << "Moving on to next test...\n";
963  return 0;
964  }
965 
966  cout << "\nInstantiating distribution utilizing TripleRand engine...\n";
967  TripleRand eng (seed);
968 
969  int nNumbers;
970  cout << "How many numbers should we generate for each mu: ";
971  cin >> nNumbers; cout << nNumbers << "\n";
972 
973  bool good = true;
974 
975  while (true) {
976  double mu;
977  cout << "Enter a value for mu: ";
978  cin >> mu; cout << mu << "\n";
979  if (mu == 0) break;
980 
981  RandPoissonQ dist (eng, mu);
982 
983  cout << "\n Sample fire(): \n";
984  double x;
985 
986  x = dist.fire();
987  cout << x;
988 
989  cout << "\n Testing operator() ... \n";
990 
991  bool this_good = poissonTest ( dist, mu, nNumbers );
992  if (!this_good) {
993  cout << "\n Poisson distribution for mu = " << mu << " is incorrect!!!\n";
994  }
995  good &= this_good;
996  } // end of the while(true)
997 
998  if (good) {
999  return 0;
1000  } else {
1001  return PoissonQBAD;
1002  }
1003 
1004 } // testRandPoissonQ()
1005 
1006 
1007 // ---------
1008 // RandPoissonT
1009 // ---------
1010 
1012 
1013  cout << "\n--------------------------------------------\n";
1014  cout << "Test of RandPoissonT distribution \n\n";
1015 
1016  long seed;
1017  cout << "Please enter an integer seed: ";
1018  cin >> seed; cout << seed << "\n";
1019  if (seed == 0) {
1020  cout << "Moving on to next test...\n";
1021  return 0;
1022  }
1023 
1024  cout << "\nInstantiating distribution utilizing TripleRand engine...\n";
1025  TripleRand eng (seed);
1026 
1027  int nNumbers;
1028  cout << "How many numbers should we generate for each mu: ";
1029  cin >> nNumbers; cout << nNumbers << "\n";
1030 
1031  bool good = true;
1032 
1033  while (true) {
1034  double mu;
1035  cout << "Enter a value for mu: ";
1036  cin >> mu; cout << mu << "\n";
1037  if (mu == 0) break;
1038 
1039  RandPoissonT dist (eng, mu);
1040 
1041  cout << "\n Sample fire(): \n";
1042  double x;
1043 
1044  x = dist.fire();
1045  cout << x;
1046 
1047  cout << "\n Testing operator() ... \n";
1048 
1049  bool this_good = poissonTest ( dist, mu, nNumbers );
1050  if (!this_good) {
1051  cout << "\n Poisson distribution for mu = " << mu << " is incorrect!!!\n";
1052  }
1053  good &= this_good;
1054  } // end of the while(true)
1055 
1056  if (good) {
1057  return 0;
1058  } else {
1059  return PoissonTBAD;
1060  }
1061 
1062 } // testRandPoissonT()
1063 
1064 
1065 // -----------
1066 // RandGeneral
1067 // -----------
1068 
1070 
1071  cout << "\n--------------------------------------------\n";
1072  cout << "Test of RandGeneral distribution (using a Gaussian shape)\n\n";
1073 
1074  bool good;
1075 
1076  long seed;
1077  cout << "Please enter an integer seed: ";
1078  cin >> seed; cout << seed << "\n";
1079  if (seed == 0) {
1080  cout << "Moving on to next test...\n";
1081  return 0;
1082  }
1083 
1084  int nNumbers;
1085  cout << "How many numbers should we generate: ";
1086  cin >> nNumbers; cout << nNumbers << "\n";
1087 
1088  double mu;
1089  double sigma;
1090  mu = .5; // Since randGeneral always ranges from 0 to 1
1091  sigma = .06;
1092 
1093  cout << "Enter sigma: ";
1094  cin >> sigma; cout << sigma << "\n";
1095  // We suggest sigma be .06. This leaves room for 8 sigma
1096  // in the distribution. If it is much smaller, the number
1097  // of bins necessary to expect a good match will increase.
1098  // If sigma is much larger, the cutoff before 5 sigma can
1099  // cause the Gaussian hypothesis to be rejected. At .14, for
1100  // example, the 4th moment is 7 sigma away from expectation.
1101 
1102  int nBins;
1103  cout << "Enter nBins for stepwise pdf test: ";
1104  cin >> nBins; cout << nBins << "\n";
1105  // We suggest at least 10000 bins; fewer would risk
1106  // false rejection because the step-function curve
1107  // does not match an actual Gaussian. At 10000 bins,
1108  // a million-hit test does not have the resolving power
1109  // to tell the boxy pdf from the true Gaussian. At 5000
1110  // bins, it does.
1111 
1112  double xBins = nBins;
1113  double* aProbFunc = new double [nBins];
1114  double x;
1115  for ( int iBin = 0; iBin < nBins; iBin++ ) {
1116  x = iBin / (xBins-1);
1117  aProbFunc [iBin] = std::exp ( - (x-mu)*(x-mu) / (2*sigma*sigma) );
1118  }
1119  // Note that this pdf is not normalized; RandGeneral does that
1120 
1121  cout << "\nInstantiating distribution utilizing Ranlux64 engine...\n";
1122  Ranlux64Engine eng (seed, 3);
1123 
1124  { // Open block for testing type 1 - step function pdf
1125 
1126  RandGeneral dist (eng, aProbFunc, nBins, 1);
1127  delete[] aProbFunc;
1128 
1129  double* garbage = new double[nBins];
1130  // We wish to verify that deleting the pdf
1131  // after instantiating the engine is fine.
1132  for ( int gBin = 0; gBin < nBins; gBin++ ) {
1133  garbage [gBin] = 1;
1134  }
1135 
1136  cout << "\n Sample fire(): \n";
1137 
1138  x = dist.fire();
1139  cout << x;
1140 
1141  cout << "\n Testing operator() ... \n";
1142 
1143  good = gaussianTest ( dist, mu, sigma, nNumbers );
1144 
1145  delete[] garbage;
1146 
1147  } // Close block for testing type 1 - step function pdf
1148  // dist goes out of scope but eng is supposed to stick around;
1149  // by closing this block we shall verify that!
1150 
1151  cout << "Enter nBins for linearized pdf test: ";
1152  cin >> nBins; cout << nBins << "\n";
1153  // We suggest at least 1000 bins; fewer would risk
1154  // false rejection because the non-smooth curve
1155  // does not match an actual Gaussian. At 1000 bins,
1156  // a million-hit test does not resolve the non-smoothness;
1157  // at 300 bins it does.
1158 
1159  xBins = nBins;
1160  aProbFunc = new double [nBins];
1161  for ( int jBin = 0; jBin < nBins; jBin++ ) {
1162  x = jBin / (xBins-1);
1163  aProbFunc [jBin] = std::exp ( - (x-mu)*(x-mu) / (2*sigma*sigma) );
1164  }
1165  // Note that this pdf is not normalized; RandGeneral does that
1166 
1167  RandGeneral dist (eng, aProbFunc, nBins, 0);
1168 
1169  cout << "\n Sample operator(): \n";
1170 
1171  x = dist();
1172  cout << x;
1173 
1174  cout << "\n Testing operator() ... \n";
1175 
1176  bool good2 = gaussianTest ( dist, mu, sigma, nNumbers );
1177  good = good && good2;
1178 
1179  if (good) {
1180  return 0;
1181  } else {
1182  return GeneralBAD;
1183  }
1184 
1185 } // testRandGeneral()
1186 
1187 
1188 
1189 
1190 // **********************
1191 //
1192 // SECTION IV - Main
1193 //
1194 // **********************
1195 
1196 int main() {
1197 
1198  int mask = 0;
1199 
1200  mask |= testRandGauss();
1201  mask |= testRandGaussQ();
1202  mask |= testRandGaussT();
1203 
1204  mask |= testRandGeneral();
1205 
1206  mask |= testRandPoisson();
1207  mask |= testRandPoissonQ();
1208  mask |= testRandPoissonT();
1209 
1210  mask |= testSkewNormal(); // k = 0 (gaussian)
1211  mask |= testSkewNormal(); // k = -2
1212  mask |= testSkewNormal(); // k = 1
1213  mask |= testSkewNormal(); // k = 5
1214 
1215  return mask > 0 ? -mask : mask;
1216 }
1217 
bool skewNormalTest(HepRandom &dist, double k, int nNumbers)
double operator()(int r)
int testRandGaussQ()
int testRandPoissonT()
bool poissonTest(RandPoisson &dist, double mu, int N)
int testRandPoisson()
int testRandGeneral()
double gammln(double xx)
Definition: RandPoisson.cc:55
double * createRefDist(poisson pdist, int N, int MINBIN, int MAXBINS, int clumping, int &firstBin, int &lastBin)
#define exit(x)
int testRandGauss()
int testRandPoissonQ()
bool gaussianTest(HepRandom &dist, double mu, double sigma, int nNumbers)
int testSkewNormal()
int testRandGaussT()
poisson(double mu)
int main()