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 
delta
HepRotation delta() setPhi()
CLHEP::RandPoissonT
Definition: Matrix/CLHEP/Random/RandPoissonT.h:41
poissonTest
bool poissonTest(RandPoisson &dist, double mu, int N)
Definition: testRandDists.cc:529
testRandPoissonT
int testRandPoissonT()
Definition: testRandDists.cc:1011
CLHEP::RandPoissonQ::fire
long fire()
Definition: RandPoissonQ.cc:134
a
@ a
Definition: testCategories.cc:125
CLHEP::RandGaussQ::fire
double fire()
CLHEP::RandSkewNormal
Definition: Matrix/CLHEP/Random/RandSkewNormal.h:34
samples
this we validated samples(it is more time consuming) with no problems. We validated the quick() routine(which differs only above mu
b
@ b
Definition: testCategories.cc:125
CLHEP::RandPoisson::fire
long fire()
Definition: RandPoisson.cc:214
CLHEP::TripleRand
Definition: Matrix/CLHEP/Random/TripleRand.h:52
mu
the goal is to keep the overall false rejection probability down at the to level For each validated we discuss which of course is by necessity relative timing We take the time for a single random via one of the fastest good and at any rate the ratios will vary by around depending on the processor and memory configuration used A timing for a distribution of units would mean no time used beyond the uniform random Summary Distribution Validated Validation Rejected Past N RandGauss RandGaussT RandGaussQ bins stepwise bins RandPoisson RandPoissonT mu< 100 N=50, 000, 000 ------- mu > RandPoissonQ mu< 100 N=50, 000, 000 -------(same as RandPoissonT) mu > RandGauss shoot() etc 2.5 units Validation tests applied and a very accurate series is substituted for the table method shoot() etc 1.7 units Validation tests applied and below about **a series approximation is but we have applied this test with N up to and never get anywhere near the rejection point Analytical considerations indicate that the validation test would not reject until O(10 **24) samples were inspected. ---------------------------------------------------------- 2. RandGeneral Method since we wish to have good resolution power even if just one of the mu values is it would be unwise to dial this down any further Validation for several values o mu)
Definition: validation.doc:263
CLHEP::RandGaussQ
Definition: Matrix/CLHEP/Random/RandGaussQ.h:32
CLHEP::RandGeneral
Definition: Matrix/CLHEP/Random/RandGeneral.h:40
poisson::operator()
double operator()(int r)
Definition: testRandDists.cc:455
CLHEP::gammln
double gammln(double xx)
Definition: RandPoisson.cc:55
testRandPoisson
int testRandPoisson()
Definition: testRandDists.cc:895
createRefDist
double * createRefDist(poisson pdist, int N, int MINBIN, int MAXBINS, int clumping, int &firstBin, int &lastBin)
Definition: testRandDists.cc:461
CLHEP::detail::n
n
Definition: Ranlux64Engine.cc:85
CLHEP::RandGaussT::fire
double fire()
CLHEP::RandSkewNormal::fire
double fire()
Definition: RandSkewNormal.cc:80
testRandPoissonQ
int testRandPoissonQ()
Definition: testRandDists.cc:953
N
the goal is to keep the overall false rejection probability down at the to level For each validated we discuss which of course is by necessity relative timing We take the time for a single random via one of the fastest good and at any rate the ratios will vary by around depending on the processor and memory configuration used A timing for a distribution of units would mean no time used beyond the uniform random Summary Distribution Validated Validation Rejected Past N RandGauss N
Definition: validation.doc:48
testRandGeneral
int testRandGeneral()
Definition: testRandDists.cc:1069
CLHEP
Definition: ClhepVersion.h:13
testRandGauss
int testRandGauss()
Definition: testRandDists.cc:684
CLHEP::RandGaussT
Definition: Matrix/CLHEP/Random/RandGaussT.h:41
poisson
Definition: testRandDists.cc:451
CLHEP::Ranlux64Engine
Definition: Matrix/CLHEP/Random/Ranlux64Engine.h:50
main
int main()
Definition: testRandDists.cc:1196
poisson::poisson
poisson(double mu)
Definition: testRandDists.cc:454
CLHEP::DualRand
Definition: Matrix/CLHEP/Random/DualRand.h:51
i
long i
Definition: JamesRandomSeeding.txt:27
testRandGaussT
int testRandGaussT()
Definition: testRandDists.cc:786
gaussianTest
bool gaussianTest(HepRandom &dist, double mu, double sigma, int nNumbers)
Definition: testRandDists.cc:182
CLHEP::RandGauss
Definition: Matrix/CLHEP/Random/RandGauss.h:42
CLHEP::RandGeneral::fire
double fire()
CLHEP::RandPoissonT::fire
long fire()
Definition: RandPoissonT.cc:73
testSkewNormal
int testSkewNormal()
Definition: testRandDists.cc:737
CLHEP::RandPoissonQ
Definition: Matrix/CLHEP/Random/RandPoissonQ.h:33
testRandGaussQ
int testRandGaussQ()
Definition: testRandDists.cc:839
exit
it has advantages For I leave the ZMthrows but substitute I replaced ZMthrow with ZMthrowA in this package when will issue its message and exit(-1). Also
CLHEP::RandPoisson
Definition: Matrix/CLHEP/Random/RandPoisson.h:42
x
any side effects of that construction would occur twice The semantics of throw x
Definition: whyZMthrowRethrows.txt:37
k
long k
Definition: JamesRandomSeeding.txt:29
CLHEP::HepRandom
Definition: Matrix/CLHEP/Random/Random.h:50
CLHEP::RandGauss::fire
double fire()
skewNormalTest
bool skewNormalTest(HepRandom &dist, double k, int nNumbers)
Definition: testRandDists.cc:350