CLHEP VERSION Reference Documentation
   
CLHEP Home Page     CLHEP Documentation     CLHEP Bug Reports

RandBinomial.cc
Go to the documentation of this file.
1// $Id: RandBinomial.cc,v 1.5 2010/06/16 17:24:53 garren Exp $
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
20#include "CLHEP/Random/defs.h"
22#include <algorithm> // for min() and max()
23#include <cmath> // for exp()
24
25namespace CLHEP {
26
27std::string RandBinomial::name() const {return "RandBinomial";}
28HepRandomEngine & RandBinomial::engine() {return *localEngine;}
29
32
33double RandBinomial::shoot( HepRandomEngine *anEngine, long n,
34 double p ) {
35 return genBinomial( anEngine, n, p );
36}
37
38double RandBinomial::shoot( long n, double p ) {
40 return genBinomial( anEngine, n, p );
41}
42
43double RandBinomial::fire( long n, double p ) {
44 return genBinomial( localEngine.get(), n, p );
45}
46
47void 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
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
62void RandBinomial::fireArray( const int size, double* vect)
63{
64 for( double* v = vect; v != vect+size; ++v )
65 *v = fire(defaultN,defaultP);
66}
67
68void 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 std::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
93static 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 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
126double 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 std::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 long int n_last = -1L, n_prev = -1L;
179 static double par,np,p0,q,p_last = -1.0, p_prev = -1.0;
180 static long b,m,nm;
181 static 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<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 = 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)) +
315 StirlingCorrection(m + 1L) + StirlingCorrection(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
330std::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";
335 t = DoubConv::dto2longs(defaultP);
336 os << defaultN << " " << defaultP << " " << t[0] << " " << t[1] << "\n";
337 os.precision(pr);
338 return os;
339#ifdef REMOVED
340 int pr=os.precision(20);
341 os << " " << name() << "\n";
342 os << defaultN << " " << defaultP << "\n";
343 os.precision(pr);
344 return os;
345#endif
346}
347
348std::istream & RandBinomial::get ( std::istream & is ) {
349 std::string inName;
350 is >> inName;
351 if (inName != name()) {
352 is.clear(std::ios::badbit | is.rdstate());
353 std::cerr << "Mismatch when expecting to read state of a "
354 << name() << " distribution\n"
355 << "Name found was " << inName
356 << "\nistream is left in the badbit state\n";
357 return is;
358 }
359 if (possibleKeywordInput(is, "Uvec", defaultN)) {
360 std::vector<unsigned long> t(2);
361 is >> defaultN >> defaultP;
362 is >> t[0] >> t[1]; defaultP = DoubConv::longs2double(t);
363 return is;
364 }
365 // is >> defaultN encompassed by possibleKeywordInput
366 is >> defaultP;
367 return is;
368}
369
370
371} // 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:106
static std::vector< unsigned long > dto2longs(double d)
Definition DoubConv.cc:90
static HepRandomEngine * getTheEngine()
Definition Random.cc:166
std::string name() const
HepRandomEngine & engine()
static double shoot()
void fireArray(const int size, double *vect)
static void shootArray(const int size, double *vect, long n=1, double p=0.5)
std::ostream & put(std::ostream &os) const
std::istream & get(std::istream &is)
#define double(obj)
void f(void g())
bool possibleKeywordInput(IS &is, const std::string &key, T &t)
incomplete * p0