39 #include "CLHEP/Units/GlobalPhysicalConstants.h"
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"
57 using std::setprecision;
58 using namespace CLHEP;
64 static const double REJECT = 4.0;
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;
98 static double c[6] = {
103 0.001208650973866179,
104 -0.000005395239384953 };
107 tmp -= (
x+.5)*std::log(tmp);
108 ser = 1.000000000190015;
109 for (
int i = 0;
i < 6;
i++) {
112 double ans = (-tmp + std::log (std::sqrt(CLHEP::twopi)*ser/
x));
117 double gser(
double a,
double x) {
118 const int ITMAX = 100;
119 const double EPS = 1.0E-8;
123 for (
int n=0;
n < ITMAX;
n++) {
127 if (std::fabs(del) < std::fabs(sum)*EPS) {
128 return sum*std::exp(-
x+
a*std::log(
x)-
gammln(
a));
131 cout <<
"Problem - inaccurate gser " <<
a <<
", " <<
x <<
"\n";
132 return sum*std::exp(-
x+
a*std::log(
x)-
gammln(
a));
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;
141 double c = 1/VERYSMALL;
144 for (
int i = 1;
i <= ITMAX;
i++) {
145 double an = -
i*(
i-
a);
148 if (std::fabs(d) < VERYSMALL) d = VERYSMALL;
150 if (std::fabs(c) < VERYSMALL) c = VERYSMALL;
154 if (std::fabs(del-1.0) < EPS) {
155 return std::exp(-
x+
a*std::log(
x)-
gammln(
a))*h;
158 cout <<
"Problem - inaccurate gcf " <<
a <<
", " <<
x <<
"\n";
159 return std::exp(-
x+
a*std::log(
x)-
gammln(
a))*h;
163 double gammp (
double a,
double x) {
183 double sigma,
int nNumbers ) {
186 double worstSigma = 0;
204 for (ciu = 0; ciu < 11; ciu++) {
209 int oldprecision = cout.precision();
214 int ipr = nNumbers / 10 + 1;
215 for (
int ifire = 0; ifire < nNumbers; ifire++) {
218 if(
x <
mu - 12.0*sigma ) {
219 cout <<
"x = " <<
x <<
"\n";
221 if ( (ifire % ipr) == 0 ) {
222 cout << ifire << endl;
230 u = (
x -
mu) / sigma;
233 if (ciu>10) ciu = 10;
237 if (ciu>10) ciu = 10;
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;
257 cout <<
"Mean (should be close to " <<
mu <<
"): " << mean << endl;
258 cout <<
"Second moment (should be close to " << sigma*sigma <<
260 cout <<
"Third moment (should be close to zero): " << u3 << endl;
261 cout <<
"Fourth moment (should be close to " << 3*sigma*sigma*sigma*sigma <<
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;
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);
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;
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";
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 "
321 << ncounts[m1] <<
" negative and " << counts[m1] <<
" positive " <<
"\n";
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";
331 if ( negSigs > worstSigma ) worstSigma = negSigs;
332 if ( posSigs > worstSigma ) worstSigma = posSigs;
335 cout <<
"\n The worst deviation encountered (out of about 25) was "
336 << worstSigma <<
" sigma \n\n";
338 cout.precision(oldprecision);
353 double worstSigma = 0;
368 int oldprecision = cout.precision();
373 double delta =
k / std::sqrt( 1 +
k*
k );
374 double mu =
delta/std::sqrt(CLHEP::halfpi);
381 int ipr = nNumbers / 10 + 1;
382 for (
int ifire = 0; ifire < nNumbers; ifire++) {
385 if(
x <
mu - 12.0 ) {
386 cout <<
"x = " <<
x <<
"\n";
388 if ( (ifire % ipr) == 0 ) {
389 cout << ifire << endl;
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;
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;
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);
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;
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";
437 cout <<
"\n The worst deviation encountered (out of about 25) was "
438 << worstSigma <<
" sigma \n\n";
440 cout.precision(oldprecision);
456 double logAnswer = -mu_ + r*std::log(mu_) -
gammln(r+1);
457 return std::exp(logAnswer);
462 int MINBIN,
int MAXBINS,
int clumping,
463 int& firstBin,
int& lastBin ) {
470 double * refdist =
new double [MAXBINS];
483 while ( c < MAXBINS ) {
484 for ( ic=0, binc=0; ic < clumping; ic++, r++ ) {
485 binc += pdist(r) *
N;
488 if (binc >= MINBIN)
break;
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";
497 refdist[firstBin] = start;
502 while ( c < MAXBINS ) {
503 for ( ic=0, binc=0; ic < clumping; ic++, r++ ) {
504 binc += pdist(r) *
N;
507 if (next < MINBIN)
break;
514 next += refdist[lastBin];
515 while ( c < MAXBINS ) {
516 for ( ic=0, binc=0; ic < clumping; ic++, r++ ) {
517 binc += pdist(r) *
N;
522 refdist[lastBin] = next;
544 int clumping = int(std::sqrt(
mu));
545 if (clumping <= 1) clumping = 2;
546 const int MINBIN = 20;
547 const int MAXBINS = 1000;
556 MINBIN, MAXBINS, 1, firstBin, lastBin);
558 MINBIN, MAXBINS, clumping, firstBin2, lastBin2);
565 double*
samples =
new double [MAXBINS];
566 double* samples2 =
new double [MAXBINS];
568 for (r = 0; r < MAXBINS; r++) {
575 for (
int i = 0;
i <
N;
i++) {
578 moment += (r -
mu)*(r -
mu);
580 if (r1 < firstBin) r1 = firstBin;
581 if (r1 > lastBin) r1 = lastBin;
584 if (r2 < firstBin2) r2 = firstBin2;
585 if (r2 > lastBin2) r2 = lastBin2;
591 for (
k = firstBin;
k <= lastBin;
k++) {
592 cout <<
k <<
" " <<
samples[
k] <<
" " << refdist[
k] <<
" " <<
596 for (
k = firstBin2;
k <= lastBin2;
k++) {
597 cout <<
k <<
" " << samples2[
k] <<
" " << refdist2[
k] <<
"\n";
604 for ( r = firstBin; r <= lastBin; r++ ) {
608 int degFreedom = (lastBin - firstBin + 1) - 1;
617 pval = 1.0 - gammp ( .5*degFreedom , .5*chi2 );
619 cout <<
"Chi^2 is " << chi2 <<
" on " << degFreedom <<
" degrees of freedom."
620 <<
" p = " << pval <<
"\n";
628 for ( r = firstBin2; r <= lastBin2; r++ ) {
629 double delta = (samples2[r] - refdist2[r]);
632 degFreedom = (lastBin2 - firstBin2 + 1) - 1;
634 pval2 = 1.0 - gammp ( .5*degFreedom , .5*chi2 );
636 cout <<
"Clumps: Chi^2 is " << chi2 <<
" on " << degFreedom <<
637 " degrees of freedom." <<
" p = " << pval2 <<
"\n";
644 double mean = sum /
N;
645 double sigma = std::sqrt( moment / (
N-1) );
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);
651 cout <<
"Mean (should be " <<
mu <<
") is " << mean <<
"\n";
652 cout <<
"Sigma (should be " << std::sqrt(
mu) <<
") is " << sigma <<
"\n";
654 cout <<
"These are " << deviationMean <<
" and " << deviationSigma <<
655 " standard deviations from expected values\n\n";
663 if ( (pval < .0001) || (pval2 < .0001) ||
664 (deviationMean > 3.5) || (deviationSigma > 3.5) ) {
666 cout <<
"REJECT this distributon!!!\n";
686 cout <<
"\n--------------------------------------------\n";
687 cout <<
"Test of RandGauss distribution \n\n";
690 cout <<
"Please enter an integer seed: ";
691 cin >> seed; cout << seed <<
"\n";
693 cout <<
"Moving on to next test...\n";
698 cout <<
"How many numbers should we generate: ";
699 cin >> nNumbers; cout << nNumbers <<
"\n";
703 cout <<
"Enter mu: ";
704 cin >>
mu; cout <<
mu <<
"\n";
706 cout <<
"Enter sigma: ";
707 cin >> sigma; cout << sigma <<
"\n";
709 cout <<
"\nInstantiating distribution utilizing TripleRand engine...\n";
713 cout <<
"\n Sample fire(): \n";
719 cout <<
"\n Testing operator() ... \n";
739 cout <<
"\n--------------------------------------------\n";
740 cout <<
"Test of SkewNormal distribution \n\n";
743 cout <<
"Please enter an integer seed: ";
744 cin >> seed; cout << seed <<
"\n";
746 cout <<
"Moving on to next test...\n";
751 cout <<
"How many numbers should we generate: ";
752 cin >> nNumbers; cout << nNumbers <<
"\n";
756 cin >>
k; cout <<
k <<
"\n";
758 cout <<
"\nInstantiating distribution utilizing TripleRand engine...\n";
762 cout <<
"\n Sample fire(): \n";
768 cout <<
"\n Testing operator() ... \n";
775 return SkewNormalBAD;
788 cout <<
"\n--------------------------------------------\n";
789 cout <<
"Test of RandGaussT distribution \n\n";
792 cout <<
"Please enter an integer seed: ";
793 cin >> seed; cout << seed <<
"\n";
795 cout <<
"Moving on to next test...\n";
800 cout <<
"How many numbers should we generate: ";
801 cin >> nNumbers; cout << nNumbers <<
"\n";
805 cout <<
"Enter mu: ";
806 cin >>
mu; cout <<
mu <<
"\n";
808 cout <<
"Enter sigma: ";
809 cin >> sigma; cout << sigma <<
"\n";
811 cout <<
"\nInstantiating distribution utilizing TripleRand engine...\n";
815 cout <<
"\n Sample fire(): \n";
821 cout <<
"\n Testing operator() ... \n";
841 cout <<
"\n--------------------------------------------\n";
842 cout <<
"Test of RandGaussQ distribution \n\n";
845 cout <<
"Please enter an integer seed: ";
846 cin >> seed; cout << seed <<
"\n";
848 cout <<
"Moving on to next test...\n";
853 cout <<
"How many numbers should we generate: ";
854 cin >> nNumbers; cout << nNumbers <<
"\n";
856 if (nNumbers >= 20000000) {
857 cout <<
"With that many samples RandGaussQ need not pass validation...\n";
862 cout <<
"Enter mu: ";
863 cin >>
mu; cout <<
mu <<
"\n";
865 cout <<
"Enter sigma: ";
866 cin >> sigma; cout << sigma <<
"\n";
868 cout <<
"\nInstantiating distribution utilizing DualRand engine...\n";
872 cout <<
"\n Sample fire(): \n";
878 cout <<
"\n Testing operator() ... \n";
897 cout <<
"\n--------------------------------------------\n";
898 cout <<
"Test of RandPoisson distribution \n\n";
901 cout <<
"Please enter an integer seed: ";
902 cin >> seed; cout << seed <<
"\n";
904 cout <<
"Moving on to next test...\n";
908 cout <<
"\nInstantiating distribution utilizing TripleRand engine...\n";
912 cout <<
"How many numbers should we generate for each mu: ";
913 cin >> nNumbers; cout << nNumbers <<
"\n";
919 cout <<
"Enter a value for mu: ";
920 cin >>
mu; cout <<
mu <<
"\n";
925 cout <<
"\n Sample fire(): \n";
931 cout <<
"\n Testing operator() ... \n";
935 cout <<
"\n Poisson distribution for mu = " <<
mu <<
" is incorrect!!!\n";
955 cout <<
"\n--------------------------------------------\n";
956 cout <<
"Test of RandPoissonQ distribution \n\n";
959 cout <<
"Please enter an integer seed: ";
960 cin >> seed; cout << seed <<
"\n";
962 cout <<
"Moving on to next test...\n";
966 cout <<
"\nInstantiating distribution utilizing TripleRand engine...\n";
970 cout <<
"How many numbers should we generate for each mu: ";
971 cin >> nNumbers; cout << nNumbers <<
"\n";
977 cout <<
"Enter a value for mu: ";
978 cin >>
mu; cout <<
mu <<
"\n";
983 cout <<
"\n Sample fire(): \n";
989 cout <<
"\n Testing operator() ... \n";
993 cout <<
"\n Poisson distribution for mu = " <<
mu <<
" is incorrect!!!\n";
1013 cout <<
"\n--------------------------------------------\n";
1014 cout <<
"Test of RandPoissonT distribution \n\n";
1017 cout <<
"Please enter an integer seed: ";
1018 cin >> seed; cout << seed <<
"\n";
1020 cout <<
"Moving on to next test...\n";
1024 cout <<
"\nInstantiating distribution utilizing TripleRand engine...\n";
1028 cout <<
"How many numbers should we generate for each mu: ";
1029 cin >> nNumbers; cout << nNumbers <<
"\n";
1035 cout <<
"Enter a value for mu: ";
1036 cin >>
mu; cout <<
mu <<
"\n";
1041 cout <<
"\n Sample fire(): \n";
1047 cout <<
"\n Testing operator() ... \n";
1051 cout <<
"\n Poisson distribution for mu = " <<
mu <<
" is incorrect!!!\n";
1071 cout <<
"\n--------------------------------------------\n";
1072 cout <<
"Test of RandGeneral distribution (using a Gaussian shape)\n\n";
1077 cout <<
"Please enter an integer seed: ";
1078 cin >> seed; cout << seed <<
"\n";
1080 cout <<
"Moving on to next test...\n";
1085 cout <<
"How many numbers should we generate: ";
1086 cin >> nNumbers; cout << nNumbers <<
"\n";
1093 cout <<
"Enter sigma: ";
1094 cin >> sigma; cout << sigma <<
"\n";
1103 cout <<
"Enter nBins for stepwise pdf test: ";
1104 cin >> nBins; cout << nBins <<
"\n";
1112 double xBins = nBins;
1113 double* aProbFunc =
new double [nBins];
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) );
1121 cout <<
"\nInstantiating distribution utilizing Ranlux64 engine...\n";
1129 double* garbage =
new double[nBins];
1132 for (
int gBin = 0; gBin < nBins; gBin++ ) {
1136 cout <<
"\n Sample fire(): \n";
1141 cout <<
"\n Testing operator() ... \n";
1151 cout <<
"Enter nBins for linearized pdf test: ";
1152 cin >> nBins; cout << nBins <<
"\n";
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) );
1169 cout <<
"\n Sample operator(): \n";
1174 cout <<
"\n Testing operator() ... \n";
1177 good = good && good2;
1215 return mask > 0 ? -mask : mask;