Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
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
21#include <algorithm> // for min() and max()
22#include <cmath> // for exp()
23
24namespace CLHEP {
25
26std::string RandBinomial::name() const {return "RandBinomial";}
27HepRandomEngine & RandBinomial::engine() {return *localEngine;}
28
30}
31
32double RandBinomial::shoot( HepRandomEngine *anEngine, long n,
33 double p ) {
34 return genBinomial( anEngine, n, p );
35}
36
37double RandBinomial::shoot( long n, double p ) {
39 return genBinomial( anEngine, n, p );
40}
41
42double RandBinomial::fire( long n, double p ) {
43 return genBinomial( localEngine.get(), n, p );
44}
45
46void 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
61void RandBinomial::fireArray( const int size, double* vect)
62{
63 for( double* v = vect; v != vect+size; ++v )
64 *v = fire(defaultN,defaultP);
65}
66
67void 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
92static 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 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
125double 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 long int n_last = -1L, n_prev = -1L;
178 static double par,np,p0,q,p_last = -1.0, p_prev = -1.0;
179 static long b,m,nm;
180 static 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<10) //Inversion Chop-down
225 {
226 double pk;
227
228 K=0;
229 pk=p0;
230 U=anEngine->flat();
231 while (U>pk)
232 {
233 ++K;
234 if (K>b)
235 {
236 U=anEngine->flat();
237 K=0;
238 pk=p0;
239 }
240 else
241 {
242 U-=pk;
243 pk=(double)(((n-K+1)*par*pk)/(K*q));
244 }
245 }
246 return ((p>0.5) ? (double)(n-K):(double)K);
247 }
248
249 for (;;)
250 {
251 V = anEngine->flat();
252 if ((U = anEngine->flat() * p4) <= p1) // triangular region
253 {
254 K=(long int) (xm - U + p1*V);
255 return ((p>0.5) ? (double)(n-K):(double)K); // immediate accept
256 }
257 if (U <= p2) // parallelogram
258 {
259 X = xl + (U - p1)/c;
260 if ((V = V*c + 1.0 - std::fabs(xm - X)/p1) >= 1.0) continue;
261 K = (long int) X;
262 }
263 else if (U <= p3) // left tail
264 {
265 if ((X = xl + std::log(V)/ll) < 0.0) continue;
266 K = (long int) X;
267 V *= (U - p2) * ll;
268 }
269 else // right tail
270 {
271 if ((K = (long int) (xr - std::log(V)/lr)) > n) continue;
272 V *= (U - p3) * lr;
273 }
274
275 // acceptance test : two cases, depending on |K - m|
276 if ((Km = std::labs(K - m)) <= DMAX_KM || Km + Km + 2L >= ss)
277 {
278
279 // computation of p(K) via recurrence relationship from the mode
280 f = 1.0; // f(m)
281 if (m < K)
282 {
283 for (i = m; i < K; )
284 {
285 if ((f *= (rc / ++i - pq)) < V) break; // multiply f
286 }
287 }
288 else
289 {
290 for (i = K; i < m; )
291 {
292 if ((V *= (rc / ++i - pq)) > f) break; // multiply V
293 }
294 }
295 if (V <= f) break; // acceptance test
296 }
297 else
298 {
299
300 // lower and upper squeeze tests, based on lower bounds for log p(K)
301 V = std::log(V);
302 T = - Km * Km / (ss + ss);
303 E = (Km / ss) * ((Km * (Km * C1_3 + C5_8) + C1_6) / ss + 0.5);
304 if (V <= T - E) break;
305 if (V <= T + E)
306 {
307 if (n != n_prev || par != p_prev)
308 {
309 n_prev = n;
310 p_prev = par;
311
312 nm = n - m + 1L;
313 ch = xm * std::log((m + 1.0)/(pq * nm)) +
314 StirlingCorrection(m + 1L) + StirlingCorrection(nm);
315 }
316 nK = n - K + 1L;
317
318 // computation of log f(K) via Stirling's formula
319 // final acceptance-rejection test
320 if (V <= ch + (n + 1.0)*std::log((double) nm / (double) nK) +
321 (K + 0.5)*std::log(nK * pq / (K + 1.0)) -
322 StirlingCorrection(K + 1L) - StirlingCorrection(nK)) break;
323 }
324 }
325 }
326 return ((p>0.5) ? (double)(n-K):(double)K);
327}
328
329std::ostream & RandBinomial::put ( std::ostream & os ) const {
330 int pr=os.precision(20);
331 std::vector<unsigned long> t(2);
332 os << " " << name() << "\n";
333 os << "Uvec" << "\n";
334 t = DoubConv::dto2longs(defaultP);
335 os << defaultN << " " << defaultP << " " << t[0] << " " << t[1] << "\n";
336 os.precision(pr);
337 return os;
338}
339
340std::istream & RandBinomial::get ( std::istream & is ) {
341 std::string inName;
342 is >> inName;
343 if (inName != name()) {
344 is.clear(std::ios::badbit | is.rdstate());
345 std::cerr << "Mismatch when expecting to read state of a "
346 << name() << " distribution\n"
347 << "Name found was " << inName
348 << "\nistream is left in the badbit state\n";
349 return is;
350 }
351 if (possibleKeywordInput(is, "Uvec", defaultN)) {
352 std::vector<unsigned long> t(2);
353 is >> defaultN >> defaultP;
354 is >> t[0] >> t[1]; defaultP = DoubConv::longs2double(t);
355 return is;
356 }
357 // is >> defaultN encompassed by possibleKeywordInput
358 is >> defaultP;
359 return is;
360}
361
362
363} // namespace CLHEP
#define C5_8
#define C5
#define C1
#define C1_6
#define C3
#define C1_3
#define DMAX_KM
#define C7
static double longs2double(const std::vector< unsigned long > &v)
Definition: DoubConv.cc:114
static std::vector< unsigned long > dto2longs(double d)
Definition: DoubConv.cc:98
static HepRandomEngine * getTheEngine()
Definition: Random.cc:165
std::string name() const
Definition: RandBinomial.cc:26
HepRandomEngine & engine()
Definition: RandBinomial.cc:27
static double shoot()
void fireArray(const int size, double *vect)
Definition: RandBinomial.cc:61
static void shootArray(const int size, double *vect, long n=1, double p=0.5)
Definition: RandBinomial.cc:46
virtual ~RandBinomial()
Definition: RandBinomial.cc:29
std::ostream & put(std::ostream &os) const
std::istream & get(std::istream &is)
Definition: DoubConv.h:17
bool possibleKeywordInput(IS &is, const std::string &key, T &t)
Definition: RandomEngine.h:167