ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RandBinomial.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file RandBinomial.cc
1 // -*- C++ -*-
2 //
3 // -----------------------------------------------------------------------
4 // HEP Random
5 // --- RandBinomial ---
6 // class implementation file
7 // -----------------------------------------------------------------------
8 
9 // =======================================================================
10 // John Marraffino - Created: 12th May 1998
11 // M Fischler - put and get to/from streams 12/10/04
12 // M Fischler - put/get to/from streams uses pairs of ulongs when
13 // + storing doubles avoid problems with precision
14 // 4/14/05
15 //
16 // =======================================================================
17 
19 #include "CLHEP/Random/DoubConv.h"
21 #include <algorithm> // for min() and max()
22 #include <cmath> // for exp()
23 
24 namespace CLHEP {
25 
26 std::string RandBinomial::name() const {return "RandBinomial";}
28 
30 }
31 
32 double RandBinomial::shoot( HepRandomEngine *anEngine, long n,
33  double p ) {
34  return genBinomial( anEngine, n, p );
35 }
36 
37 double RandBinomial::shoot( long n, double p ) {
39  return genBinomial( anEngine, n, p );
40 }
41 
42 double RandBinomial::fire( long n, double p ) {
43  return genBinomial( localEngine.get(), n, p );
44 }
45 
46 void RandBinomial::shootArray( const int size, double* vect,
47  long n, double p )
48 {
49  for( double* v = vect; v != vect+size; ++v )
50  *v = shoot(n,p);
51 }
52 
54  const int size, double* vect,
55  long n, double p )
56 {
57  for( double* v = vect; v != vect+size; ++v )
58  *v = shoot(anEngine,n,p);
59 }
60 
61 void RandBinomial::fireArray( const int size, double* vect)
62 {
63  for( double* v = vect; v != vect+size; ++v )
64  *v = fire(defaultN,defaultP);
65 }
66 
67 void RandBinomial::fireArray( const int size, double* vect,
68  long n, double p )
69 {
70  for( double* v = vect; v != vect+size; ++v )
71  *v = fire(n,p);
72 }
73 
74 /*************************************************************************
75  * *
76  * StirlingCorrection() *
77  * *
78  * Correction term of the Stirling approximation for log(k!) *
79  * (series in 1/k, or table values for small k) *
80  * with long int parameter k *
81  * *
82  *************************************************************************
83  * *
84  * log k! = (k + 1/2)log(k + 1) - (k + 1) + (1/2)log(2Pi) + *
85  * StirlingCorrection(k + 1) *
86  * *
87  * log k! = (k + 1/2)log(k) - k + (1/2)log(2Pi) + *
88  * StirlingCorrection(k) *
89  * *
90  *************************************************************************/
91 
92 static double StirlingCorrection(long int k)
93 {
94  #define C1 8.33333333333333333e-02 // +1/12
95  #define C3 -2.77777777777777778e-03 // -1/360
96  #define C5 7.93650793650793651e-04 // +1/1260
97  #define C7 -5.95238095238095238e-04 // -1/1680
98 
99  static const double c[31] = { 0.0,
100  8.106146679532726e-02, 4.134069595540929e-02,
101  2.767792568499834e-02, 2.079067210376509e-02,
102  1.664469118982119e-02, 1.387612882307075e-02,
103  1.189670994589177e-02, 1.041126526197209e-02,
104  9.255462182712733e-03, 8.330563433362871e-03,
105  7.573675487951841e-03, 6.942840107209530e-03,
106  6.408994188004207e-03, 5.951370112758848e-03,
107  5.554733551962801e-03, 5.207655919609640e-03,
108  4.901395948434738e-03, 4.629153749334029e-03,
109  4.385560249232324e-03, 4.166319691996922e-03,
110  3.967954218640860e-03, 3.787618068444430e-03,
111  3.622960224683090e-03, 3.472021382978770e-03,
112  3.333155636728090e-03, 3.204970228055040e-03,
113  3.086278682608780e-03, 2.976063983550410e-03,
114  2.873449362352470e-03, 2.777674929752690e-03,
115  };
116  double r, rr;
117 
118  if (k > 30L) {
119  r = 1.0 / (double) k; rr = r * r;
120  return(r*(C1 + rr*(C3 + rr*(C5 + rr*C7))));
121  }
122  else return(c[k]);
123 }
124 
125 double RandBinomial::genBinomial( HepRandomEngine *anEngine, long n, double p )
126 {
127 /******************************************************************
128  * *
129  * Binomial-Distribution - Acceptance Rejection/Inversion *
130  * *
131  ******************************************************************
132  * *
133  * Acceptance Rejection method combined with Inversion for *
134  * generating Binomial random numbers with parameters *
135  * n (number of trials) and p (probability of success). *
136  * For min(n*p,n*(1-p)) < 10 the Inversion method is applied: *
137  * The random numbers are generated via sequential search, *
138  * starting at the lowest index k=0. The cumulative probabilities *
139  * are avoided by using the technique of chop-down. *
140  * For min(n*p,n*(1-p)) >= 10 Acceptance Rejection is used: *
141  * The algorithm is based on a hat-function which is uniform in *
142  * the centre region and exponential in the tails. *
143  * A triangular immediate acceptance region in the centre speeds *
144  * up the generation of binomial variates. *
145  * If candidate k is near the mode, f(k) is computed recursively *
146  * starting at the mode m. *
147  * The acceptance test by Stirling's formula is modified *
148  * according to W. Hoermann (1992): The generation of binomial *
149  * random variates, to appear in J. Statist. Comput. Simul. *
150  * If p < .5 the algorithm is applied to parameters n, p. *
151  * Otherwise p is replaced by 1-p, and k is replaced by n - k. *
152  * *
153  ******************************************************************
154  * *
155  * FUNCTION: - btpec samples a random number from the binomial *
156  * distribution with parameters n and p and is *
157  * valid for n*min(p,1-p) > 0. *
158  * REFERENCE: - V. Kachitvichyanukul, B.W. Schmeiser (1988): *
159  * Binomial random variate generation, *
160  * Communications of the ACM 31, 216-222. *
161  * SUBPROGRAMS: - StirlingCorrection() *
162  * ... Correction term of the Stirling *
163  * approximation for log(k!) *
164  * (series in 1/k or table values *
165  * for small k) with long int k *
166  * - anEngine ... Pointer to a (0,1)-Uniform *
167  * engine *
168  * *
169  * Implemented by H. Zechner and P. Busswald, September 1992 *
170  ******************************************************************/
171 
172 #define C1_3 0.33333333333333333
173 #define C5_8 0.62500000000000000
174 #define C1_6 0.16666666666666667
175 #define DMAX_KM 20L
176 
177  static CLHEP_THREAD_LOCAL long int n_last = -1L, n_prev = -1L;
178  static CLHEP_THREAD_LOCAL double par,np,p0,q,p_last = -1.0, p_prev = -1.0;
179  static CLHEP_THREAD_LOCAL long b,m,nm;
180  static CLHEP_THREAD_LOCAL double pq, rc, ss, xm, xl, xr, ll, lr, c,
181  p1, p2, p3, p4, ch;
182 
183  long bh,i, K, Km, nK;
184  double f, rm, U, V, X, T, E;
185 
186  if (n != n_last || p != p_last) // set-up
187  {
188  n_last = n;
189  p_last = p;
190  par=std::min(p,1.0-p);
191  q=1.0-par;
192  np = n*par;
193 
194 // Check for invalid input values
195 
196  if( np <= 0.0 ) return (-1.0);
197 
198  rm = np + par;
199  m = (long int) rm; // mode, integer
200  if (np<10)
201  {
202  p0=std::exp(n*std::log(q)); // Chop-down
203  bh=(long int)(np+10.0*std::sqrt(np*q));
204  b=std::min(n,bh);
205  }
206  else
207  {
208  rc = (n + 1.0) * (pq = par / q); // recurr. relat.
209  ss = np * q; // variance
210  i = (long int) (2.195*std::sqrt(ss) - 4.6*q); // i = p1 - 0.5
211  xm = m + 0.5;
212  xl = (double) (m - i); // limit left
213  xr = (double) (m + i + 1L); // limit right
214  f = (rm - xl) / (rm - xl*par); ll = f * (1.0 + 0.5*f);
215  f = (xr - rm) / (xr * q); lr = f * (1.0 + 0.5*f);
216  c = 0.134 + 20.5/(15.3 + (double) m); // parallelogram
217  // height
218  p1 = i + 0.5;
219  p2 = p1 * (1.0 + c + c); // probabilities
220  p3 = p2 + c/ll; // of regions 1-4
221  p4 = p3 + c/lr;
222  }
223  }
224  if( np <= 0.0 ) return (-1.0);
225  if (np<10) //Inversion Chop-down
226  {
227  double pk;
228 
229  K=0;
230  pk=p0;
231  U=anEngine->flat();
232  while (U>pk)
233  {
234  ++K;
235  if (K>b)
236  {
237  U=anEngine->flat();
238  K=0;
239  pk=p0;
240  }
241  else
242  {
243  U-=pk;
244  pk=(double)(((n-K+1)*par*pk)/(K*q));
245  }
246  }
247  return ((p>0.5) ? (double)(n-K):(double)K);
248  }
249 
250  for (;;)
251  {
252  V = anEngine->flat();
253  if ((U = anEngine->flat() * p4) <= p1) // triangular region
254  {
255  K=(long int) (xm - U + p1*V);
256  return ((p>0.5) ? (double)(n-K):(double)K); // immediate accept
257  }
258  if (U <= p2) // parallelogram
259  {
260  X = xl + (U - p1)/c;
261  if ((V = V*c + 1.0 - std::fabs(xm - X)/p1) >= 1.0) continue;
262  K = (long int) X;
263  }
264  else if (U <= p3) // left tail
265  {
266  if ((X = xl + std::log(V)/ll) < 0.0) continue;
267  K = (long int) X;
268  V *= (U - p2) * ll;
269  }
270  else // right tail
271  {
272  if ((K = (long int) (xr - std::log(V)/lr)) > n) continue;
273  V *= (U - p3) * lr;
274  }
275 
276  // acceptance test : two cases, depending on |K - m|
277  if ((Km = std::labs(K - m)) <= DMAX_KM || Km + Km + 2L >= ss)
278  {
279 
280  // computation of p(K) via recurrence relationship from the mode
281  f = 1.0; // f(m)
282  if (m < K)
283  {
284  for (i = m; i < K; )
285  {
286  if ((f *= (rc / ++i - pq)) < V) break; // multiply f
287  }
288  }
289  else
290  {
291  for (i = K; i < m; )
292  {
293  if ((V *= (rc / ++i - pq)) > f) break; // multiply V
294  }
295  }
296  if (V <= f) break; // acceptance test
297  }
298  else
299  {
300 
301  // lower and upper squeeze tests, based on lower bounds for log p(K)
302  V = std::log(V);
303  T = - Km * Km / (ss + ss);
304  E = (Km / ss) * ((Km * (Km * C1_3 + C5_8) + C1_6) / ss + 0.5);
305  if (V <= T - E) break;
306  if (V <= T + E)
307  {
308  if (n != n_prev || par != p_prev)
309  {
310  n_prev = n;
311  p_prev = par;
312 
313  nm = n - m + 1L;
314  ch = xm * std::log((m + 1.0)/(pq * nm)) +
316  }
317  nK = n - K + 1L;
318 
319  // computation of log f(K) via Stirling's formula
320  // final acceptance-rejection test
321  if (V <= ch + (n + 1.0)*std::log((double) nm / (double) nK) +
322  (K + 0.5)*std::log(nK * pq / (K + 1.0)) -
323  StirlingCorrection(K + 1L) - StirlingCorrection(nK)) break;
324  }
325  }
326  }
327  return ((p>0.5) ? (double)(n-K):(double)K);
328 }
329 
330 std::ostream & RandBinomial::put ( std::ostream & os ) const {
331  int pr=os.precision(20);
332  std::vector<unsigned long> t(2);
333  os << " " << name() << "\n";
334  os << "Uvec" << "\n";
336  os << defaultN << " " << defaultP << " " << t[0] << " " << t[1] << "\n";
337  os.precision(pr);
338  return os;
339 }
340 
341 std::istream & RandBinomial::get ( std::istream & is ) {
342  std::string inName;
343  is >> inName;
344  if (inName != name()) {
345  is.clear(std::ios::badbit | is.rdstate());
346  std::cerr << "Mismatch when expecting to read state of a "
347  << name() << " distribution\n"
348  << "Name found was " << inName
349  << "\nistream is left in the badbit state\n";
350  return is;
351  }
352  if (possibleKeywordInput(is, "Uvec", defaultN)) {
353  std::vector<unsigned long> t(2);
354  is >> defaultN >> defaultP;
355  is >> t[0] >> t[1]; defaultP = DoubConv::longs2double(t);
356  return is;
357  }
358  // is >> defaultN encompassed by possibleKeywordInput
359  is >> defaultP;
360  return is;
361 }
362 
363 
364 } // namespace CLHEP