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

RandPoissonQ.cc
Go to the documentation of this file.
1// $Id: RandPoissonQ.cc,v 1.7 2010/06/16 17:24:53 garren Exp $
2// -*- C++ -*-
3//
4// -----------------------------------------------------------------------
5// HEP Random
6// --- RandPoissonQ ---
7// class implementation file
8// -----------------------------------------------------------------------
9
10// =======================================================================
11// M. Fischler - Implemented new, much faster table-driven algorithm
12// applicable for mu < 100
13// - Implemented "quick()" methods, shich are the same as the
14// new methods for mu < 100 and are a skew-corrected gaussian
15// approximation for large mu.
16// M. Fischler - Removed mean=100 from the table-driven set, since it
17// uses a value just off the end of the table. (April 2004)
18// M Fischler - put and get to/from streams 12/15/04
19// M Fischler - Utilize RandGaussQ rather than RandGauss, as clearly
20// intended by the inclusion of RandGaussQ.h. Using RandGauss
21// introduces a subtle trap in that the state of RandPoissonQ
22// can never be properly captured without also saveing the
23// state of RandGauss! RandGaussQ is, on the other hand,
24// stateless except for the engine used.
25// M Fisculer - Modified use of wrong engine when shoot (anEngine, mean)
26// is called. This flaw was preventing any hope of proper
27// saving and restoring in the instance cases.
28// M Fischler - fireArray using defaultMean 2/10/05
29// M Fischler - put/get to/from streams uses pairs of ulongs when
30// + storing doubles avoid problems with precision
31// 4/14/05
32// M Fisculer - Modified use of shoot (mean) instead of
33// shoot(getLocalEngine(), mean) when fire(mean) is called.
34// This flaw was causing bad "cross-talk" between modules
35// in CMS, where one used its own engine, and the other
36// used the static generator. 10/18/07
37//
38// =======================================================================
39
40#include "CLHEP/Random/defs.h"
44#include "CLHEP/Random/Stat.h"
45#include <cmath> // for std::pow()
46
47namespace CLHEP {
48
49std::string RandPoissonQ::name() const {return "RandPoissonQ";}
51
52// Initialization of static data: Note that this is all const static data,
53// so that saveEngineStatus properly saves all needed information.
54
55 // The following MUST MATCH the corresponding values used (in
56 // poissonTables.cc) when poissonTables.cdat was created.
57
58const double RandPoissonQ::FIRST_MU = 10;// lowest mu value in table
59const double RandPoissonQ::LAST_MU = 95;// highest mu value
60const double RandPoissonQ::S = 5; // Spacing between mu values
61const int RandPoissonQ::BELOW = 30; // Starting point for N is at mu - BELOW
62const int RandPoissonQ::ENTRIES = 51; // Number of entries in each mu row
63
64const double RandPoissonQ::MAXIMUM_POISSON_DEVIATE = 2.0E9;
65 // Careful -- this is NOT the maximum number that can be held in
66 // a long. It actually should be some large number of sigma below
67 // that.
68
69 // Here comes the big (9K bytes) table, kept in a file of
70 // ENTRIES * (FIRST_MU - LAST_MU + 1)/S doubles
71
72static const double poissonTables [ 51 * ( (95-10)/5 + 1 ) ] = {
73#include "poissonTables.cdat"
74};
75
76
77//
78// Constructors and destructors:
79//
80
83
84void RandPoissonQ::setupForDefaultMu() {
85
86 // The following are useful for quick approximation, for large mu
87
88 double sig2 = defaultMean * (.9998654 - .08346/defaultMean);
89 sigma = std::sqrt(sig2);
90 // sigma for the Guassian which approximates the Poisson -- naively
91 // std::sqrt (defaultMean).
92 //
93 // The multiplier corrects for fact that discretization of the form
94 // [gaussian+.5] increases the second moment by a small amount.
95
96 double t = 1./(sig2);
97
98 a2 = t/6 + t*t/324;
99 a1 = std::sqrt (1-2*a2*a2*sig2);
100 a0 = defaultMean + .5 - sig2 * a2;
101
102 // The formula will be a0 + a1*x + a2*x*x where x has 2nd moment of sigma.
103 // The coeffeicients are chosen to match the first THREE moments of the
104 // true Poisson distribution.
105 //
106 // Actually, if the correction for discretization were not needed, then
107 // a2 could be taken one order higher by adding t*t*t/5832. However,
108 // the discretization correction is not perfect, leading to inaccuracy
109 // on the order to 1/mu**2, so adding a third term is overkill.
110
111} // setupForDefaultMu()
112
113
114//
115// fire, quick, operator(), and shoot methods:
116//
117
118long RandPoissonQ::shoot(double xm) {
119 return shoot(getTheEngine(), xm);
120}
121
123 return (double) fire();
124}
125
126double RandPoissonQ::operator()( double mean ) {
127 return (double) fire(mean);
128}
129
130long RandPoissonQ::fire(double mean) {
131 return shoot(getLocalEngine(), mean);
132}
133
135 if ( defaultMean < LAST_MU + S ) {
136 return poissonDeviateSmall ( getLocalEngine(), defaultMean );
137 } else {
138 return poissonDeviateQuick ( getLocalEngine(), a0, a1, a2, sigma );
139 }
140} // fire()
141
142long RandPoissonQ::shoot(HepRandomEngine* anEngine, double mean) {
143
144 // The following variables, static to this method, apply to the
145 // last time a large mean was supplied; they obviate certain calculations
146 // if consecutive calls use the same mean.
147
148 static double lastLargeMean = -1.; // Mean from previous shoot
149 // requiring poissonDeviateQuick()
150 static double lastA0;
151 static double lastA1;
152 static double lastA2;
153 static double lastSigma;
154
155 if ( mean < LAST_MU + S ) {
156 return poissonDeviateSmall ( anEngine, mean );
157 } else {
158 if ( mean != lastLargeMean ) {
159 // Compute the coefficients defining the quadratic transformation from a
160 // Gaussian to a Poisson for this mean. Also save these for next time.
161 double sig2 = mean * (.9998654 - .08346/mean);
162 lastSigma = std::sqrt(sig2);
163 double t = 1./sig2;
164 lastA2 = t*(1./6.) + t*t*(1./324.);
165 lastA1 = std::sqrt (1-2*lastA2*lastA2*sig2);
166 lastA0 = mean + .5 - sig2 * lastA2;
167 }
168 return poissonDeviateQuick ( anEngine, lastA0, lastA1, lastA2, lastSigma );
169 }
170
171} // shoot (anEngine, mean)
172
173void RandPoissonQ::shootArray(const int size, long* vect, double m) {
174 for( long* v = vect; v != vect + size; ++v )
175 *v = shoot(m);
176 // Note: We could test for m > 100, and if it is, precompute a0, a1, a2,
177 // and sigma and call the appropriate form of poissonDeviateQuick.
178 // But since those are cached anyway, not much time would be saved.
179}
180
181void RandPoissonQ::fireArray(const int size, long* vect, double m) {
182 for( long* v = vect; v != vect + size; ++v )
183 *v = fire( m );
184}
185
186void RandPoissonQ::fireArray(const int size, long* vect) {
187 for( long* v = vect; v != vect + size; ++v )
188 *v = fire( defaultMean );
189}
190
191
192// Quick Poisson deviate algorithm used by quick for large mu:
193
194long RandPoissonQ::poissonDeviateQuick ( HepRandomEngine *e, double mu ) {
195
196 // Compute the coefficients defining the quadratic transformation from a
197 // Gaussian to a Poisson:
198
199 double sig2 = mu * (.9998654 - .08346/mu);
200 double sig = std::sqrt(sig2);
201 // The multiplier corrects for fact that discretization of the form
202 // [gaussian+.5] increases the second moment by a small amount.
203
204 double t = 1./sig2;
205
206 double sa2 = t*(1./6.) + t*t*(1./324.);
207 double sa1 = std::sqrt (1-2*sa2*sa2*sig2);
208 double sa0 = mu + .5 - sig2 * sa2;
209
210 // The formula will be sa0 + sa1*x + sa2*x*x where x has sigma of sq.
211 // The coeffeicients are chosen to match the first THREE moments of the
212 // true Poisson distribution.
213
214 return poissonDeviateQuick ( e, sa0, sa1, sa2, sig );
215}
216
217
218long RandPoissonQ::poissonDeviateQuick ( HepRandomEngine *e,
219 double A0, double A1, double A2, double sig) {
220//
221// Quick Poisson deviate algorithm used by quick for large mu:
222//
223// The principle: For very large mu, a poisson distribution can be approximated
224// by a gaussian: return the integer part of mu + .5 + g where g is a unit
225// normal. However, this yelds a miserable approximation at values as
226// "large" as 100. The primary problem is that the poisson distribution is
227// supposed to have a skew of 1/mu**2, and the zero skew of the Guassian
228// leads to errors of order as big as 1/mu**2.
229//
230// We substitute for the gaussian a quadratic function of that gaussian random.
231// The expression looks very nearly like mu + .5 - 1/6 + g + g**2/(6*mu).
232// The small positive quadratic term causes the resulting variate to have
233// a positive skew; the -1/6 constant term is there to correct for this bias
234// in the mean. By adjusting these two and the linear term, we can match the
235// first three moments to high accuracy in 1/mu.
236//
237// The sigma used is not precisely std::sqrt(mu) since a rounded-off Gaussian
238// has a second moment which is slightly larger than that of the Gaussian.
239// To compensate, sig is multiplied by a factor which is slightly less than 1.
240
241 // double g = RandGauss::shootQuick( e ); // TEMPORARY MOD:
242 double g = RandGaussQ::shoot( e ); // Unit normal
243 g *= sig;
244 double p = A2*g*g + A1*g + A0;
245 if ( p < 0 ) return 0; // Shouldn't ever possibly happen since
246 // mean should not be less than 100, but
247 // we check due to paranoia.
249
250 return long(p);
251
252} // poissonDeviateQuick ()
253
254
255
256long RandPoissonQ::poissonDeviateSmall (HepRandomEngine * e, double mean) {
257 long N1;
258 long N2;
259 // The following are for later use to form a secondary random s:
260 double rRange; // This will hold the interval between cdf for the
261 // computed N1 and cdf for N1+1.
262 double rRemainder = 0; // This will hold the length into that interval.
263
264 // Coming in, mean should not be more than LAST_MU + S. However, we will
265 // be paranoid and test for this:
266
267 if ( mean > LAST_MU + S ) {
268 return RandPoisson::shoot(e, mean);
269 }
270
271 if (mean <= 0) {
272 return 0; // Perhaps we ought to balk harder here!
273 }
274
275 // >>> 1 <<<
276 // Generate the first random, which we always will need.
277
278 double r = e->flat();
279
280 // >>> 2 <<<
281 // For small mean, below the start of the tables,
282 // do the series for cdf directly.
283
284 // In this case, since we know the series will terminate relatively quickly,
285 // almost alwaye use a precomputed 1/N array without fear of overrunning it.
286
287 static const double oneOverN[50] =
288 { 0, 1., 1/2., 1/3., 1/4., 1/5., 1/6., 1/7., 1/8., 1/9.,
289 1/10., 1/11., 1/12., 1/13., 1/14., 1/15., 1/16., 1/17., 1/18., 1/19.,
290 1/20., 1/21., 1/22., 1/23., 1/24., 1/25., 1/26., 1/27., 1/28., 1/29.,
291 1/30., 1/31., 1/32., 1/33., 1/34., 1/35., 1/36., 1/37., 1/38., 1/39.,
292 1/40., 1/41., 1/42., 1/43., 1/44., 1/45., 1/46., 1/47., 1/48., 1/49. };
293
294
295 if ( mean < FIRST_MU ) {
296
297 long N = 0;
298 double term = std::exp(-mean);
299 double cdf = term;
300
301 if ( r < (1 - 1.0E-9) ) {
302 //
303 // **** This is a normal path: ****
304 //
305 // Except when r is very close to 1, it is certain that we will exceed r
306 // before the 30-th term in the series, so a simple while loop is OK.
307 const double* oneOverNptr = oneOverN;
308 while( cdf <= r ) {
309 N++ ;
310 oneOverNptr++;
311 term *= ( mean * (*oneOverNptr) );
312 cdf += term;
313 }
314 return N;
315 //
316 // **** ****
317 //
318 } else { // r is almost 1...
319 // For r very near to 1 we would have to check that we don't fall
320 // off the end of the table of 1/N. Since this is very rare, we just
321 // ignore the table and do the identical while loop, using explicit
322 // division.
323 double cdf0;
324 while ( cdf <= r ) {
325 N++ ;
326 term *= ( mean / N );
327 cdf0 = cdf;
328 cdf += term;
329 if (cdf == cdf0) break; // Can't happen, but just in case...
330 }
331 return N;
332 } // end of if ( r compared to (1 - 1.0E-9) )
333
334 } // End of the code for mean < FIRST_MU
335
336 // >>> 3 <<<
337 // Find the row of the tables corresponding to the highest tabulated mu
338 // which is no greater than our actual mean.
339
340 int rowNumber = int((mean - FIRST_MU)/S);
341 const double * cdfs = &poissonTables [rowNumber*ENTRIES];
342 double mu = FIRST_MU + rowNumber*S;
343 double deltaMu = mean - mu;
344 int Nmin = int(mu - BELOW);
345 if (Nmin < 1) Nmin = 1;
346 int Nmax = Nmin + (ENTRIES - 1);
347
348
349 // >>> 4 <<<
350 // If r is less that the smallest entry in the row, then
351 // generate the deviate directly from the series.
352
353 if ( r < cdfs[0] ) {
354
355 // In this case, we are tempted to use the actual mean, and not
356 // generate a second deviate to account for the leftover part mean - mu.
357 // That would be an error, generating a distribution with enough excess
358 // at Nmin + (mean-mu)/2 to be detectable in 4,000,000 trials.
359
360 // Since this case is very rare (never more than .2% of the r values)
361 // and can happen where N will be large (up to 65 for the mu=95 row)
362 // we use explicit division so as to avoid having to worry about running
363 // out of oneOverN table.
364
365 long N = 0;
366 double term = std::exp(-mu);
367 double cdf = term;
368 double cdf0;
369
370 while(cdf <= r) {
371 N++ ;
372 term *= ( mu / N );
373 cdf0 = cdf;
374 cdf += term;
375 if (cdf == cdf0) break; // Can't happen, but just in case...
376 }
377 N1 = N;
378 // std::cout << r << " " << N << " ";
379 // DBG_small = true;
380 rRange = 0; // In this case there is always a second r needed
381
382 } // end of small-r case
383
384
385 // >>> 5 <<<
386 // Assuming r lies within the scope of the row for this mu, find the
387 // largest entry not greater than r. N1 is the N corresponding to the
388 // index a.
389
390 else if ( r < cdfs[ENTRIES-1] ) { // r is also >= cdfs[0]
391
392 //
393 // **** This is the normal code path ****
394 //
395
396 int a = 0; // Highest value of index such that cdfs[a]
397 // is known NOT to be greater than r.
398 int b = ENTRIES - 1; // Lowest value of index such that cdfs[b] is
399 // known to exeed r.
400
401 while (b != (a+1) ) {
402 int c = (a+b+1)>>1;
403 if (r > cdfs[c]) {
404 a = c;
405 } else {
406 b = c;
407 }
408 }
409
410 N1 = Nmin + a;
411 rRange = cdfs[a+1] - cdfs[a];
412 rRemainder = r - cdfs[a];
413
414 //
415 // **** ****
416 //
417
418 } // end of medium-r (normal) case
419
420
421 // >>> 6 <<<
422 // If r exceeds the greatest entry in the table for this mu, then start
423 // from that cdf, and use the series to compute from there until r is
424 // exceeded.
425
426 else { // if ( r >= cdfs[ENTRIES-1] ) {
427
428 // Here, division must be done explicitly, and we must also protect against
429 // roundoff preventing termination.
430
431 //
432 //+++ cdfs[ENTRIES-1] is std::exp(-mu) sum (mu**m/m! , m=0 to Nmax)
433 //+++ (where Nmax = mu - BELOW + ENTRIES - 1)
434 //+++ cdfs[ENTRIES-1]-cdfs[ENTRIES-2] is std::exp(-mu) mu**(Nmax)/(Nmax)!
435 //+++ If the sum up to k-1 <= r < sum up to k, then N = k-1
436 //+++ Consider k = Nmax in the above statement:
437 //+++ If cdfs[ENTRIES-2] <= r < cdfs[ENTRIES-1], N would be Nmax-1
438 //+++ But here r >= cdfs[ENTRIES-1] so N >= Nmax
439 //
440
441 // Erroneous:
442 //+++ cdfs[ENTRIES-1] is std::exp(-mu) sum (mu**m/m! , m=0 to Nmax-1)
443 //+++ cdfs[ENTRIES-1]-cdfs[ENTRIES-2] is std::exp(-mu) mu**(Nmax-1)/(Nmax-1)!
444 //+++ If a sum up to k-1 <= r < sum up to k, then N = k-1
445 //+++ So if cdfs[ENTRIES-1] were > r, N would be Nmax-1 (or less)
446 //+++ But here r >= cdfs[ENTRIES-1] so N >= Nmax
447 //
448
449 // std::cout << "r = " << r << " mu = " << mu << "\n";
450 long N = Nmax -1;
451 double cdf = cdfs[ENTRIES-1];
452 double term = cdf - cdfs[ENTRIES-2];
453 double cdf0;
454 while(cdf <= r) {
455 N++ ;
456 // std::cout << " N " << N << " term " <<
457 // term << " cdf " << cdf << "\n";
458 term *= ( mu / N );
459 cdf0 = cdf;
460 cdf += term;
461 if (cdf == cdf0) break; // If term gets so small cdf stops increasing,
462 // terminate using that value of N since we
463 // would never reach r.
464 }
465 N1 = N;
466 rRange = 0; // We can't validly omit the second true random
467
468 // N = Nmax -1;
469 // cdf = cdfs[ENTRIES-1];
470 // term = cdf - cdfs[ENTRIES-2];
471 // for (int isxz=0; isxz < 100; isxz++) {
472 // N++ ;
473 // term *= ( mu / N );
474 // cdf0 = cdf;
475 // cdf += term;
476 // }
477 // std::cout.precision(20);
478 // std::cout << "Final sum is " << cdf << "\n";
479
480 } // end of large-r case
481
482
483
484 // >>> 7 <<<
485 // Form a second random, s, based on the position of r within the range
486 // of this table entry to the next entry.
487
488 // However, if this range is very small, then we lose too many bits of
489 // randomness. In that situation, we generate a second random for s.
490
491 double s;
492
493 static const double MINRANGE = .01; // Sacrifice up to two digits of
494 // randomness when using r to produce
495 // a second random s. Leads to up to
496 // .09 extra randoms each time.
497
498 if ( rRange > MINRANGE ) {
499 //
500 // **** This path taken 90% of the time ****
501 //
502 s = rRemainder / rRange;
503 } else {
504 s = e->flat(); // extra true random needed about one time in 10.
505 }
506
507 // >>> 8 <<<
508 // Use the direct summation method to form a second poisson deviate N2
509 // from deltaMu and s.
510
511 N2 = 0;
512 double term = std::exp(-deltaMu);
513 double cdf = term;
514
515 if ( s < (1 - 1.0E-10) ) {
516 //
517 // This is the normal path:
518 //
519 const double* oneOverNptr = oneOverN;
520 while( cdf <= s ) {
521 N2++ ;
522 oneOverNptr++;
523 term *= ( deltaMu * (*oneOverNptr) );
524 cdf += term;
525 }
526 } else { // s is almost 1...
527 while( cdf <= s ) {
528 N2++ ;
529 term *= ( deltaMu / N2 );
530 cdf += term;
531 }
532 } // end of if ( s compared to (1 - 1.0E-10) )
533
534 // >>> 9 <<<
535 // The result is the sum of those two deviates
536
537 // if (DBG_small) {
538 // std::cout << N2 << " " << N1+N2 << "\n";
539 // DBG_small = false;
540 // }
541
542 return N1 + N2;
543
544} // poissonDeviate()
545
546std::ostream & RandPoissonQ::put ( std::ostream & os ) const {
547 int pr=os.precision(20);
548 std::vector<unsigned long> t(2);
549 os << " " << name() << "\n";
550 os << "Uvec" << "\n";
551 t = DoubConv::dto2longs(a0);
552 os << a0 << " " << t[0] << " " << t[1] << "\n";
553 t = DoubConv::dto2longs(a1);
554 os << a1 << " " << t[0] << " " << t[1] << "\n";
555 t = DoubConv::dto2longs(a2);
556 os << a2 << " " << t[0] << " " << t[1] << "\n";
557 t = DoubConv::dto2longs(sigma);
558 os << sigma << " " << t[0] << " " << t[1] << "\n";
560 os.precision(pr);
561 return os;
562#ifdef REMOVED
563 int pr=os.precision(20);
564 os << " " << name() << "\n";
565 os << a0 << " " << a1 << " " << a2 << "\n";
566 os << sigma << "\n";
568 os.precision(pr);
569 return os;
570#endif
571}
572
573std::istream & RandPoissonQ::get ( std::istream & is ) {
574 std::string inName;
575 is >> inName;
576 if (inName != name()) {
577 is.clear(std::ios::badbit | is.rdstate());
578 std::cerr << "Mismatch when expecting to read state of a "
579 << name() << " distribution\n"
580 << "Name found was " << inName
581 << "\nistream is left in the badbit state\n";
582 return is;
583 }
584 if (possibleKeywordInput(is, "Uvec", a0)) {
585 std::vector<unsigned long> t(2);
586 is >> a0 >> t[0] >> t[1]; a0 = DoubConv::longs2double(t);
587 is >> a1 >> t[0] >> t[1]; a1 = DoubConv::longs2double(t);
588 is >> a2 >> t[0] >> t[1]; a2 = DoubConv::longs2double(t);
589 is >> sigma >> t[0] >> t[1]; sigma = DoubConv::longs2double(t);
591 return is;
592 }
593 // is >> a0 encompassed by possibleKeywordInput
594 is >> a1 >> a2 >> sigma;
596 return is;
597}
598
599} // namespace CLHEP
600
static double longs2double(const std::vector< unsigned long > &v)
Definition DoubConv.cc:106
static std::vector< unsigned long > dto2longs(double d)
Definition DoubConv.cc:90
static HepRandomEngine * getTheEngine()
Definition Random.cc:166
static double shoot()
static long shoot(double m=1.0)
std::ostream & put(std::ostream &os) const
void fireArray(const int size, long *vect)
static void shootArray(const int size, long *vect, double m=1.0)
std::istream & get(std::istream &is)
HepRandomEngine & engine()
std::string name() const
std::ostream & put(std::ostream &os) const
HepRandomEngine * getLocalEngine()
static long shoot(double m=1.0)
std::istream & get(std::istream &is)
HepRandomEngine & engine()
bool possibleKeywordInput(IS &is, const std::string &key, T &t)