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
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
53using std::cin;
54using std::cout;
55using std::cerr;
56using std::endl;
57using std::setprecision;
58using namespace CLHEP;
59//#ifndef _WIN32
60//using std::exp;
61//#endif
62
63// Tolerance of deviation from expected results
64static const double REJECT = 4.0;
65
66// Mask bits to form a word indicating which if any dists were "bad"
67static const int GaussBAD = 1 << 0;
68static const int GeneralBAD = 1 << 1;
69static const int PoissonBAD = 1 << 2;
70static const int GaussQBAD = 1 << 3;
71static const int GaussTBAD = 1 << 4;
72static const int PoissonQBAD = 1 << 5;
73static const int PoissonTBAD = 1 << 6;
74static const int SkewNormalBAD = 1 << 7;
75
76
77// **********************
78//
79// SECTION I - General tools for the various tests
80//
81// **********************
82
83static
84double 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
116static
117double 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
135static
136double 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
162static
163double 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
182bool 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
350bool 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
451class 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
461double* 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
529bool 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
1196int 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
poisson(double mu)
double operator()(int r)
#define exit(x)
double gammln(double xx)
int testRandGaussQ()
int testRandPoissonT()
int testRandPoissonQ()
bool gaussianTest(HepRandom &dist, double mu, double sigma, int nNumbers)
bool poissonTest(RandPoisson &dist, double mu, int N)
int testRandPoisson()
int testRandGauss()
double * createRefDist(poisson pdist, int N, int MINBIN, int MAXBINS, int clumping, int &firstBin, int &lastBin)
int testRandGeneral()
bool skewNormalTest(HepRandom &dist, double k, int nNumbers)
int testSkewNormal()
int testRandGaussT()
int main()