CLHEP 2.0.4.7 Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // $Id: RandBinomial.cc,v 1.3.4.2 2005/04/15 16:32:53 garren Exp $ 00002 // -*- C++ -*- 00003 // 00004 // ----------------------------------------------------------------------- 00005 // HEP Random 00006 // --- RandBinomial --- 00007 // class implementation file 00008 // ----------------------------------------------------------------------- 00009 00010 // ======================================================================= 00011 // John Marraffino - Created: 12th May 1998 00012 // M Fischler - put and get to/from streams 12/10/04 00013 // M Fischler - put/get to/from streams uses pairs of ulongs when 00014 // + storing doubles avoid problems with precision 00015 // 4/14/05 00016 // 00017 // ======================================================================= 00018 00019 #include "CLHEP/Random/RandBinomial.h" 00020 #include "CLHEP/Random/defs.h" 00021 #include "CLHEP/Random/DoubConv.hh" 00022 #include <algorithm> // for min() and max() 00023 #include <cmath> // for exp() 00024 00025 using namespace std; 00026 00027 namespace CLHEP { 00028 00029 std::string RandBinomial::name() const {return "RandBinomial";} 00030 HepRandomEngine & RandBinomial::engine() {return *localEngine;} 00031 00032 RandBinomial::~RandBinomial() { 00033 if ( deleteEngine ) delete localEngine; 00034 } 00035 00036 RandBinomial::RandBinomial(const RandBinomial& right) 00037 : defaultN(right.defaultN), defaultP(right.defaultP) 00038 {;} 00039 00040 double RandBinomial::shoot( HepRandomEngine *anEngine, long n, 00041 double p ) { 00042 return genBinomial( anEngine, n, p ); 00043 } 00044 00045 double RandBinomial::shoot( long n, double p ) { 00046 HepRandomEngine *anEngine = HepRandom::getTheEngine(); 00047 return genBinomial( anEngine, n, p ); 00048 } 00049 00050 double RandBinomial::fire( long n, double p ) { 00051 return genBinomial( localEngine, n, p ); 00052 } 00053 00054 void RandBinomial::shootArray( const int size, double* vect, 00055 long n, double p ) 00056 { 00057 int i; 00058 00059 for (i=0; i<size; ++i) 00060 vect[i] = shoot(n,p); 00061 } 00062 00063 void RandBinomial::shootArray( HepRandomEngine* anEngine, 00064 const int size, double* vect, 00065 long n, double p ) 00066 { 00067 int i; 00068 00069 for (i=0; i<size; ++i) 00070 vect[i] = shoot(anEngine,n,p); 00071 } 00072 00073 void RandBinomial::fireArray( const int size, double* vect) 00074 { 00075 int i; 00076 00077 for (i=0; i<size; ++i) 00078 vect[i] = fire(defaultN,defaultP); 00079 } 00080 00081 void RandBinomial::fireArray( const int size, double* vect, 00082 long n, double p ) 00083 { 00084 int i; 00085 00086 for (i=0; i<size; ++i) 00087 vect[i] = fire(n,p); 00088 } 00089 00090 /************************************************************************* 00091 * * 00092 * StirlingCorrection() * 00093 * * 00094 * Correction term of the Stirling approximation for log(k!) * 00095 * (series in 1/k, or table values for small k) * 00096 * with long int parameter k * 00097 * * 00098 ************************************************************************* 00099 * * 00100 * log k! = (k + 1/2)log(k + 1) - (k + 1) + (1/2)log(2Pi) + * 00101 * StirlingCorrection(k + 1) * 00102 * * 00103 * log k! = (k + 1/2)log(k) - k + (1/2)log(2Pi) + * 00104 * StirlingCorrection(k) * 00105 * * 00106 *************************************************************************/ 00107 00108 static double StirlingCorrection(long int k) 00109 { 00110 #define C1 8.33333333333333333e-02 // +1/12 00111 #define C3 -2.77777777777777778e-03 // -1/360 00112 #define C5 7.93650793650793651e-04 // +1/1260 00113 #define C7 -5.95238095238095238e-04 // -1/1680 00114 00115 static double c[31] = { 0.0, 00116 8.106146679532726e-02, 4.134069595540929e-02, 00117 2.767792568499834e-02, 2.079067210376509e-02, 00118 1.664469118982119e-02, 1.387612882307075e-02, 00119 1.189670994589177e-02, 1.041126526197209e-02, 00120 9.255462182712733e-03, 8.330563433362871e-03, 00121 7.573675487951841e-03, 6.942840107209530e-03, 00122 6.408994188004207e-03, 5.951370112758848e-03, 00123 5.554733551962801e-03, 5.207655919609640e-03, 00124 4.901395948434738e-03, 4.629153749334029e-03, 00125 4.385560249232324e-03, 4.166319691996922e-03, 00126 3.967954218640860e-03, 3.787618068444430e-03, 00127 3.622960224683090e-03, 3.472021382978770e-03, 00128 3.333155636728090e-03, 3.204970228055040e-03, 00129 3.086278682608780e-03, 2.976063983550410e-03, 00130 2.873449362352470e-03, 2.777674929752690e-03, 00131 }; 00132 double r, rr; 00133 00134 if (k > 30L) { 00135 r = 1.0 / (double) k; rr = r * r; 00136 return(r*(C1 + rr*(C3 + rr*(C5 + rr*C7)))); 00137 } 00138 else return(c[k]); 00139 } 00140 00141 double RandBinomial::genBinomial( HepRandomEngine *anEngine, long n, double p ) 00142 { 00143 /****************************************************************** 00144 * * 00145 * Binomial-Distribution - Acceptance Rejection/Inversion * 00146 * * 00147 ****************************************************************** 00148 * * 00149 * Acceptance Rejection method combined with Inversion for * 00150 * generating Binomial random numbers with parameters * 00151 * n (number of trials) and p (probability of success). * 00152 * For min(n*p,n*(1-p)) < 10 the Inversion method is applied: * 00153 * The random numbers are generated via sequential search, * 00154 * starting at the lowest index k=0. The cumulative probabilities * 00155 * are avoided by using the technique of chop-down. * 00156 * For min(n*p,n*(1-p)) >= 10 Acceptance Rejection is used: * 00157 * The algorithm is based on a hat-function which is uniform in * 00158 * the centre region and exponential in the tails. * 00159 * A triangular immediate acceptance region in the centre speeds * 00160 * up the generation of binomial variates. * 00161 * If candidate k is near the mode, f(k) is computed recursively * 00162 * starting at the mode m. * 00163 * The acceptance test by Stirling's formula is modified * 00164 * according to W. Hoermann (1992): The generation of binomial * 00165 * random variates, to appear in J. Statist. Comput. Simul. * 00166 * If p < .5 the algorithm is applied to parameters n, p. * 00167 * Otherwise p is replaced by 1-p, and k is replaced by n - k. * 00168 * * 00169 ****************************************************************** 00170 * * 00171 * FUNCTION: - btpec samples a random number from the binomial * 00172 * distribution with parameters n and p and is * 00173 * valid for n*min(p,1-p) > 0. * 00174 * REFERENCE: - V. Kachitvichyanukul, B.W. Schmeiser (1988): * 00175 * Binomial random variate generation, * 00176 * Communications of the ACM 31, 216-222. * 00177 * SUBPROGRAMS: - StirlingCorrection() * 00178 * ... Correction term of the Stirling * 00179 * approximation for log(k!) * 00180 * (series in 1/k or table values * 00181 * for small k) with long int k * 00182 * - anEngine ... Pointer to a (0,1)-Uniform * 00183 * engine * 00184 * * 00185 * Implemented by H. Zechner and P. Busswald, September 1992 * 00186 ******************************************************************/ 00187 00188 #define C1_3 0.33333333333333333 00189 #define C5_8 0.62500000000000000 00190 #define C1_6 0.16666666666666667 00191 #define DMAX_KM 20L 00192 00193 static long int n_last = -1L, n_prev = -1L; 00194 static double par,np,p0,q,p_last = -1.0, p_prev = -1.0; 00195 static long b,m,nm; 00196 static double pq, rc, ss, xm, xl, xr, ll, lr, c, 00197 p1, p2, p3, p4, ch; 00198 00199 long bh,i, K, Km, nK; 00200 double f, rm, U, V, X, T, E; 00201 00202 if (n != n_last || p != p_last) // set-up 00203 { 00204 n_last = n; 00205 p_last = p; 00206 par=min(p,1.0-p); 00207 q=1.0-par; 00208 np = n*par; 00209 00210 // Check for invalid input values 00211 00212 if( np <= 0.0 ) return (-1.0); 00213 00214 rm = np + par; 00215 m = (long int) rm; // mode, integer 00216 if (np<10) 00217 { 00218 p0=exp(n*log(q)); // Chop-down 00219 bh=(long int)(np+10.0*sqrt(np*q)); 00220 b=min(n,bh); 00221 } 00222 else 00223 { 00224 rc = (n + 1.0) * (pq = par / q); // recurr. relat. 00225 ss = np * q; // variance 00226 i = (long int) (2.195*sqrt(ss) - 4.6*q); // i = p1 - 0.5 00227 xm = m + 0.5; 00228 xl = (double) (m - i); // limit left 00229 xr = (double) (m + i + 1L); // limit right 00230 f = (rm - xl) / (rm - xl*par); ll = f * (1.0 + 0.5*f); 00231 f = (xr - rm) / (xr * q); lr = f * (1.0 + 0.5*f); 00232 c = 0.134 + 20.5/(15.3 + (double) m); // parallelogram 00233 // height 00234 p1 = i + 0.5; 00235 p2 = p1 * (1.0 + c + c); // probabilities 00236 p3 = p2 + c/ll; // of regions 1-4 00237 p4 = p3 + c/lr; 00238 } 00239 } 00240 if (np<10) //Inversion Chop-down 00241 { 00242 double pk; 00243 00244 K=0; 00245 pk=p0; 00246 U=anEngine->flat(); 00247 while (U>pk) 00248 { 00249 ++K; 00250 if (K>b) 00251 { 00252 U=anEngine->flat(); 00253 K=0; 00254 pk=p0; 00255 } 00256 else 00257 { 00258 U-=pk; 00259 pk=(double)(((n-K+1)*par*pk)/(K*q)); 00260 } 00261 } 00262 return ((p>0.5) ? (double)(n-K):(double)K); 00263 } 00264 00265 for (;;) 00266 { 00267 V = anEngine->flat(); 00268 if ((U = anEngine->flat() * p4) <= p1) // triangular region 00269 { 00270 K=(long int) (xm - U + p1*V); 00271 return ((p>0.5) ? (double)(n-K):(double)K); // immediate accept 00272 } 00273 if (U <= p2) // parallelogram 00274 { 00275 X = xl + (U - p1)/c; 00276 if ((V = V*c + 1.0 - fabs(xm - X)/p1) >= 1.0) continue; 00277 K = (long int) X; 00278 } 00279 else if (U <= p3) // left tail 00280 { 00281 if ((X = xl + log(V)/ll) < 0.0) continue; 00282 K = (long int) X; 00283 V *= (U - p2) * ll; 00284 } 00285 else // right tail 00286 { 00287 if ((K = (long int) (xr - log(V)/lr)) > n) continue; 00288 V *= (U - p3) * lr; 00289 } 00290 00291 // acceptance test : two cases, depending on |K - m| 00292 if ((Km = labs(K - m)) <= DMAX_KM || Km + Km + 2L >= ss) 00293 { 00294 00295 // computation of p(K) via recurrence relationship from the mode 00296 f = 1.0; // f(m) 00297 if (m < K) 00298 { 00299 for (i = m; i < K; ) 00300 { 00301 if ((f *= (rc / ++i - pq)) < V) break; // multiply f 00302 } 00303 } 00304 else 00305 { 00306 for (i = K; i < m; ) 00307 { 00308 if ((V *= (rc / ++i - pq)) > f) break; // multiply V 00309 } 00310 } 00311 if (V <= f) break; // acceptance test 00312 } 00313 else 00314 { 00315 00316 // lower and upper squeeze tests, based on lower bounds for log p(K) 00317 V = log(V); 00318 T = - Km * Km / (ss + ss); 00319 E = (Km / ss) * ((Km * (Km * C1_3 + C5_8) + C1_6) / ss + 0.5); 00320 if (V <= T - E) break; 00321 if (V <= T + E) 00322 { 00323 if (n != n_prev || par != p_prev) 00324 { 00325 n_prev = n; 00326 p_prev = par; 00327 00328 nm = n - m + 1L; 00329 ch = xm * log((m + 1.0)/(pq * nm)) + 00330 StirlingCorrection(m + 1L) + StirlingCorrection(nm); 00331 } 00332 nK = n - K + 1L; 00333 00334 // computation of log f(K) via Stirling's formula 00335 // final acceptance-rejection test 00336 if (V <= ch + (n + 1.0)*log((double) nm / (double) nK) + 00337 (K + 0.5)*log(nK * pq / (K + 1.0)) - 00338 StirlingCorrection(K + 1L) - StirlingCorrection(nK)) break; 00339 } 00340 } 00341 } 00342 return ((p>0.5) ? (double)(n-K):(double)K); 00343 } 00344 00345 std::ostream & RandBinomial::put ( std::ostream & os ) const { 00346 int pr=os.precision(20); 00347 std::vector<unsigned long> t(2); 00348 os << " " << name() << "\n"; 00349 os << "Uvec" << "\n"; 00350 t = DoubConv::dto2longs(defaultP); 00351 os << defaultN << " " << defaultP << " " << t[0] << " " << t[1] << "\n"; 00352 os.precision(pr); 00353 return os; 00354 #ifdef REMOVED 00355 int pr=os.precision(20); 00356 os << " " << name() << "\n"; 00357 os << defaultN << " " << defaultP << "\n"; 00358 os.precision(pr); 00359 return os; 00360 #endif 00361 } 00362 00363 std::istream & RandBinomial::get ( std::istream & is ) { 00364 std::string inName; 00365 is >> inName; 00366 if (inName != name()) { 00367 is.clear(std::ios::badbit | is.rdstate()); 00368 std::cerr << "Mismatch when expecting to read state of a " 00369 << name() << " distribution\n" 00370 << "Name found was " << inName 00371 << "\nistream is left in the badbit state\n"; 00372 return is; 00373 } 00374 if (possibleKeywordInput(is, "Uvec", defaultN)) { 00375 std::vector<unsigned long> t(2); 00376 is >> defaultN >> defaultP; 00377 is >> t[0] >> t[1]; defaultP = DoubConv::longs2double(t); 00378 return is; 00379 } 00380 // is >> defaultN encompassed by possibleKeywordInput 00381 is >> defaultP; 00382 return is; 00383 } 00384 00385 00386 } // namespace CLHEP