CLHEP 2.4.6.4
C++ Class Library for High Energy Physics
Loading...
Searching...
No Matches
testRandDists.cc File Reference
#include "CLHEP/Units/GlobalPhysicalConstants.h"
#include "CLHEP/Random/Randomize.h"
#include "CLHEP/Random/RandGaussQ.h"
#include "CLHEP/Random/RandGaussT.h"
#include "CLHEP/Random/RandPoissonQ.h"
#include "CLHEP/Random/RandPoissonT.h"
#include "CLHEP/Random/RandSkewNormal.h"
#include "CLHEP/Random/defs.h"
#include <iostream>
#include <iomanip>
#include <cmath>
#include <stdlib.h>
#include <cstdlib>

Go to the source code of this file.

Classes

class  poisson
 

Functions

bool gaussianTest (HepRandom &dist, double mu, double sigma, int nNumbers)
 
bool skewNormalTest (HepRandom &dist, double k, int nNumbers)
 
doublecreateRefDist (poisson pdist, int N1, int MINBIN, int MAXBINS, int clumping, int &firstBin, int &lastBin)
 
bool poissonTest (RandPoisson &dist, double mu, int N2)
 
int testRandGauss ()
 
int testSkewNormal ()
 
int testRandGaussT ()
 
int testRandGaussQ ()
 
int testRandPoisson ()
 
int testRandPoissonQ ()
 
int testRandPoissonT ()
 
int testRandGeneral ()
 
int main ()
 

Function Documentation

◆ createRefDist()

double * createRefDist ( poisson  pdist,
int  N1,
int  MINBIN,
int  MAXBINS,
int  clumping,
int &  firstBin,
int &  lastBin 
)

Definition at line 459 of file testRandDists.cc.

461 {
462
463 // Create the reference distribution -- making sure there are more than
464 // 20 points at each value. The entire tail will be rolled up into one
465 // value (at each end). We shall end up with some range of bins starting
466 // at 0 or more, and ending at MAXBINS-1 or less.
467
468 double * refdist = new double [MAXBINS];
469
470 int c = 0; // c is the number of the clump, that is, the member number
471 // of the refdist array.
472 int ic = 0; // ic is the number within the clump; mod clumping
473 int r = 0; // r is the value of the variate.
474
475 // Determine the first bin: at least 20 entries must be at the level
476 // of that bin (so that we won't immediately dip belpw 20) but the number
477 // to enter is cumulative up to that bin.
478
479 double start = 0;
480 double binc;
481 while ( c < MAXBINS ) {
482 for ( ic=0, binc=0; ic < clumping; ic++, r++ ) {
483 binc += pdist(r) * N1;
484 }
485 start += binc;
486 if (binc >= MINBIN) break;
487 c++;
488 if ( c > MAXBINS/3 ) {
489 cout << "The number of samples supplied " << N1 <<
490 " is too small to set up a chi^2 to test this distribution.\n";
491 exit(-1);
492 }
493 }
494 firstBin = c;
495 refdist[firstBin] = start;
496 c++;
497
498 // Fill all the other bins until one has less than 20 items.
499 double next = 0;
500 while ( c < MAXBINS ) {
501 for ( ic=0, binc=0; ic < clumping; ic++, r++ ) {
502 binc += pdist(r) * N1;
503 }
504 next = binc;
505 if (next < MINBIN) break;
506 refdist[c] = next;
507 c++;
508 }
509
510 // Shove all the remaining items into last bin.
511 lastBin = c-1;
512 next += refdist[lastBin];
513 while ( c < MAXBINS ) {
514 for ( ic=0, binc=0; ic < clumping; ic++, r++ ) {
515 binc += pdist(r) * N1;
516 }
517 next += binc;
518 c++;
519 }
520 refdist[lastBin] = next;
521
522 return refdist;
523
524} // createRefDist()
#define exit(x)

Referenced by poissonTest().

◆ gaussianTest()

bool gaussianTest ( HepRandom dist,
double  mu,
double  sigma,
int  nNumbers 
)

Definition at line 180 of file testRandDists.cc.

181 {
182
183 bool good = true;
184 double worstSigma = 0;
185
186// We will accumulate mean and moments up to the sixth,
187// The second moment should be sigma**2, the fourth 3 sigma**4,
188// the sixth 15 sigma**6. The expected variance in these is
189// (for the m-th moment with m even) (2m-1)!! (m-1)!!**2 / n
190// (for the m-th moment with m odd) (2m-1)!! m!!**2 / n
191// We also do a histogram with bins every half sigma.
192
193 double sumx = 0;
194 double sumx2 = 0;
195 double sumx3 = 0;
196 double sumx4 = 0;
197 double sumx5 = 0;
198 double sumx6 = 0;
199 int counts[11];
200 int ncounts[11];
201 int ciu;
202 for (ciu = 0; ciu < 11; ciu++) {
203 counts[ciu] = 0;
204 ncounts[ciu] = 0;
205 }
206
207 long oldprecision = cout.precision();
208 cout.precision(5);
209 // hack so that gcc 4.3 puts x and u into memory instead of a register
210 volatile double x;
211 volatile double u;
212 int ipr = nNumbers / 10 + 1;
213 for (int ifire = 0; ifire < nNumbers; ifire++) {
214 x = dist(); // We avoid fire() because that is not virtual
215 // in HepRandom.
216 if( x < mu - 12.0*sigma ) {
217 cout << "x = " << x << "\n";
218 }
219 if ( (ifire % ipr) == 0 ) {
220 cout << ifire << endl;
221 }
222 sumx += x;
223 sumx2 += x*x;
224 sumx3 += x*x*x;
225 sumx4 += x*x*x*x;
226 sumx5 += x*x*x*x*x;
227 sumx6 += x*x*x*x*x*x;
228 u = (x - mu) / sigma;
229 if ( u >= 0 ) {
230 ciu = (int)(2*u);
231 if (ciu>10) ciu = 10;
232 counts[ciu]++;
233 } else {
234 ciu = (int)(2*(-u));
235 if (ciu>10) ciu = 10;
236 ncounts[ciu]++;
237 }
238 }
239
240 double mean = sumx / nNumbers;
241 double u2 = sumx2/nNumbers - mean*mean;
242 double u3 = sumx3/nNumbers - 3*sumx2*mean/nNumbers + 2*mean*mean*mean;
243 double u4 = sumx4/nNumbers - 4*sumx3*mean/nNumbers
244 + 6*sumx2*mean*mean/nNumbers - 3*mean*mean*mean*mean;
245 double u5 = sumx5/nNumbers - 5*sumx4*mean/nNumbers
246 + 10*sumx3*mean*mean/nNumbers
247 - 10*sumx2*mean*mean*mean/nNumbers
248 + 4*mean*mean*mean*mean*mean;
249 double u6 = sumx6/nNumbers - 6*sumx5*mean/nNumbers
250 + 15*sumx4*mean*mean/nNumbers
251 - 20*sumx3*mean*mean*mean/nNumbers
252 + 15*sumx2*mean*mean*mean*mean/nNumbers
253 - 5*mean*mean*mean*mean*mean*mean;
254
255 cout << "Mean (should be close to " << mu << "): " << mean << endl;
256 cout << "Second moment (should be close to " << sigma*sigma <<
257 "): " << u2 << endl;
258 cout << "Third moment (should be close to zero): " << u3 << endl;
259 cout << "Fourth moment (should be close to " << 3*sigma*sigma*sigma*sigma <<
260 "): " << u4 << endl;
261 cout << "Fifth moment (should be close to zero): " << u5 << endl;
262 cout << "Sixth moment (should be close to "
263 << 15*sigma*sigma*sigma*sigma*sigma*sigma
264 << "): " << u6 << endl;
265
266 // For large N, the variance squared in the scaled 2nd, 3rd 4th 5th and
267 // 6th moments are roughly 2/N, 6/N, 96/N, 720/N and 10170/N respectively.
268 // Based on this, we can judge how many sigma a result represents:
269
270 double del1 = std::sqrt ( (double) nNumbers ) * std::abs(mean - mu) / sigma;
271 double del2 = std::sqrt ( nNumbers/2.0 ) * std::abs(u2 - sigma*sigma) / (sigma*sigma);
272 double del3 = std::sqrt ( nNumbers/6.0 ) * std::abs(u3) / (sigma*sigma*sigma);
273 double sigma4 = sigma*sigma*sigma*sigma;
274 double del4 = std::sqrt ( nNumbers/96.0 ) * std::abs(u4 - 3 * sigma4) / sigma4;
275 double del5 = std::sqrt ( nNumbers/720.0 ) * std::abs(u5) / (sigma*sigma4);
276 double del6 = std::sqrt ( nNumbers/10170.0 ) * std::abs(u6 - 15*sigma4*sigma*sigma)
277 / (sigma4*sigma*sigma);
278
279 cout << " These represent " <<
280 del1 << ", " << del2 << ", " << del3 << ", \n"
281 <<" " << del4 << ", " << del5 << ", " << del6
282 <<"\n standard deviations from expectations\n";
283 if ( del1 > worstSigma ) worstSigma = del1;
284 if ( del2 > worstSigma ) worstSigma = del2;
285 if ( del3 > worstSigma ) worstSigma = del3;
286 if ( del4 > worstSigma ) worstSigma = del4;
287 if ( del5 > worstSigma ) worstSigma = del5;
288 if ( del6 > worstSigma ) worstSigma = del6;
289
290 if ( del1 > REJECT || del2 > REJECT || del3 > REJECT
291 || del4 > REJECT || del5 > REJECT || del6 > REJECT ) {
292 cout << "REJECT hypothesis that this distribution is correct!!\n";
293 good = false;
294 }
295
296 // The variance of the bin counts is given by a Poisson estimate (std::sqrt(npq)).
297
298 double table[11] = { // Table of integrated density in each range:
299 .191462, // 0.0 - 0.5 sigma
300 .149882, // 0.5 - 1.0 sigma
301 .091848, // 1.0 - 1.5 sigma
302 .044057, // 1.5 - 2.0 sigma
303 .016540, // 2.0 - 2.5 sigma
304 .004860, // 2.5 - 3.0 sigma
305 .001117, // 3.0 - 3.5 sigma
306 .000201, // 3.5 - 4.0 sigma
307 2.83E-5, // 4.0 - 4.5 sigma
308 3.11E-6, // 4.5 - 5.0 sigma
309 3.87E-7 // 5.0 sigma and up
310 };
311
312 for (int m1 = 0; m1 < 11; m1++) {
313 double expect = table[m1]*nNumbers;
314 double sig = std::sqrt ( table[m1] * (1.0-table[m1]) * nNumbers );
315 cout.precision(oldprecision);
316 cout << "Between " << m1/2.0 << " sigma and "
317 << m1/2.0+.5 << " sigma (should be about " << expect << "):\n "
318 << " "
319 << ncounts[m1] << " negative and " << counts[m1] << " positive " << "\n";
320 cout.precision(5);
321 double negSigs = std::abs ( ncounts[m1] - expect ) / sig;
322 double posSigs = std::abs ( counts[m1] - expect ) / sig;
323 cout << " These represent " <<
324 negSigs << " and " << posSigs << " sigma from expectations\n";
325 if ( negSigs > REJECT || posSigs > REJECT ) {
326 cout << "REJECT hypothesis that this distribution is correct!!\n";
327 good = false;
328 }
329 if ( negSigs > worstSigma ) worstSigma = negSigs;
330 if ( posSigs > worstSigma ) worstSigma = posSigs;
331 }
332
333 cout << "\n The worst deviation encountered (out of about 25) was "
334 << worstSigma << " sigma \n\n";
335
336 cout.precision(oldprecision);
337
338 return good;
339
340} // gaussianTest()

Referenced by testRandGauss(), testRandGaussQ(), testRandGaussT(), and testRandGeneral().

◆ main()

int main ( )

Definition at line 1194 of file testRandDists.cc.

1194 {
1195
1196 int mask = 0;
1197
1198 mask |= testRandGauss();
1199 mask |= testRandGaussQ();
1200 mask |= testRandGaussT();
1201
1202 mask |= testRandGeneral();
1203
1204 mask |= testRandPoisson();
1205 mask |= testRandPoissonQ();
1206 mask |= testRandPoissonT();
1207
1208 mask |= testSkewNormal(); // k = 0 (gaussian)
1209 mask |= testSkewNormal(); // k = -2
1210 mask |= testSkewNormal(); // k = 1
1211 mask |= testSkewNormal(); // k = 5
1212
1213 return mask > 0 ? -mask : mask;
1214}
int testRandGaussQ()
int testRandPoissonT()
int testRandPoissonQ()
int testRandPoisson()
int testRandGauss()
int testRandGeneral()
int testSkewNormal()
int testRandGaussT()

◆ poissonTest()

bool poissonTest ( RandPoisson dist,
double  mu,
int  N2 
)

Definition at line 527 of file testRandDists.cc.

527 {
528
529// Three tests will be done:
530//
531// A chi-squared test will be used to test the hypothesis that the
532// generated distribution of N2 numbers matches the proper Poisson distribution.
533//
534// The same test will be applied to the distribution of numbers "clumping"
535// together std::sqrt(mu) bins. This will detect small deviations over several
536// touching bins, when mu is not small.
537//
538// The mean and second moment are checked against their theoretical values.
539
540 bool good = true;
541
542 int clumping = int(std::sqrt(mu));
543 if (clumping <= 1) clumping = 2;
544 const int MINBIN = 20;
545 const int MAXBINS = 1000;
546 int firstBin;
547 int lastBin;
548 int firstBin2;
549 int lastBin2;
550
551 poisson pdist(mu);
552
553 double* refdist = createRefDist( pdist, N2,
554 MINBIN, MAXBINS, 1, firstBin, lastBin);
555 double* refdist2 = createRefDist( pdist, N2,
556 MINBIN, MAXBINS, clumping, firstBin2, lastBin2);
557
558 // Now roll the random dists, treating the tails in the same way as we go.
559
560 double sum = 0;
561 double moment = 0;
562
563 double* samples = new double [MAXBINS];
564 double* samples2 = new double [MAXBINS];
565 int r;
566 for (r = 0; r < MAXBINS; r++) {
567 samples[r] = 0;
568 samples2[r] = 0;
569 }
570
571 int r1;
572 int r2;
573 for (int i = 0; i < N2; i++) {
574 r = (int)dist.fire();
575 sum += r;
576 moment += (r - mu)*(r - mu);
577 r1 = r;
578 if (r1 < firstBin) r1 = firstBin;
579 if (r1 > lastBin) r1 = lastBin;
580 samples[r1] += 1;
581 r2 = r/clumping;
582 if (r2 < firstBin2) r2 = firstBin2;
583 if (r2 > lastBin2) r2 = lastBin2;
584 samples2[r2] += 1;
585 }
586
587// #ifdef DIAGNOSTIC
588 int k;
589 for (k = firstBin; k <= lastBin; k++) {
590 cout << k << " " << samples[k] << " " << refdist[k] << " " <<
591 (samples[k]-refdist[k])*(samples[k]-refdist[k])/refdist[k] << "\n";
592 }
593 cout << "----\n";
594 for (k = firstBin2; k <= lastBin2; k++) {
595 cout << k << " " << samples2[k] << " " << refdist2[k] << "\n";
596 }
597// #endif // DIAGNOSTIC
598
599 // Now find chi^2 for samples[] to apply the first test
600
601 double chi2 = 0;
602 for ( r = firstBin; r <= lastBin; r++ ) {
603 double delta = (samples[r] - refdist[r]);
604 chi2 += delta*delta/refdist[r];
605 }
606 int degFreedom = (lastBin - firstBin + 1) - 1;
607
608 // and finally, p. Since we only care about it for small values,
609 // and never care about it past the 10% level, we can use the approximations
610 // CL(chi^2,n) = 1/std::sqrt(CLHEP::twopi) * ErrIntC ( y ) with
611 // y = std::sqrt(2*chi2) - std::sqrt(2*n-1)
612 // errIntC (y) = std::exp((-y^2)/2)/(y*std::sqrt(CLHEP::twopi))
613
614 double pval;
615 pval = 1.0 - gammp ( .5*degFreedom , .5*chi2 );
616
617 cout << "Chi^2 is " << chi2 << " on " << degFreedom << " degrees of freedom."
618 << " p = " << pval << "\n";
619
620 delete[] refdist;
621 delete[] samples;
622
623 // Repeat the chi^2 and p for the clumped sample, to apply the second test
624
625 chi2 = 0;
626 for ( r = firstBin2; r <= lastBin2; r++ ) {
627 double delta = (samples2[r] - refdist2[r]);
628 chi2 += delta*delta/refdist2[r];
629 }
630 degFreedom = (lastBin2 - firstBin2 + 1) - 1;
631 double pval2;
632 pval2 = 1.0 - gammp ( .5*degFreedom , .5*chi2 );
633
634 cout << "Clumps: Chi^2 is " << chi2 << " on " << degFreedom <<
635 " degrees of freedom." << " p = " << pval2 << "\n";
636
637 delete[] refdist2;
638 delete[] samples2;
639
640 // Check out the mean and sigma to apply the third test
641
642 double mean = sum / N2;
643 double sigma = std::sqrt( moment / (N2-1) );
644
645 double deviationMean = std::fabs(mean - mu)/(std::sqrt(mu/N2));
646 double expectedSigma2Variance = (2*N2*mu*mu/(N2-1) + mu) / N2;
647 double deviationSigma = std::fabs(sigma*sigma-mu)/std::sqrt(expectedSigma2Variance);
648
649 cout << "Mean (should be " << mu << ") is " << mean << "\n";
650 cout << "Sigma (should be " << std::sqrt(mu) << ") is " << sigma << "\n";
651
652 cout << "These are " << deviationMean << " and " << deviationSigma <<
653 " standard deviations from expected values\n\n";
654
655 // If either p-value for the chi-squared tests is less that .0001, or
656 // either the mean or sigma are more than 3.5 standard deviations off,
657 // then reject the validation. This would happen by chance one time
658 // in 2000. Since we will be validating for several values of mu, the
659 // net chance of false rejection remains acceptable.
660
661 if ( (pval < .0001) || (pval2 < .0001) ||
662 (deviationMean > 3.5) || (deviationSigma > 3.5) ) {
663 good = false;
664 cout << "REJECT this distributon!!!\n";
665 }
666
667 return good;
668
669} // poissonTest()
double * createRefDist(poisson pdist, int N1, int MINBIN, int MAXBINS, int clumping, int &firstBin, int &lastBin)

Referenced by testRandPoisson(), testRandPoissonQ(), and testRandPoissonT().

◆ skewNormalTest()

bool skewNormalTest ( HepRandom dist,
double  k,
int  nNumbers 
)

Definition at line 348 of file testRandDists.cc.

348 {
349
350 bool good = true;
351 double worstSigma = 0;
352
353// We will accumulate mean and moments up to the sixth,
354// The second moment should be sigma**2, the fourth 3 sigma**4.
355// The expected variance in these is
356// (for the m-th moment with m even) (2m-1)!! (m-1)!!**2 / n
357// (for the m-th moment with m odd) (2m-1)!! m!!**2 / n
358
359 double sumx = 0;
360 double sumx2 = 0;
361 double sumx3 = 0;
362 double sumx4 = 0;
363 double sumx5 = 0;
364 double sumx6 = 0;
365
366 long oldprecision = cout.precision();
367 cout.precision(5);
368 // hack so that gcc 4.3 puts x into memory instead of a register
369 volatile double x;
370 // calculate mean and sigma
371 double delta = k / std::sqrt( 1 + k*k );
372 double mu = delta/std::sqrt(CLHEP::halfpi);
373 double mom2 = 1.;
374 double mom3 = 3*delta*(1-(delta*delta)/3.)/std::sqrt(CLHEP::halfpi);
375 double mom4 = 3.;
376 double mom5 = 15*delta*(1-2.*(delta*delta)/3.+(delta*delta*delta*delta)/5.)/std::sqrt(CLHEP::halfpi);
377 double mom6 = 15.;
378
379 int ipr = nNumbers / 10 + 1;
380 for (int ifire = 0; ifire < nNumbers; ifire++) {
381 x = dist(); // We avoid fire() because that is not virtual
382 // in HepRandom.
383 if( x < mu - 12.0 ) {
384 cout << "x = " << x << "\n";
385 }
386 if ( (ifire % ipr) == 0 ) {
387 cout << ifire << endl;
388 }
389 sumx += x;
390 sumx2 += x*x;
391 sumx3 += x*x*x;
392 sumx4 += x*x*x*x;
393 sumx5 += x*x*x*x*x;
394 sumx6 += x*x*x*x*x*x;
395 }
396
397 double mean = sumx / nNumbers;
398 double u2 = sumx2/nNumbers;
399 double u3 = sumx3/nNumbers;
400 double u4 = sumx4/nNumbers;
401 double u5 = sumx5/nNumbers;
402 double u6 = sumx6/nNumbers;
403
404 cout << "Mean (should be close to " << mu << "): " << mean << endl;
405 cout << "Second moment (should be close to " << mom2 << "): " << u2 << endl;
406 cout << "Third moment (should be close to " << mom3 << "): " << u3 << endl;
407 cout << "Fourth moment (should be close to " << mom4 << "): " << u4 << endl;
408 cout << "Fifth moment (should be close to " << mom5 << "): " << u5 << endl;
409 cout << "Sixth moment (should be close to " << mom6 << "): " << u6 << endl;
410
411 double del1 = std::sqrt ( (double) nNumbers ) * std::abs(mean - mu);
412 double del2 = std::sqrt ( nNumbers/2.0 ) * std::abs(u2 - mom2);
413 double del3 = std::sqrt ( nNumbers/(15.-mom3*mom3) ) * std::abs(u3 - mom3 );
414 double del4 = std::sqrt ( nNumbers/96.0 ) * std::abs(u4 - mom4);
415 double del5 = std::sqrt ( nNumbers/(945.-mom5*mom5) ) * std::abs(u5 - mom5 );
416 double del6 = std::sqrt ( nNumbers/10170.0 ) * std::abs(u6 - mom6);
417
418 cout << " These represent " <<
419 del1 << ", " << del2 << ", " << del3 << ", \n"
420 <<" " << del4 << ", " << del5 << ", " << del6
421 <<"\n standard deviations from expectations\n";
422 if ( del1 > worstSigma ) worstSigma = del1;
423 if ( del2 > worstSigma ) worstSigma = del2;
424 if ( del3 > worstSigma ) worstSigma = del3;
425 if ( del4 > worstSigma ) worstSigma = del4;
426 if ( del5 > worstSigma ) worstSigma = del5;
427 if ( del6 > worstSigma ) worstSigma = del6;
428
429 if ( del1 > REJECT || del2 > REJECT || del3 > REJECT ||
430 del4 > REJECT || del5 > REJECT || del6 > REJECT ) {
431 cout << "REJECT hypothesis that this distribution is correct!!\n";
432 good = false;
433 }
434
435 cout << "\n The worst deviation encountered (out of about 25) was "
436 << worstSigma << " sigma \n\n";
437
438 cout.precision(oldprecision);
439
440 return good;
441
442} // skewNormalTest()

Referenced by testSkewNormal().

◆ testRandGauss()

int testRandGauss ( )

Definition at line 682 of file testRandDists.cc.

682 {
683
684 cout << "\n--------------------------------------------\n";
685 cout << "Test of RandGauss distribution \n\n";
686
687 long seed;
688 cout << "Please enter an integer seed: ";
689 cin >> seed; cout << seed << "\n";
690 if (seed == 0) {
691 cout << "Moving on to next test...\n";
692 return 0;
693 }
694
695 int nNumbers;
696 cout << "How many numbers should we generate: ";
697 cin >> nNumbers; cout << nNumbers << "\n";
698
699 double mu;
700 double sigma;
701 cout << "Enter mu: ";
702 cin >> mu; cout << mu << "\n";
703
704 cout << "Enter sigma: ";
705 cin >> sigma; cout << sigma << "\n";
706
707 cout << "\nInstantiating distribution utilizing TripleRand engine...\n";
708 TripleRand eng (seed);
709 RandGauss dist (eng, mu, sigma);
710
711 cout << "\n Sample fire(): \n";
712 double x;
713
714 x = dist.fire();
715 cout << x;
716
717 cout << "\n Testing operator() ... \n";
718
719 bool good = gaussianTest ( dist, mu, sigma, nNumbers );
720
721 if (good) {
722 return 0;
723 } else {
724 return GaussBAD;
725 }
726
727} // testRandGauss()
bool gaussianTest(HepRandom &dist, double mu, double sigma, int nNumbers)

Referenced by main().

◆ testRandGaussQ()

int testRandGaussQ ( )

Definition at line 837 of file testRandDists.cc.

837 {
838
839 cout << "\n--------------------------------------------\n";
840 cout << "Test of RandGaussQ distribution \n\n";
841
842 long seed;
843 cout << "Please enter an integer seed: ";
844 cin >> seed; cout << seed << "\n";
845 if (seed == 0) {
846 cout << "Moving on to next test...\n";
847 return 0;
848 }
849
850 int nNumbers;
851 cout << "How many numbers should we generate: ";
852 cin >> nNumbers; cout << nNumbers << "\n";
853
854 if (nNumbers >= 20000000) {
855 cout << "With that many samples RandGaussQ need not pass validation...\n";
856 }
857
858 double mu;
859 double sigma;
860 cout << "Enter mu: ";
861 cin >> mu; cout << mu << "\n";
862
863 cout << "Enter sigma: ";
864 cin >> sigma; cout << sigma << "\n";
865
866 cout << "\nInstantiating distribution utilizing DualRand engine...\n";
867 DualRand eng (seed);
868 RandGaussQ dist (eng, mu, sigma);
869
870 cout << "\n Sample fire(): \n";
871 double x;
872
873 x = dist.fire();
874 cout << x;
875
876 cout << "\n Testing operator() ... \n";
877
878 bool good = gaussianTest ( dist, mu, sigma, nNumbers );
879
880 if (good) {
881 return 0;
882 } else {
883 return GaussQBAD;
884 }
885
886} // testRandGaussQ()

Referenced by main().

◆ testRandGaussT()

int testRandGaussT ( )

Definition at line 784 of file testRandDists.cc.

784 {
785
786 cout << "\n--------------------------------------------\n";
787 cout << "Test of RandGaussT distribution \n\n";
788
789 long seed;
790 cout << "Please enter an integer seed: ";
791 cin >> seed; cout << seed << "\n";
792 if (seed == 0) {
793 cout << "Moving on to next test...\n";
794 return 0;
795 }
796
797 int nNumbers;
798 cout << "How many numbers should we generate: ";
799 cin >> nNumbers; cout << nNumbers << "\n";
800
801 double mu;
802 double sigma;
803 cout << "Enter mu: ";
804 cin >> mu; cout << mu << "\n";
805
806 cout << "Enter sigma: ";
807 cin >> sigma; cout << sigma << "\n";
808
809 cout << "\nInstantiating distribution utilizing TripleRand engine...\n";
810 TripleRand eng (seed);
811 RandGaussT dist (eng, mu, sigma);
812
813 cout << "\n Sample fire(): \n";
814 double x;
815
816 x = dist.fire();
817 cout << x;
818
819 cout << "\n Testing operator() ... \n";
820
821 bool good = gaussianTest ( dist, mu, sigma, nNumbers );
822
823 if (good) {
824 return 0;
825 } else {
826 return GaussTBAD;
827 }
828
829} // testRandGaussT()

Referenced by main().

◆ testRandGeneral()

int testRandGeneral ( )

Definition at line 1067 of file testRandDists.cc.

1067 {
1068
1069 cout << "\n--------------------------------------------\n";
1070 cout << "Test of RandGeneral distribution (using a Gaussian shape)\n\n";
1071
1072 bool good;
1073
1074 long seed;
1075 cout << "Please enter an integer seed: ";
1076 cin >> seed; cout << seed << "\n";
1077 if (seed == 0) {
1078 cout << "Moving on to next test...\n";
1079 return 0;
1080 }
1081
1082 int nNumbers;
1083 cout << "How many numbers should we generate: ";
1084 cin >> nNumbers; cout << nNumbers << "\n";
1085
1086 double mu;
1087 double sigma;
1088 mu = .5; // Since randGeneral always ranges from 0 to 1
1089 sigma = .06;
1090
1091 cout << "Enter sigma: ";
1092 cin >> sigma; cout << sigma << "\n";
1093 // We suggest sigma be .06. This leaves room for 8 sigma
1094 // in the distribution. If it is much smaller, the number
1095 // of bins necessary to expect a good match will increase.
1096 // If sigma is much larger, the cutoff before 5 sigma can
1097 // cause the Gaussian hypothesis to be rejected. At .14, for
1098 // example, the 4th moment is 7 sigma away from expectation.
1099
1100 int nBins;
1101 cout << "Enter nBins for stepwise pdf test: ";
1102 cin >> nBins; cout << nBins << "\n";
1103 // We suggest at least 10000 bins; fewer would risk
1104 // false rejection because the step-function curve
1105 // does not match an actual Gaussian. At 10000 bins,
1106 // a million-hit test does not have the resolving power
1107 // to tell the boxy pdf from the true Gaussian. At 5000
1108 // bins, it does.
1109
1110 double xBins = nBins;
1111 double* aProbFunc = new double [nBins];
1112 double x;
1113 for ( int iBin = 0; iBin < nBins; iBin++ ) {
1114 x = iBin / (xBins-1);
1115 aProbFunc [iBin] = std::exp ( - (x-mu)*(x-mu) / (2*sigma*sigma) );
1116 }
1117 // Note that this pdf is not normalized; RandGeneral does that
1118
1119 cout << "\nInstantiating distribution utilizing Ranlux64 engine...\n";
1120 Ranlux64Engine eng (seed, 3);
1121
1122 { // Open block for testing type 1 - step function pdf
1123
1124 RandGeneral dist (eng, aProbFunc, nBins, 1);
1125 delete[] aProbFunc;
1126
1127 double* garbage = new double[nBins];
1128 // We wish to verify that deleting the pdf
1129 // after instantiating the engine is fine.
1130 for ( int gBin = 0; gBin < nBins; gBin++ ) {
1131 garbage [gBin] = 1;
1132 }
1133
1134 cout << "\n Sample fire(): \n";
1135
1136 x = dist.fire();
1137 cout << x;
1138
1139 cout << "\n Testing operator() ... \n";
1140
1141 good = gaussianTest ( dist, mu, sigma, nNumbers );
1142
1143 delete[] garbage;
1144
1145 } // Close block for testing type 1 - step function pdf
1146 // dist goes out of scope but eng is supposed to stick around;
1147 // by closing this block we shall verify that!
1148
1149 cout << "Enter nBins for linearized pdf test: ";
1150 cin >> nBins; cout << nBins << "\n";
1151 // We suggest at least 1000 bins; fewer would risk
1152 // false rejection because the non-smooth curve
1153 // does not match an actual Gaussian. At 1000 bins,
1154 // a million-hit test does not resolve the non-smoothness;
1155 // at 300 bins it does.
1156
1157 xBins = nBins;
1158 aProbFunc = new double [nBins];
1159 for ( int jBin = 0; jBin < nBins; jBin++ ) {
1160 x = jBin / (xBins-1);
1161 aProbFunc [jBin] = std::exp ( - (x-mu)*(x-mu) / (2*sigma*sigma) );
1162 }
1163 // Note that this pdf is not normalized; RandGeneral does that
1164
1165 RandGeneral dist (eng, aProbFunc, nBins, 0);
1166
1167 cout << "\n Sample operator(): \n";
1168
1169 x = dist();
1170 cout << x;
1171
1172 cout << "\n Testing operator() ... \n";
1173
1174 bool good2 = gaussianTest ( dist, mu, sigma, nNumbers );
1175 good = good && good2;
1176
1177 if (good) {
1178 return 0;
1179 } else {
1180 return GeneralBAD;
1181 }
1182
1183} // testRandGeneral()

Referenced by main().

◆ testRandPoisson()

int testRandPoisson ( )

Definition at line 893 of file testRandDists.cc.

893 {
894
895 cout << "\n--------------------------------------------\n";
896 cout << "Test of RandPoisson distribution \n\n";
897
898 long seed;
899 cout << "Please enter an integer seed: ";
900 cin >> seed; cout << seed << "\n";
901 if (seed == 0) {
902 cout << "Moving on to next test...\n";
903 return 0;
904 }
905
906 cout << "\nInstantiating distribution utilizing TripleRand engine...\n";
907 TripleRand eng (seed);
908
909 int nNumbers;
910 cout << "How many numbers should we generate for each mu: ";
911 cin >> nNumbers; cout << nNumbers << "\n";
912
913 bool good = true;
914
915 while (true) {
916 double mu;
917 cout << "Enter a value for mu: ";
918 cin >> mu; cout << mu << "\n";
919 if (mu == 0) break;
920
921 RandPoisson dist (eng, mu);
922
923 cout << "\n Sample fire(): \n";
924 double x;
925
926 x = dist.fire();
927 cout << x;
928
929 cout << "\n Testing operator() ... \n";
930
931 bool this_good = poissonTest ( dist, mu, nNumbers );
932 if (!this_good) {
933 cout << "\n Poisson distribution for mu = " << mu << " is incorrect!!!\n";
934 }
935 good &= this_good;
936 } // end of the while(true)
937
938 if (good) {
939 return 0;
940 } else {
941 return PoissonBAD;
942 }
943
944} // testRandPoisson()
bool poissonTest(RandPoisson &dist, double mu, int N2)

Referenced by main().

◆ testRandPoissonQ()

int testRandPoissonQ ( )

Definition at line 951 of file testRandDists.cc.

951 {
952
953 cout << "\n--------------------------------------------\n";
954 cout << "Test of RandPoissonQ distribution \n\n";
955
956 long seed;
957 cout << "Please enter an integer seed: ";
958 cin >> seed; cout << seed << "\n";
959 if (seed == 0) {
960 cout << "Moving on to next test...\n";
961 return 0;
962 }
963
964 cout << "\nInstantiating distribution utilizing TripleRand engine...\n";
965 TripleRand eng (seed);
966
967 int nNumbers;
968 cout << "How many numbers should we generate for each mu: ";
969 cin >> nNumbers; cout << nNumbers << "\n";
970
971 bool good = true;
972
973 while (true) {
974 double mu;
975 cout << "Enter a value for mu: ";
976 cin >> mu; cout << mu << "\n";
977 if (mu == 0) break;
978
979 RandPoissonQ dist (eng, mu);
980
981 cout << "\n Sample fire(): \n";
982 double x;
983
984 x = dist.fire();
985 cout << x;
986
987 cout << "\n Testing operator() ... \n";
988
989 bool this_good = poissonTest ( dist, mu, nNumbers );
990 if (!this_good) {
991 cout << "\n Poisson distribution for mu = " << mu << " is incorrect!!!\n";
992 }
993 good &= this_good;
994 } // end of the while(true)
995
996 if (good) {
997 return 0;
998 } else {
999 return PoissonQBAD;
1000 }
1001
1002} // testRandPoissonQ()

Referenced by main().

◆ testRandPoissonT()

int testRandPoissonT ( )

Definition at line 1009 of file testRandDists.cc.

1009 {
1010
1011 cout << "\n--------------------------------------------\n";
1012 cout << "Test of RandPoissonT distribution \n\n";
1013
1014 long seed;
1015 cout << "Please enter an integer seed: ";
1016 cin >> seed; cout << seed << "\n";
1017 if (seed == 0) {
1018 cout << "Moving on to next test...\n";
1019 return 0;
1020 }
1021
1022 cout << "\nInstantiating distribution utilizing TripleRand engine...\n";
1023 TripleRand eng (seed);
1024
1025 int nNumbers;
1026 cout << "How many numbers should we generate for each mu: ";
1027 cin >> nNumbers; cout << nNumbers << "\n";
1028
1029 bool good = true;
1030
1031 while (true) {
1032 double mu;
1033 cout << "Enter a value for mu: ";
1034 cin >> mu; cout << mu << "\n";
1035 if (mu == 0) break;
1036
1037 RandPoissonT dist (eng, mu);
1038
1039 cout << "\n Sample fire(): \n";
1040 double x;
1041
1042 x = dist.fire();
1043 cout << x;
1044
1045 cout << "\n Testing operator() ... \n";
1046
1047 bool this_good = poissonTest ( dist, mu, nNumbers );
1048 if (!this_good) {
1049 cout << "\n Poisson distribution for mu = " << mu << " is incorrect!!!\n";
1050 }
1051 good &= this_good;
1052 } // end of the while(true)
1053
1054 if (good) {
1055 return 0;
1056 } else {
1057 return PoissonTBAD;
1058 }
1059
1060} // testRandPoissonT()

Referenced by main().

◆ testSkewNormal()

int testSkewNormal ( )

Definition at line 735 of file testRandDists.cc.

735 {
736
737 cout << "\n--------------------------------------------\n";
738 cout << "Test of SkewNormal distribution \n\n";
739
740 long seed;
741 cout << "Please enter an integer seed: ";
742 cin >> seed; cout << seed << "\n";
743 if (seed == 0) {
744 cout << "Moving on to next test...\n";
745 return 0;
746 }
747
748 int nNumbers;
749 cout << "How many numbers should we generate: ";
750 cin >> nNumbers; cout << nNumbers << "\n";
751
752 double k;
753 cout << "Enter k: ";
754 cin >> k; cout << k << "\n";
755
756 cout << "\nInstantiating distribution utilizing TripleRand engine...\n";
757 TripleRand eng (seed);
758 RandSkewNormal dist (eng, k);
759
760 cout << "\n Sample fire(): \n";
761 double x;
762
763 x = dist.fire();
764 cout << x;
765
766 cout << "\n Testing operator() ... \n";
767
768 bool good = skewNormalTest ( dist, k, nNumbers );
769
770 if (good) {
771 return 0;
772 } else {
773 return SkewNormalBAD;
774 }
775
776} // testSkewNormal()
bool skewNormalTest(HepRandom &dist, double k, int nNumbers)

Referenced by main().