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