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