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