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