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 
19 #include "CLHEP/Random/RandBinomial.h"
20 #include "CLHEP/Random/defs.h"
21 #include "CLHEP/Random/DoubConv.hh"
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 
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 ) {
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 
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 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 
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 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 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 
330 std::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 
348 std::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
X
Definition: testSharedPtrBasic.cc:28
double
#define double(obj)
Definition: excDblThrow.cc:32
CLHEP::HepRandomEngine
Definition: Matrix/CLHEP/Random/RandomEngine.h:55
b
@ b
Definition: testCategories.cc:125
C3
#define C3
CLHEP::RandBinomial::shoot
static double shoot()
is
HepRotation and so forth isNear() norm2() rectify() static Rotation row1 row4(To avoid bloat in the code pulled in for programs which don 't use all these features, we split the implementation .cc files. Only isNear() goes into the original Rotation.cc) --------------------------------------- HepAxisAngle and HepEulerAngles classes --------------------------------------- These classes are very useful and simple structures for holding the result of a nice intuituve decomposition of a rotation there is no longer much content in the distinct ZOOM PhysicsVectors library The only content left in the library is the object files representing the various Exception objects When we build the CLHEP classes for the ZOOM we will set up so as to use ZOOM SpaceVector is(but we can disable namespace usage and most of our users do so at this point). What I do is leave Hep3Vector in the global namespace
C7
#define C7
C5
#define C5
CLHEP::RandBinomial::engine
HepRandomEngine & engine()
Definition: RandBinomial.cc:28
CLHEP::DoubConv::longs2double
static double longs2double(const std::vector< unsigned long > &v)
Definition: DoubConv.cc:106
C1_3
#define C1_3
CLHEP::RandBinomial::fire
double fire()
C1
#define C1
CLHEP::detail::n
n
Definition: Ranlux64Engine.cc:85
size
user code seldom needs to call this function directly ZMerrno whether or not they are still recorded ZMerrno size() Return the(integer) number of ZMthrow 'n exceptions currently recorded. 5) ZMerrno.clear() Set an internal counter to zero. This counter is available(see next function) to user code to track ZMthrow 'n exceptions that have occurred during any arbitrary time interval. 6) ZMerrno.countSinceCleared() Return the(integer) number of ZMthrow 'n exceptions that have been recorded via ZMerrno.write()
CLHEP::RandBinomial::name
std::string name() const
Definition: RandBinomial.cc:27
f
void f(void g())
Definition: excDblThrow.cc:38
CLHEP
Definition: ClhepVersion.h:13
v
they are gone ZOOM Features Discontinued The following features of the ZOOM package were felt to be extreme overkill These have been after checking that no existing user code was utilizing as in SpaceVector v
Definition: keyMergeIssues.doc:324
CLHEP::RandBinomial::shootArray
static void shootArray(const int size, double *vect, long n=1, double p=0.5)
Definition: RandBinomial.cc:47
n_constructors::p0
incomplete * p0
Definition: testSharedPtr.cc:393
C1_6
#define C1_6
CLHEP::RandBinomial::get
std::istream & get(std::istream &is)
Definition: RandBinomial.cc:348
CLHEP::RandBinomial::fireArray
void fireArray(const int size, double *vect)
Definition: RandBinomial.cc:62
DMAX_KM
#define DMAX_KM
C5_8
#define C5_8
CLHEP::possibleKeywordInput
bool possibleKeywordInput(IS &is, const std::string &key, T &t)
Definition: Matrix/CLHEP/Random/RandomEngine.h:168
CLHEP::RandBinomial::put
std::ostream & put(std::ostream &os) const
Definition: RandBinomial.cc:330
i
long i
Definition: JamesRandomSeeding.txt:27
CLHEP::DoubConv::dto2longs
static std::vector< unsigned long > dto2longs(double d)
Definition: DoubConv.cc:90
CLHEP::HepRandom::getTheEngine
static HepRandomEngine * getTheEngine()
Definition: Random.cc:166
k
long k
Definition: JamesRandomSeeding.txt:29
n_constructors::m
int m
Definition: testSharedPtr.cc:370
CLHEP::RandBinomial::~RandBinomial
virtual ~RandBinomial()
Definition: RandBinomial.cc:30