Parent Directory
| 
Revision Log
Revision 1.1 - (view) (download) (as text)
| 1 : | efrank | 1.1 | #include "Python.h" /* for throwing exceptions */ |
| 2 : | #include "Numeric/ranlib.h" | ||
| 3 : | #include <stdio.h> | ||
| 4 : | #include <math.h> | ||
| 5 : | #include <stdlib.h> | ||
| 6 : | #define ABS(x) ((x) >= 0 ? (x) : -(x)) | ||
| 7 : | #ifndef min | ||
| 8 : | #define min(a,b) ((a) <= (b) ? (a) : (b)) | ||
| 9 : | #endif | ||
| 10 : | #ifndef max | ||
| 11 : | #define max(a,b) ((a) >= (b) ? (a) : (b)) | ||
| 12 : | #endif | ||
| 13 : | void ftnstop(char*); | ||
| 14 : | float genbet(float aa,float bb) | ||
| 15 : | |||
| 16 : | #ifdef _MSC_VER /** disable extremly common MSVC warnings **/ | ||
| 17 : | #pragma warning(disable:4305) | ||
| 18 : | #pragma warning(disable:4244) | ||
| 19 : | #endif | ||
| 20 : | |||
| 21 : | |||
| 22 : | /* | ||
| 23 : | ********************************************************************** | ||
| 24 : | float genbet(float aa,float bb) | ||
| 25 : | GeNerate BETa random deviate | ||
| 26 : | Function | ||
| 27 : | Returns a single random deviate from the beta distribution with | ||
| 28 : | parameters A and B. The density of the beta is | ||
| 29 : | x^(a-1) * (1-x)^(b-1) / B(a,b) for 0 < x < 1 | ||
| 30 : | Arguments | ||
| 31 : | aa --> First parameter of the beta distribution | ||
| 32 : | |||
| 33 : | bb --> Second parameter of the beta distribution | ||
| 34 : | |||
| 35 : | Method | ||
| 36 : | R. C. H. Cheng | ||
| 37 : | Generating Beta Variatew with Nonintegral Shape Parameters | ||
| 38 : | Communications of the ACM, 21:317-322 (1978) | ||
| 39 : | (Algorithms BB and BC) | ||
| 40 : | ********************************************************************** | ||
| 41 : | */ | ||
| 42 : | { | ||
| 43 : | #define expmax 89.0 | ||
| 44 : | #define infnty 1.0E38 | ||
| 45 : | static float olda = -1.0; | ||
| 46 : | static float oldb = -1.0; | ||
| 47 : | static float genbet,a,alpha,b,beta,delta,gamma,k1,k2,r,s,t,u1,u2,v,w,y,z; | ||
| 48 : | static long qsame; | ||
| 49 : | |||
| 50 : | qsame = olda == aa && oldb == bb; | ||
| 51 : | if(qsame) goto S20; | ||
| 52 : | if(!(aa <= 0.0 || bb <= 0.0)) goto S10; | ||
| 53 : | fputs(" AA or BB <= 0 in GENBET - Abort!\n",stderr); | ||
| 54 : | fprintf(stderr," AA: %16.6E BB %16.6E\n",aa,bb); | ||
| 55 : | PyErr_SetString (PyExc_ValueError, "Described above."); | ||
| 56 : | return 0.0; | ||
| 57 : | S10: | ||
| 58 : | olda = aa; | ||
| 59 : | oldb = bb; | ||
| 60 : | S20: | ||
| 61 : | if(!(min(aa,bb) > 1.0)) goto S100; | ||
| 62 : | /* | ||
| 63 : | Alborithm BB | ||
| 64 : | Initialize | ||
| 65 : | */ | ||
| 66 : | if(qsame) goto S30; | ||
| 67 : | a = min(aa,bb); | ||
| 68 : | b = max(aa,bb); | ||
| 69 : | alpha = a+b; | ||
| 70 : | beta = sqrt((alpha-2.0)/(2.0*a*b-alpha)); | ||
| 71 : | gamma = a+1.0/beta; | ||
| 72 : | S30: | ||
| 73 : | S40: | ||
| 74 : | u1 = ranf(); | ||
| 75 : | /* | ||
| 76 : | Step 1 | ||
| 77 : | */ | ||
| 78 : | u2 = ranf(); | ||
| 79 : | v = beta*log(u1/(1.0-u1)); | ||
| 80 : | if(!(v > expmax)) goto S50; | ||
| 81 : | w = infnty; | ||
| 82 : | goto S60; | ||
| 83 : | S50: | ||
| 84 : | w = a*exp(v); | ||
| 85 : | S60: | ||
| 86 : | z = pow(u1,2.0)*u2; | ||
| 87 : | r = gamma*v-1.3862944; | ||
| 88 : | s = a+r-w; | ||
| 89 : | /* | ||
| 90 : | Step 2 | ||
| 91 : | */ | ||
| 92 : | if(s+2.609438 >= 5.0*z) goto S70; | ||
| 93 : | /* | ||
| 94 : | Step 3 | ||
| 95 : | */ | ||
| 96 : | t = log(z); | ||
| 97 : | if(s > t) goto S70; | ||
| 98 : | /* | ||
| 99 : | Step 4 | ||
| 100 : | */ | ||
| 101 : | if(r+alpha*log(alpha/(b+w)) < t) goto S40; | ||
| 102 : | S70: | ||
| 103 : | /* | ||
| 104 : | Step 5 | ||
| 105 : | */ | ||
| 106 : | if(!(aa == a)) goto S80; | ||
| 107 : | genbet = w/(b+w); | ||
| 108 : | goto S90; | ||
| 109 : | S80: | ||
| 110 : | genbet = b/(b+w); | ||
| 111 : | S90: | ||
| 112 : | goto S230; | ||
| 113 : | S100: | ||
| 114 : | /* | ||
| 115 : | Algorithm BC | ||
| 116 : | Initialize | ||
| 117 : | */ | ||
| 118 : | if(qsame) goto S110; | ||
| 119 : | a = max(aa,bb); | ||
| 120 : | b = min(aa,bb); | ||
| 121 : | alpha = a+b; | ||
| 122 : | beta = 1.0/b; | ||
| 123 : | delta = 1.0+a-b; | ||
| 124 : | k1 = delta*(1.38889E-2+4.16667E-2*b)/(a*beta-0.777778); | ||
| 125 : | k2 = 0.25+(0.5+0.25/delta)*b; | ||
| 126 : | S110: | ||
| 127 : | S120: | ||
| 128 : | u1 = ranf(); | ||
| 129 : | /* | ||
| 130 : | Step 1 | ||
| 131 : | */ | ||
| 132 : | u2 = ranf(); | ||
| 133 : | if(u1 >= 0.5) goto S130; | ||
| 134 : | /* | ||
| 135 : | Step 2 | ||
| 136 : | */ | ||
| 137 : | y = u1*u2; | ||
| 138 : | z = u1*y; | ||
| 139 : | if(0.25*u2+z-y >= k1) goto S120; | ||
| 140 : | goto S170; | ||
| 141 : | S130: | ||
| 142 : | /* | ||
| 143 : | Step 3 | ||
| 144 : | */ | ||
| 145 : | z = pow(u1,2.0)*u2; | ||
| 146 : | if(!(z <= 0.25)) goto S160; | ||
| 147 : | v = beta*log(u1/(1.0-u1)); | ||
| 148 : | if(!(v > expmax)) goto S140; | ||
| 149 : | w = infnty; | ||
| 150 : | goto S150; | ||
| 151 : | S140: | ||
| 152 : | w = a*exp(v); | ||
| 153 : | S150: | ||
| 154 : | goto S200; | ||
| 155 : | S160: | ||
| 156 : | if(z >= k2) goto S120; | ||
| 157 : | S170: | ||
| 158 : | /* | ||
| 159 : | Step 4 | ||
| 160 : | Step 5 | ||
| 161 : | */ | ||
| 162 : | v = beta*log(u1/(1.0-u1)); | ||
| 163 : | if(!(v > expmax)) goto S180; | ||
| 164 : | w = infnty; | ||
| 165 : | goto S190; | ||
| 166 : | S180: | ||
| 167 : | w = a*exp(v); | ||
| 168 : | S190: | ||
| 169 : | if(alpha*(log(alpha/(b+w))+v)-1.3862944 < log(z)) goto S120; | ||
| 170 : | S200: | ||
| 171 : | /* | ||
| 172 : | Step 6 | ||
| 173 : | */ | ||
| 174 : | if(!(a == aa)) goto S210; | ||
| 175 : | genbet = w/(b+w); | ||
| 176 : | goto S220; | ||
| 177 : | S210: | ||
| 178 : | genbet = b/(b+w); | ||
| 179 : | S230: | ||
| 180 : | S220: | ||
| 181 : | return genbet; | ||
| 182 : | #undef expmax | ||
| 183 : | #undef infnty | ||
| 184 : | } | ||
| 185 : | float genchi(float df) | ||
| 186 : | /* | ||
| 187 : | ********************************************************************** | ||
| 188 : | float genchi(float df) | ||
| 189 : | Generate random value of CHIsquare variable | ||
| 190 : | Function | ||
| 191 : | Generates random deviate from the distribution of a chisquare | ||
| 192 : | with DF degrees of freedom random variable. | ||
| 193 : | Arguments | ||
| 194 : | df --> Degrees of freedom of the chisquare | ||
| 195 : | (Must be positive) | ||
| 196 : | |||
| 197 : | Method | ||
| 198 : | Uses relation between chisquare and gamma. | ||
| 199 : | ********************************************************************** | ||
| 200 : | */ | ||
| 201 : | { | ||
| 202 : | static float genchi; | ||
| 203 : | |||
| 204 : | if(!(df <= 0.0)) goto S10; | ||
| 205 : | fputs("DF <= 0 in GENCHI - ABORT\n",stderr); | ||
| 206 : | fprintf(stderr,"Value of DF: %16.6E\n",df); | ||
| 207 : | PyErr_SetString (PyExc_ValueError, "Described above."); | ||
| 208 : | return 0.0; | ||
| 209 : | S10: | ||
| 210 : | genchi = 2.0*gengam(1.0,df/2.0); | ||
| 211 : | return genchi; | ||
| 212 : | } | ||
| 213 : | float genexp(float av) | ||
| 214 : | /* | ||
| 215 : | ********************************************************************** | ||
| 216 : | float genexp(float av) | ||
| 217 : | GENerate EXPonential random deviate | ||
| 218 : | Function | ||
| 219 : | Generates a single random deviate from an exponential | ||
| 220 : | distribution with mean AV. | ||
| 221 : | Arguments | ||
| 222 : | av --> The mean of the exponential distribution from which | ||
| 223 : | a random deviate is to be generated. | ||
| 224 : | Method | ||
| 225 : | Renames SEXPO from TOMS as slightly modified by BWB to use RANF | ||
| 226 : | instead of SUNIF. | ||
| 227 : | For details see: | ||
| 228 : | Ahrens, J.H. and Dieter, U. | ||
| 229 : | Computer Methods for Sampling From the | ||
| 230 : | Exponential and Normal Distributions. | ||
| 231 : | Comm. ACM, 15,10 (Oct. 1972), 873 - 882. | ||
| 232 : | ********************************************************************** | ||
| 233 : | */ | ||
| 234 : | { | ||
| 235 : | static float genexp; | ||
| 236 : | |||
| 237 : | genexp = sexpo()*av; | ||
| 238 : | return genexp; | ||
| 239 : | } | ||
| 240 : | float genf(float dfn,float dfd) | ||
| 241 : | /* | ||
| 242 : | ********************************************************************** | ||
| 243 : | float genf(float dfn,float dfd) | ||
| 244 : | GENerate random deviate from the F distribution | ||
| 245 : | Function | ||
| 246 : | Generates a random deviate from the F (variance ratio) | ||
| 247 : | distribution with DFN degrees of freedom in the numerator | ||
| 248 : | and DFD degrees of freedom in the denominator. | ||
| 249 : | Arguments | ||
| 250 : | dfn --> Numerator degrees of freedom | ||
| 251 : | (Must be positive) | ||
| 252 : | dfd --> Denominator degrees of freedom | ||
| 253 : | (Must be positive) | ||
| 254 : | Method | ||
| 255 : | Directly generates ratio of chisquare variates | ||
| 256 : | ********************************************************************** | ||
| 257 : | */ | ||
| 258 : | { | ||
| 259 : | static float genf,xden,xnum; | ||
| 260 : | |||
| 261 : | if(!(dfn <= 0.0 || dfd <= 0.0)) goto S10; | ||
| 262 : | fputs("Degrees of freedom nonpositive in GENF - abort!\n",stderr); | ||
| 263 : | fprintf(stderr,"DFN value: %16.6EDFD value: %16.6E\n",dfn,dfd); | ||
| 264 : | PyErr_SetString (PyExc_ValueError, "Described above."); | ||
| 265 : | return 0.0; | ||
| 266 : | S10: | ||
| 267 : | xnum = genchi(dfn)/dfn; | ||
| 268 : | /* | ||
| 269 : | GENF = ( GENCHI( DFN ) / DFN ) / ( GENCHI( DFD ) / DFD ) | ||
| 270 : | */ | ||
| 271 : | xden = genchi(dfd)/dfd; | ||
| 272 : | if(!(xden <= 9.999999999998E-39*xnum)) goto S20; | ||
| 273 : | fputs(" GENF - generated numbers would cause overflow",stderr); | ||
| 274 : | fprintf(stderr," Numerator %16.6E Denominator %16.6E\n",xnum,xden); | ||
| 275 : | fputs(" GENF returning 1.0E38",stderr); | ||
| 276 : | genf = 1.0E38; | ||
| 277 : | goto S30; | ||
| 278 : | S20: | ||
| 279 : | genf = xnum/xden; | ||
| 280 : | S30: | ||
| 281 : | return genf; | ||
| 282 : | } | ||
| 283 : | float gengam(float a,float r) | ||
| 284 : | /* | ||
| 285 : | ********************************************************************** | ||
| 286 : | float gengam(float a,float r) | ||
| 287 : | GENerates random deviates from GAMma distribution | ||
| 288 : | Function | ||
| 289 : | Generates random deviates from the gamma distribution whose | ||
| 290 : | density is | ||
| 291 : | (A**R)/Gamma(R) * X**(R-1) * Exp(-A*X) | ||
| 292 : | Arguments | ||
| 293 : | a --> Location parameter of Gamma distribution | ||
| 294 : | r --> Shape parameter of Gamma distribution | ||
| 295 : | Method | ||
| 296 : | Renames SGAMMA from TOMS as slightly modified by BWB to use RANF | ||
| 297 : | instead of SUNIF. | ||
| 298 : | For details see: | ||
| 299 : | (Case R >= 1.0) | ||
| 300 : | Ahrens, J.H. and Dieter, U. | ||
| 301 : | Generating Gamma Variates by a | ||
| 302 : | Modified Rejection Technique. | ||
| 303 : | Comm. ACM, 25,1 (Jan. 1982), 47 - 54. | ||
| 304 : | Algorithm GD | ||
| 305 : | (Case 0.0 <= R <= 1.0) | ||
| 306 : | Ahrens, J.H. and Dieter, U. | ||
| 307 : | Computer Methods for Sampling from Gamma, | ||
| 308 : | Beta, Poisson and Binomial Distributions. | ||
| 309 : | Computing, 12 (1974), 223-246/ | ||
| 310 : | Adapted algorithm GS. | ||
| 311 : | ********************************************************************** | ||
| 312 : | */ | ||
| 313 : | { | ||
| 314 : | static float gengam; | ||
| 315 : | |||
| 316 : | gengam = sgamma(r); | ||
| 317 : | gengam /= a; | ||
| 318 : | return gengam; | ||
| 319 : | } | ||
| 320 : | void genmn(float *parm,float *x,float *work) | ||
| 321 : | /* | ||
| 322 : | ********************************************************************** | ||
| 323 : | void genmn(float *parm,float *x,float *work) | ||
| 324 : | GENerate Multivariate Normal random deviate | ||
| 325 : | Arguments | ||
| 326 : | parm --> Parameters needed to generate multivariate normal | ||
| 327 : | deviates (MEANV and Cholesky decomposition of | ||
| 328 : | COVM). Set by a previous call to SETGMN. | ||
| 329 : | 1 : 1 - size of deviate, P | ||
| 330 : | 2 : P + 1 - mean vector | ||
| 331 : | P+2 : P*(P+3)/2 + 1 - upper half of cholesky | ||
| 332 : | decomposition of cov matrix | ||
| 333 : | x <-- Vector deviate generated. | ||
| 334 : | work <--> Scratch array | ||
| 335 : | Method | ||
| 336 : | 1) Generate P independent standard normal deviates - Ei ~ N(0,1) | ||
| 337 : | 2) Using Cholesky decomposition find A s.t. trans(A)*A = COVM | ||
| 338 : | 3) trans(A)E + MEANV ~ N(MEANV,COVM) | ||
| 339 : | ********************************************************************** | ||
| 340 : | */ | ||
| 341 : | { | ||
| 342 : | static long i,icount,j,p,D1,D2,D3,D4; | ||
| 343 : | static float ae; | ||
| 344 : | |||
| 345 : | p = (long) (*parm); | ||
| 346 : | /* | ||
| 347 : | Generate P independent normal deviates - WORK ~ N(0,1) | ||
| 348 : | */ | ||
| 349 : | for(i=1; i<=p; i++) *(work+i-1) = snorm(); | ||
| 350 : | for(i=1,D3=1,D4=(p-i+D3)/D3; D4>0; D4--,i+=D3) { | ||
| 351 : | /* | ||
| 352 : | PARM (P+2 : P*(P+3)/2 + 1) contains A, the Cholesky | ||
| 353 : | decomposition of the desired covariance matrix. | ||
| 354 : | trans(A)(1,1) = PARM(P+2) | ||
| 355 : | trans(A)(2,1) = PARM(P+3) | ||
| 356 : | trans(A)(2,2) = PARM(P+2+P) | ||
| 357 : | trans(A)(3,1) = PARM(P+4) | ||
| 358 : | trans(A)(3,2) = PARM(P+3+P) | ||
| 359 : | trans(A)(3,3) = PARM(P+2-1+2P) ... | ||
| 360 : | trans(A)*WORK + MEANV ~ N(MEANV,COVM) | ||
| 361 : | */ | ||
| 362 : | icount = 0; | ||
| 363 : | ae = 0.0; | ||
| 364 : | for(j=1,D1=1,D2=(i-j+D1)/D1; D2>0; D2--,j+=D1) { | ||
| 365 : | icount += (j-1); | ||
| 366 : | ae += (*(parm+i+(j-1)*p-icount+p)**(work+j-1)); | ||
| 367 : | } | ||
| 368 : | *(x+i-1) = ae+*(parm+i); | ||
| 369 : | } | ||
| 370 : | } | ||
| 371 : | void genmul(long n,float *p,long ncat,long *ix) | ||
| 372 : | /* | ||
| 373 : | ********************************************************************** | ||
| 374 : | |||
| 375 : | void genmul(int n,float *p,int ncat,int *ix) | ||
| 376 : | GENerate an observation from the MULtinomial distribution | ||
| 377 : | Arguments | ||
| 378 : | N --> Number of events that will be classified into one of | ||
| 379 : | the categories 1..NCAT | ||
| 380 : | P --> Vector of probabilities. P(i) is the probability that | ||
| 381 : | an event will be classified into category i. Thus, P(i) | ||
| 382 : | must be [0,1]. Only the first NCAT-1 P(i) must be defined | ||
| 383 : | since P(NCAT) is 1.0 minus the sum of the first | ||
| 384 : | NCAT-1 P(i). | ||
| 385 : | NCAT --> Number of categories. Length of P and IX. | ||
| 386 : | IX <-- Observation from multinomial distribution. All IX(i) | ||
| 387 : | will be nonnegative and their sum will be N. | ||
| 388 : | Method | ||
| 389 : | Algorithm from page 559 of | ||
| 390 : | |||
| 391 : | Devroye, Luc | ||
| 392 : | |||
| 393 : | Non-Uniform Random Variate Generation. Springer-Verlag, | ||
| 394 : | New York, 1986. | ||
| 395 : | |||
| 396 : | ********************************************************************** | ||
| 397 : | */ | ||
| 398 : | { | ||
| 399 : | static float prob,ptot,sum; | ||
| 400 : | static long i,icat,ntot; | ||
| 401 : | if(n < 0) ftnstop("N < 0 in GENMUL"); | ||
| 402 : | if(ncat <= 1) ftnstop("NCAT <= 1 in GENMUL"); | ||
| 403 : | ptot = 0.0F; | ||
| 404 : | for(i=0; i<ncat-1; i++) { | ||
| 405 : | if(*(p+i) < 0.0F) ftnstop("Some P(i) < 0 in GENMUL"); | ||
| 406 : | if(*(p+i) > 1.0F) ftnstop("Some P(i) > 1 in GENMUL"); | ||
| 407 : | ptot += *(p+i); | ||
| 408 : | } | ||
| 409 : | if(ptot > 0.99999F) ftnstop("Sum of P(i) > 1 in GENMUL"); | ||
| 410 : | /* | ||
| 411 : | Initialize variables | ||
| 412 : | */ | ||
| 413 : | ntot = n; | ||
| 414 : | sum = 1.0F; | ||
| 415 : | for(i=0; i<ncat; i++) ix[i] = 0; | ||
| 416 : | /* | ||
| 417 : | Generate the observation | ||
| 418 : | */ | ||
| 419 : | for(icat=0; icat<ncat-1; icat++) { | ||
| 420 : | prob = *(p+icat)/sum; | ||
| 421 : | *(ix+icat) = ignbin(ntot,prob); | ||
| 422 : | ntot -= *(ix+icat); | ||
| 423 : | if(ntot <= 0) return; | ||
| 424 : | sum -= *(p+icat); | ||
| 425 : | } | ||
| 426 : | *(ix+ncat-1) = ntot; | ||
| 427 : | /* | ||
| 428 : | Finished | ||
| 429 : | */ | ||
| 430 : | return; | ||
| 431 : | } | ||
| 432 : | float gennch(float df,float xnonc) | ||
| 433 : | /* | ||
| 434 : | ********************************************************************** | ||
| 435 : | float gennch(float df,float xnonc) | ||
| 436 : | Generate random value of Noncentral CHIsquare variable | ||
| 437 : | Function | ||
| 438 : | Generates random deviate from the distribution of a noncentral | ||
| 439 : | chisquare with DF degrees of freedom and noncentrality parameter | ||
| 440 : | xnonc. | ||
| 441 : | Arguments | ||
| 442 : | df --> Degrees of freedom of the chisquare | ||
| 443 : | (Must be > 1.0) | ||
| 444 : | xnonc --> Noncentrality parameter of the chisquare | ||
| 445 : | (Must be >= 0.0) | ||
| 446 : | Method | ||
| 447 : | Uses fact that noncentral chisquare is the sum of a chisquare | ||
| 448 : | deviate with DF-1 degrees of freedom plus the square of a normal | ||
| 449 : | deviate with mean XNONC and standard deviation 1. | ||
| 450 : | ********************************************************************** | ||
| 451 : | */ | ||
| 452 : | { | ||
| 453 : | static float gennch; | ||
| 454 : | |||
| 455 : | if(!(df <= 1.0 || xnonc < 0.0)) goto S10; | ||
| 456 : | fputs("DF <= 1 or XNONC < 0 in GENNCH - ABORT\n",stderr); | ||
| 457 : | fprintf(stderr,"Value of DF: %16.6E Value of XNONC%16.6E\n",df,xnonc); | ||
| 458 : | PyErr_SetString (PyExc_ValueError, "Described above."); | ||
| 459 : | return 0.0; | ||
| 460 : | S10: | ||
| 461 : | gennch = genchi(df-1.0)+pow(gennor(sqrt(xnonc),1.0),2.0); | ||
| 462 : | return gennch; | ||
| 463 : | } | ||
| 464 : | float gennf(float dfn,float dfd,float xnonc) | ||
| 465 : | /* | ||
| 466 : | ********************************************************************** | ||
| 467 : | float gennf(float dfn,float dfd,float xnonc) | ||
| 468 : | GENerate random deviate from the Noncentral F distribution | ||
| 469 : | Function | ||
| 470 : | Generates a random deviate from the noncentral F (variance ratio) | ||
| 471 : | distribution with DFN degrees of freedom in the numerator, and DFD | ||
| 472 : | degrees of freedom in the denominator, and noncentrality parameter | ||
| 473 : | XNONC. | ||
| 474 : | Arguments | ||
| 475 : | dfn --> Numerator degrees of freedom | ||
| 476 : | (Must be >= 1.0) | ||
| 477 : | dfd --> Denominator degrees of freedom | ||
| 478 : | (Must be positive) | ||
| 479 : | xnonc --> Noncentrality parameter | ||
| 480 : | (Must be nonnegative) | ||
| 481 : | Method | ||
| 482 : | Directly generates ratio of noncentral numerator chisquare variate | ||
| 483 : | to central denominator chisquare variate. | ||
| 484 : | ********************************************************************** | ||
| 485 : | */ | ||
| 486 : | { | ||
| 487 : | static float gennf,xden,xnum; | ||
| 488 : | static long qcond; | ||
| 489 : | |||
| 490 : | qcond = dfn <= 1.0 || dfd <= 0.0 || xnonc < 0.0; | ||
| 491 : | if(!qcond) goto S10; | ||
| 492 : | fputs("In GENNF - Either (1) Numerator DF <= 1.0 or\n",stderr); | ||
| 493 : | fputs("(2) Denominator DF < 0.0 or \n",stderr); | ||
| 494 : | fputs("(3) Noncentrality parameter < 0.0\n",stderr); | ||
| 495 : | fprintf(stderr, | ||
| 496 : | "DFN value: %16.6EDFD value: %16.6EXNONC value: \n%16.6E\n",dfn,dfd, | ||
| 497 : | xnonc); | ||
| 498 : | PyErr_SetString (PyExc_ValueError, "Described above."); | ||
| 499 : | return 0.0; | ||
| 500 : | S10: | ||
| 501 : | xnum = gennch(dfn,xnonc)/dfn; | ||
| 502 : | /* | ||
| 503 : | GENNF = ( GENNCH( DFN, XNONC ) / DFN ) / ( GENCHI( DFD ) / DFD ) | ||
| 504 : | */ | ||
| 505 : | xden = genchi(dfd)/dfd; | ||
| 506 : | if(!(xden <= 9.999999999998E-39*xnum)) goto S20; | ||
| 507 : | fputs(" GENNF - generated numbers would cause overflow",stderr); | ||
| 508 : | fprintf(stderr," Numerator %16.6E Denominator %16.6E\n",xnum,xden); | ||
| 509 : | fputs(" GENNF returning 1.0E38",stderr); | ||
| 510 : | gennf = 1.0E38; | ||
| 511 : | goto S30; | ||
| 512 : | S20: | ||
| 513 : | gennf = xnum/xden; | ||
| 514 : | S30: | ||
| 515 : | return gennf; | ||
| 516 : | } | ||
| 517 : | float gennor(float av,float sd) | ||
| 518 : | /* | ||
| 519 : | ********************************************************************** | ||
| 520 : | float gennor(float av,float sd) | ||
| 521 : | GENerate random deviate from a NORmal distribution | ||
| 522 : | Function | ||
| 523 : | Generates a single random deviate from a normal distribution | ||
| 524 : | with mean, AV, and standard deviation, SD. | ||
| 525 : | Arguments | ||
| 526 : | av --> Mean of the normal distribution. | ||
| 527 : | sd --> Standard deviation of the normal distribution. | ||
| 528 : | Method | ||
| 529 : | Renames SNORM from TOMS as slightly modified by BWB to use RANF | ||
| 530 : | instead of SUNIF. | ||
| 531 : | For details see: | ||
| 532 : | Ahrens, J.H. and Dieter, U. | ||
| 533 : | Extensions of Forsythe's Method for Random | ||
| 534 : | Sampling from the Normal Distribution. | ||
| 535 : | Math. Comput., 27,124 (Oct. 1973), 927 - 937. | ||
| 536 : | ********************************************************************** | ||
| 537 : | */ | ||
| 538 : | { | ||
| 539 : | static float gennor; | ||
| 540 : | |||
| 541 : | gennor = sd*snorm()+av; | ||
| 542 : | return gennor; | ||
| 543 : | } | ||
| 544 : | void genprm(long *iarray,int larray) | ||
| 545 : | /* | ||
| 546 : | ********************************************************************** | ||
| 547 : | void genprm(long *iarray,int larray) | ||
| 548 : | GENerate random PeRMutation of iarray | ||
| 549 : | Arguments | ||
| 550 : | iarray <--> On output IARRAY is a random permutation of its | ||
| 551 : | value on input | ||
| 552 : | larray <--> Length of IARRAY | ||
| 553 : | ********************************************************************** | ||
| 554 : | */ | ||
| 555 : | { | ||
| 556 : | static long i,itmp,iwhich,D1,D2; | ||
| 557 : | |||
| 558 : | for(i=1,D1=1,D2=(larray-i+D1)/D1; D2>0; D2--,i+=D1) { | ||
| 559 : | iwhich = ignuin(i,larray); | ||
| 560 : | itmp = *(iarray+iwhich-1); | ||
| 561 : | *(iarray+iwhich-1) = *(iarray+i-1); | ||
| 562 : | *(iarray+i-1) = itmp; | ||
| 563 : | } | ||
| 564 : | } | ||
| 565 : | float genunf(float low,float high) | ||
| 566 : | /* | ||
| 567 : | ********************************************************************** | ||
| 568 : | float genunf(float low,float high) | ||
| 569 : | GeNerate Uniform Real between LOW and HIGH | ||
| 570 : | Function | ||
| 571 : | Generates a real uniformly distributed between LOW and HIGH. | ||
| 572 : | Arguments | ||
| 573 : | low --> Low bound (exclusive) on real value to be generated | ||
| 574 : | high --> High bound (exclusive) on real value to be generated | ||
| 575 : | ********************************************************************** | ||
| 576 : | */ | ||
| 577 : | { | ||
| 578 : | static float genunf; | ||
| 579 : | |||
| 580 : | if(!(low > high)) goto S10; | ||
| 581 : | fprintf(stderr,"LOW > HIGH in GENUNF: LOW %16.6E HIGH: %16.6E\n",low,high); | ||
| 582 : | fputs("Abort\n",stderr); | ||
| 583 : | PyErr_SetString (PyExc_ValueError, "Described above."); | ||
| 584 : | return 0.0; | ||
| 585 : | S10: | ||
| 586 : | genunf = low+(high-low)*ranf(); | ||
| 587 : | return genunf; | ||
| 588 : | } | ||
| 589 : | void gscgn(long getset,long *g) | ||
| 590 : | /* | ||
| 591 : | ********************************************************************** | ||
| 592 : | void gscgn(long getset,long *g) | ||
| 593 : | Get/Set GeNerator | ||
| 594 : | Gets or returns in G the number of the current generator | ||
| 595 : | Arguments | ||
| 596 : | getset --> 0 Get | ||
| 597 : | 1 Set | ||
| 598 : | g <-- Number of the current random number generator (1..32) | ||
| 599 : | ********************************************************************** | ||
| 600 : | */ | ||
| 601 : | { | ||
| 602 : | #define numg 32L | ||
| 603 : | static long curntg = 1; | ||
| 604 : | if(getset == 0) *g = curntg; | ||
| 605 : | else { | ||
| 606 : | if(*g < 0 || *g > numg) { | ||
| 607 : | fputs(" Generator number out of range in GSCGN\n",stderr); | ||
| 608 : | PyErr_SetString (PyExc_ValueError, "Described above."); | ||
| 609 : | return; | ||
| 610 : | } | ||
| 611 : | curntg = *g; | ||
| 612 : | } | ||
| 613 : | #undef numg | ||
| 614 : | } | ||
| 615 : | void gsrgs(long getset,long *qvalue) | ||
| 616 : | /* | ||
| 617 : | ********************************************************************** | ||
| 618 : | void gsrgs(long getset,long *qvalue) | ||
| 619 : | Get/Set Random Generators Set | ||
| 620 : | Gets or sets whether random generators set (initialized). | ||
| 621 : | Initially (data statement) state is not set | ||
| 622 : | If getset is 1 state is set to qvalue | ||
| 623 : | If getset is 0 state returned in qvalue | ||
| 624 : | ********************************************************************** | ||
| 625 : | */ | ||
| 626 : | { | ||
| 627 : | static long qinit = 0; | ||
| 628 : | |||
| 629 : | if(getset == 0) *qvalue = qinit; | ||
| 630 : | else qinit = *qvalue; | ||
| 631 : | } | ||
| 632 : | void gssst(long getset,long *qset) | ||
| 633 : | /* | ||
| 634 : | ********************************************************************** | ||
| 635 : | void gssst(long getset,long *qset) | ||
| 636 : | Get or Set whether Seed is Set | ||
| 637 : | Initialize to Seed not Set | ||
| 638 : | If getset is 1 sets state to Seed Set | ||
| 639 : | If getset is 0 returns T in qset if Seed Set | ||
| 640 : | Else returns F in qset | ||
| 641 : | ********************************************************************** | ||
| 642 : | */ | ||
| 643 : | { | ||
| 644 : | static long qstate = 0; | ||
| 645 : | if(getset != 0) qstate = 1; | ||
| 646 : | else *qset = qstate; | ||
| 647 : | } | ||
| 648 : | long ignbin(long n,float pp) | ||
| 649 : | /* | ||
| 650 : | ********************************************************************** | ||
| 651 : | long ignbin(long n,float pp) | ||
| 652 : | GENerate BINomial random deviate | ||
| 653 : | Function | ||
| 654 : | Generates a single random deviate from a binomial | ||
| 655 : | distribution whose number of trials is N and whose | ||
| 656 : | probability of an event in each trial is P. | ||
| 657 : | Arguments | ||
| 658 : | n --> The number of trials in the binomial distribution | ||
| 659 : | from which a random deviate is to be generated. | ||
| 660 : | p --> The probability of an event in each trial of the | ||
| 661 : | binomial distribution from which a random deviate | ||
| 662 : | is to be generated. | ||
| 663 : | ignbin <-- A random deviate yielding the number of events | ||
| 664 : | from N independent trials, each of which has | ||
| 665 : | a probability of event P. | ||
| 666 : | Method | ||
| 667 : | This is algorithm BTPE from: | ||
| 668 : | Kachitvichyanukul, V. and Schmeiser, B. W. | ||
| 669 : | Binomial Random Variate Generation. | ||
| 670 : | Communications of the ACM, 31, 2 | ||
| 671 : | (February, 1988) 216. | ||
| 672 : | ********************************************************************** | ||
| 673 : | SUBROUTINE BTPEC(N,PP,ISEED,JX) | ||
| 674 : | BINOMIAL RANDOM VARIATE GENERATOR | ||
| 675 : | MEAN .LT. 30 -- INVERSE CDF | ||
| 676 : | MEAN .GE. 30 -- ALGORITHM BTPE: ACCEPTANCE-REJECTION VIA | ||
| 677 : | FOUR REGION COMPOSITION. THE FOUR REGIONS ARE A TRIANGLE | ||
| 678 : | (SYMMETRIC IN THE CENTER), A PAIR OF PARALLELOGRAMS (ABOVE | ||
| 679 : | THE TRIANGLE), AND EXPONENTIAL LEFT AND RIGHT TAILS. | ||
| 680 : | BTPE REFERS TO BINOMIAL-TRIANGLE-PARALLELOGRAM-EXPONENTIAL. | ||
| 681 : | BTPEC REFERS TO BTPE AND "COMBINED." THUS BTPE IS THE | ||
| 682 : | RESEARCH AND BTPEC IS THE IMPLEMENTATION OF A COMPLETE | ||
| 683 : | USABLE ALGORITHM. | ||
| 684 : | REFERENCE: VORATAS KACHITVICHYANUKUL AND BRUCE SCHMEISER, | ||
| 685 : | "BINOMIAL RANDOM VARIATE GENERATION," | ||
| 686 : | COMMUNICATIONS OF THE ACM, FORTHCOMING | ||
| 687 : | WRITTEN: SEPTEMBER 1980. | ||
| 688 : | LAST REVISED: MAY 1985, JULY 1987 | ||
| 689 : | REQUIRED SUBPROGRAM: RAND() -- A UNIFORM (0,1) RANDOM NUMBER | ||
| 690 : | GENERATOR | ||
| 691 : | ARGUMENTS | ||
| 692 : | N : NUMBER OF BERNOULLI TRIALS (INPUT) | ||
| 693 : | PP : PROBABILITY OF SUCCESS IN EACH TRIAL (INPUT) | ||
| 694 : | ISEED: RANDOM NUMBER SEED (INPUT AND OUTPUT) | ||
| 695 : | JX: RANDOMLY GENERATED OBSERVATION (OUTPUT) | ||
| 696 : | VARIABLES | ||
| 697 : | PSAVE: VALUE OF PP FROM THE LAST CALL TO BTPEC | ||
| 698 : | NSAVE: VALUE OF N FROM THE LAST CALL TO BTPEC | ||
| 699 : | XNP: VALUE OF THE MEAN FROM THE LAST CALL TO BTPEC | ||
| 700 : | P: PROBABILITY USED IN THE GENERATION PHASE OF BTPEC | ||
| 701 : | FFM: TEMPORARY VARIABLE EQUAL TO XNP + P | ||
| 702 : | M: INTEGER VALUE OF THE CURRENT MODE | ||
| 703 : | FM: FLOATING POINT VALUE OF THE CURRENT MODE | ||
| 704 : | XNPQ: TEMPORARY VARIABLE USED IN SETUP AND SQUEEZING STEPS | ||
| 705 : | P1: AREA OF THE TRIANGLE | ||
| 706 : | C: HEIGHT OF THE PARALLELOGRAMS | ||
| 707 : | XM: CENTER OF THE TRIANGLE | ||
| 708 : | XL: LEFT END OF THE TRIANGLE | ||
| 709 : | XR: RIGHT END OF THE TRIANGLE | ||
| 710 : | AL: TEMPORARY VARIABLE | ||
| 711 : | XLL: RATE FOR THE LEFT EXPONENTIAL TAIL | ||
| 712 : | XLR: RATE FOR THE RIGHT EXPONENTIAL TAIL | ||
| 713 : | P2: AREA OF THE PARALLELOGRAMS | ||
| 714 : | P3: AREA OF THE LEFT EXPONENTIAL TAIL | ||
| 715 : | P4: AREA OF THE RIGHT EXPONENTIAL TAIL | ||
| 716 : | U: A U(0,P4) RANDOM VARIATE USED FIRST TO SELECT ONE OF THE | ||
| 717 : | FOUR REGIONS AND THEN CONDITIONALLY TO GENERATE A VALUE | ||
| 718 : | FROM THE REGION | ||
| 719 : | V: A U(0,1) RANDOM NUMBER USED TO GENERATE THE RANDOM VALUE | ||
| 720 : | (REGION 1) OR TRANSFORMED INTO THE VARIATE TO ACCEPT OR | ||
| 721 : | REJECT THE CANDIDATE VALUE | ||
| 722 : | IX: INTEGER CANDIDATE VALUE | ||
| 723 : | X: PRELIMINARY CONTINUOUS CANDIDATE VALUE IN REGION 2 LOGIC | ||
| 724 : | AND A FLOATING POINT IX IN THE ACCEPT/REJECT LOGIC | ||
| 725 : | K: ABSOLUTE VALUE OF (IX-M) | ||
| 726 : | F: THE HEIGHT OF THE SCALED DENSITY FUNCTION USED IN THE | ||
| 727 : | ACCEPT/REJECT DECISION WHEN BOTH M AND IX ARE SMALL | ||
| 728 : | ALSO USED IN THE INVERSE TRANSFORMATION | ||
| 729 : | R: THE RATIO P/Q | ||
| 730 : | G: CONSTANT USED IN CALCULATION OF PROBABILITY | ||
| 731 : | MP: MODE PLUS ONE, THE LOWER INDEX FOR EXPLICIT CALCULATION | ||
| 732 : | OF F WHEN IX IS GREATER THAN M | ||
| 733 : | IX1: CANDIDATE VALUE PLUS ONE, THE LOWER INDEX FOR EXPLICIT | ||
| 734 : | CALCULATION OF F WHEN IX IS LESS THAN M | ||
| 735 : | I: INDEX FOR EXPLICIT CALCULATION OF F FOR BTPE | ||
| 736 : | AMAXP: MAXIMUM ERROR OF THE LOGARITHM OF NORMAL BOUND | ||
| 737 : | YNORM: LOGARITHM OF NORMAL BOUND | ||
| 738 : | ALV: NATURAL LOGARITHM OF THE ACCEPT/REJECT VARIATE V | ||
| 739 : | X1,F1,Z,W,Z2,X2,F2, AND W2 ARE TEMPORARY VARIABLES TO BE | ||
| 740 : | USED IN THE FINAL ACCEPT/REJECT TEST | ||
| 741 : | QN: PROBABILITY OF NO SUCCESS IN N TRIALS | ||
| 742 : | REMARK | ||
| 743 : | IX AND JX COULD LOGICALLY BE THE SAME VARIABLE, WHICH WOULD | ||
| 744 : | SAVE A MEMORY POSITION AND A LINE OF CODE. HOWEVER, SOME | ||
| 745 : | COMPILERS (E.G.,CDC MNF) OPTIMIZE BETTER WHEN THE ARGUMENTS | ||
| 746 : | ARE NOT INVOLVED. | ||
| 747 : | ISEED NEEDS TO BE DOUBLE PRECISION IF THE IMSL ROUTINE | ||
| 748 : | GGUBFS IS USED TO GENERATE UNIFORM RANDOM NUMBER, OTHERWISE | ||
| 749 : | TYPE OF ISEED SHOULD BE DICTATED BY THE UNIFORM GENERATOR | ||
| 750 : | ********************************************************************** | ||
| 751 : | *****DETERMINE APPROPRIATE ALGORITHM AND WHETHER SETUP IS NECESSARY | ||
| 752 : | */ | ||
| 753 : | { | ||
| 754 : | static float psave = -1.0; | ||
| 755 : | static long nsave = -1; | ||
| 756 : | static long ignbin,i,ix,ix1,k,m,mp,T1; | ||
| 757 : | static float al,alv,amaxp,c,f,f1,f2,ffm,fm,g,p,p1,p2,p3,p4,q,qn,r,u,v,w,w2,x,x1, | ||
| 758 : | x2,xl,xll,xlr,xm,xnp,xnpq,xr,ynorm,z,z2; | ||
| 759 : | |||
| 760 : | if(pp != psave) goto S10; | ||
| 761 : | if(n != nsave) goto S20; | ||
| 762 : | if(xnp < 30.0) goto S150; | ||
| 763 : | goto S30; | ||
| 764 : | S10: | ||
| 765 : | /* | ||
| 766 : | *****SETUP, PERFORM ONLY WHEN PARAMETERS CHANGE | ||
| 767 : | */ | ||
| 768 : | psave = pp; | ||
| 769 : | p = min(psave,1.0-psave); | ||
| 770 : | q = 1.0-p; | ||
| 771 : | S20: | ||
| 772 : | xnp = n*p; | ||
| 773 : | nsave = n; | ||
| 774 : | if(xnp < 30.0) goto S140; | ||
| 775 : | ffm = xnp+p; | ||
| 776 : | m = ffm; | ||
| 777 : | fm = m; | ||
| 778 : | xnpq = xnp*q; | ||
| 779 : | p1 = (long) (2.195*sqrt(xnpq)-4.6*q)+0.5; | ||
| 780 : | xm = fm+0.5; | ||
| 781 : | xl = xm-p1; | ||
| 782 : | xr = xm+p1; | ||
| 783 : | c = 0.134+20.5/(15.3+fm); | ||
| 784 : | al = (ffm-xl)/(ffm-xl*p); | ||
| 785 : | xll = al*(1.0+0.5*al); | ||
| 786 : | al = (xr-ffm)/(xr*q); | ||
| 787 : | xlr = al*(1.0+0.5*al); | ||
| 788 : | p2 = p1*(1.0+c+c); | ||
| 789 : | p3 = p2+c/xll; | ||
| 790 : | p4 = p3+c/xlr; | ||
| 791 : | S30: | ||
| 792 : | /* | ||
| 793 : | *****GENERATE VARIATE | ||
| 794 : | */ | ||
| 795 : | u = ranf()*p4; | ||
| 796 : | v = ranf(); | ||
| 797 : | /* | ||
| 798 : | TRIANGULAR REGION | ||
| 799 : | */ | ||
| 800 : | if(u > p1) goto S40; | ||
| 801 : | ix = xm-p1*v+u; | ||
| 802 : | goto S170; | ||
| 803 : | S40: | ||
| 804 : | /* | ||
| 805 : | PARALLELOGRAM REGION | ||
| 806 : | */ | ||
| 807 : | if(u > p2) goto S50; | ||
| 808 : | x = xl+(u-p1)/c; | ||
| 809 : | v = v*c+1.0-ABS(xm-x)/p1; | ||
| 810 : | if(v > 1.0 || v <= 0.0) goto S30; | ||
| 811 : | ix = x; | ||
| 812 : | goto S70; | ||
| 813 : | S50: | ||
| 814 : | /* | ||
| 815 : | LEFT TAIL | ||
| 816 : | */ | ||
| 817 : | if(u > p3) goto S60; | ||
| 818 : | ix = xl+log(v)/xll; | ||
| 819 : | if(ix < 0) goto S30; | ||
| 820 : | v *= ((u-p2)*xll); | ||
| 821 : | goto S70; | ||
| 822 : | S60: | ||
| 823 : | /* | ||
| 824 : | RIGHT TAIL | ||
| 825 : | */ | ||
| 826 : | ix = xr-log(v)/xlr; | ||
| 827 : | if(ix > n) goto S30; | ||
| 828 : | v *= ((u-p3)*xlr); | ||
| 829 : | S70: | ||
| 830 : | /* | ||
| 831 : | *****DETERMINE APPROPRIATE WAY TO PERFORM ACCEPT/REJECT TEST | ||
| 832 : | */ | ||
| 833 : | k = ABS(ix-m); | ||
| 834 : | if(k > 20 && k < xnpq/2-1) goto S130; | ||
| 835 : | /* | ||
| 836 : | EXPLICIT EVALUATION | ||
| 837 : | */ | ||
| 838 : | f = 1.0; | ||
| 839 : | r = p/q; | ||
| 840 : | g = (n+1)*r; | ||
| 841 : | T1 = m-ix; | ||
| 842 : | if(T1 < 0) goto S80; | ||
| 843 : | else if(T1 == 0) goto S120; | ||
| 844 : | else goto S100; | ||
| 845 : | S80: | ||
| 846 : | mp = m+1; | ||
| 847 : | for(i=mp; i<=ix; i++) f *= (g/i-r); | ||
| 848 : | goto S120; | ||
| 849 : | S100: | ||
| 850 : | ix1 = ix+1; | ||
| 851 : | for(i=ix1; i<=m; i++) f /= (g/i-r); | ||
| 852 : | S120: | ||
| 853 : | if(v <= f) goto S170; | ||
| 854 : | goto S30; | ||
| 855 : | S130: | ||
| 856 : | /* | ||
| 857 : | SQUEEZING USING UPPER AND LOWER BOUNDS ON ALOG(F(X)) | ||
| 858 : | */ | ||
| 859 : | amaxp = k/xnpq*((k*(k/3.0+0.625)+0.1666666666666)/xnpq+0.5); | ||
| 860 : | ynorm = -(k*k/(2.0*xnpq)); | ||
| 861 : | alv = log(v); | ||
| 862 : | if(alv < ynorm-amaxp) goto S170; | ||
| 863 : | if(alv > ynorm+amaxp) goto S30; | ||
| 864 : | /* | ||
| 865 : | STIRLING'S FORMULA TO MACHINE ACCURACY FOR | ||
| 866 : | THE FINAL ACCEPTANCE/REJECTION TEST | ||
| 867 : | */ | ||
| 868 : | x1 = ix+1.0; | ||
| 869 : | f1 = fm+1.0; | ||
| 870 : | z = n+1.0-fm; | ||
| 871 : | w = n-ix+1.0; | ||
| 872 : | z2 = z*z; | ||
| 873 : | x2 = x1*x1; | ||
| 874 : | f2 = f1*f1; | ||
| 875 : | w2 = w*w; | ||
| 876 : | if(alv <= xm*log(f1/x1)+(n-m+0.5)*log(z/w)+(ix-m)*log(w*p/(x1*q))+(13860.0- | ||
| 877 : | (462.0-(132.0-(99.0-140.0/f2)/f2)/f2)/f2)/f1/166320.0+(13860.0-(462.0- | ||
| 878 : | (132.0-(99.0-140.0/z2)/z2)/z2)/z2)/z/166320.0+(13860.0-(462.0-(132.0- | ||
| 879 : | (99.0-140.0/x2)/x2)/x2)/x2)/x1/166320.0+(13860.0-(462.0-(132.0-(99.0 | ||
| 880 : | -140.0/w2)/w2)/w2)/w2)/w/166320.0) goto S170; | ||
| 881 : | goto S30; | ||
| 882 : | S140: | ||
| 883 : | /* | ||
| 884 : | INVERSE CDF LOGIC FOR MEAN LESS THAN 30 | ||
| 885 : | */ | ||
| 886 : | qn = pow(q,(double)n); | ||
| 887 : | r = p/q; | ||
| 888 : | g = r*(n+1); | ||
| 889 : | S150: | ||
| 890 : | ix = 0; | ||
| 891 : | f = qn; | ||
| 892 : | u = ranf(); | ||
| 893 : | S160: | ||
| 894 : | if(u < f) goto S170; | ||
| 895 : | if(ix > 110) goto S150; | ||
| 896 : | u -= f; | ||
| 897 : | ix += 1; | ||
| 898 : | f *= (g/ix-r); | ||
| 899 : | goto S160; | ||
| 900 : | S170: | ||
| 901 : | if(psave > 0.5) ix = n-ix; | ||
| 902 : | ignbin = ix; | ||
| 903 : | return ignbin; | ||
| 904 : | } | ||
| 905 : | long ignnbn(long n,float p) | ||
| 906 : | /* | ||
| 907 : | ********************************************************************** | ||
| 908 : | |||
| 909 : | long ignnbn(long n,float p) | ||
| 910 : | GENerate Negative BiNomial random deviate | ||
| 911 : | Function | ||
| 912 : | Generates a single random deviate from a negative binomial | ||
| 913 : | distribution. | ||
| 914 : | Arguments | ||
| 915 : | N --> The number of trials in the negative binomial distribution | ||
| 916 : | from which a random deviate is to be generated. | ||
| 917 : | P --> The probability of an event. | ||
| 918 : | Method | ||
| 919 : | Algorithm from page 480 of | ||
| 920 : | |||
| 921 : | Devroye, Luc | ||
| 922 : | |||
| 923 : | Non-Uniform Random Variate Generation. Springer-Verlag, | ||
| 924 : | New York, 1986. | ||
| 925 : | ********************************************************************** | ||
| 926 : | */ | ||
| 927 : | { | ||
| 928 : | static long ignnbn; | ||
| 929 : | static float y,a,r; | ||
| 930 : | /* | ||
| 931 : | .. | ||
| 932 : | .. Executable Statements .. | ||
| 933 : | */ | ||
| 934 : | /* | ||
| 935 : | Check Arguments | ||
| 936 : | */ | ||
| 937 : | if(n < 0) ftnstop("N < 0 in IGNNBN"); | ||
| 938 : | if(p <= 0.0F) ftnstop("P <= 0 in IGNNBN"); | ||
| 939 : | if(p >= 1.0F) ftnstop("P >= 1 in IGNNBN"); | ||
| 940 : | /* | ||
| 941 : | Generate Y, a random gamma (n,(1-p)/p) variable | ||
| 942 : | */ | ||
| 943 : | r = (float)n; | ||
| 944 : | a = p/(1.0F-p); | ||
| 945 : | y = gengam(a,r); | ||
| 946 : | /* | ||
| 947 : | Generate a random Poisson(y) variable | ||
| 948 : | */ | ||
| 949 : | ignnbn = ignpoi(y); | ||
| 950 : | return ignnbn; | ||
| 951 : | } | ||
| 952 : | long ignpoi(float mu) | ||
| 953 : | /* | ||
| 954 : | ********************************************************************** | ||
| 955 : | long ignpoi(float mu) | ||
| 956 : | GENerate POIsson random deviate | ||
| 957 : | Function | ||
| 958 : | Generates a single random deviate from a Poisson | ||
| 959 : | distribution with mean AV. | ||
| 960 : | Arguments | ||
| 961 : | av --> The mean of the Poisson distribution from which | ||
| 962 : | a random deviate is to be generated. | ||
| 963 : | genexp <-- The random deviate. | ||
| 964 : | Method | ||
| 965 : | Renames KPOIS from TOMS as slightly modified by BWB to use RANF | ||
| 966 : | instead of SUNIF. | ||
| 967 : | For details see: | ||
| 968 : | Ahrens, J.H. and Dieter, U. | ||
| 969 : | Computer Generation of Poisson Deviates | ||
| 970 : | From Modified Normal Distributions. | ||
| 971 : | ACM Trans. Math. Software, 8, 2 | ||
| 972 : | (June 1982),163-179 | ||
| 973 : | ********************************************************************** | ||
| 974 : | ********************************************************************** | ||
| 975 : | |||
| 976 : | |||
| 977 : | P O I S S O N DISTRIBUTION | ||
| 978 : | |||
| 979 : | |||
| 980 : | ********************************************************************** | ||
| 981 : | ********************************************************************** | ||
| 982 : | |||
| 983 : | FOR DETAILS SEE: | ||
| 984 : | |||
| 985 : | AHRENS, J.H. AND DIETER, U. | ||
| 986 : | COMPUTER GENERATION OF POISSON DEVIATES | ||
| 987 : | FROM MODIFIED NORMAL DISTRIBUTIONS. | ||
| 988 : | ACM TRANS. MATH. SOFTWARE, 8,2 (JUNE 1982), 163 - 179. | ||
| 989 : | |||
| 990 : | (SLIGHTLY MODIFIED VERSION OF THE PROGRAM IN THE ABOVE ARTICLE) | ||
| 991 : | |||
| 992 : | ********************************************************************** | ||
| 993 : | INTEGER FUNCTION IGNPOI(IR,MU) | ||
| 994 : | INPUT: IR=CURRENT STATE OF BASIC RANDOM NUMBER GENERATOR | ||
| 995 : | MU=MEAN MU OF THE POISSON DISTRIBUTION | ||
| 996 : | OUTPUT: IGNPOI=SAMPLE FROM THE POISSON-(MU)-DISTRIBUTION | ||
| 997 : | MUPREV=PREVIOUS MU, MUOLD=MU AT LAST EXECUTION OF STEP P OR B. | ||
| 998 : | TABLES: COEFFICIENTS A0-A7 FOR STEP F. FACTORIALS FACT | ||
| 999 : | COEFFICIENTS A(K) - FOR PX = FK*V*V*SUM(A(K)*V**K)-DEL | ||
| 1000 : | SEPARATION OF CASES A AND B | ||
| 1001 : | */ | ||
| 1002 : | { | ||
| 1003 : | extern float fsign( float num, float sign ); | ||
| 1004 : | static float a0 = -0.5; | ||
| 1005 : | static float a1 = 0.3333333; | ||
| 1006 : | static float a2 = -0.2500068; | ||
| 1007 : | static float a3 = 0.2000118; | ||
| 1008 : | static float a4 = -0.1661269; | ||
| 1009 : | static float a5 = 0.1421878; | ||
| 1010 : | static float a6 = -0.1384794; | ||
| 1011 : | static float a7 = 0.125006; | ||
| 1012 : | static float muold = 0.0; | ||
| 1013 : | static float muprev = 0.0; | ||
| 1014 : | static float fact[10] = { | ||
| 1015 : | 1.0,1.0,2.0,6.0,24.0,120.0,720.0,5040.0,40320.0,362880.0 | ||
| 1016 : | }; | ||
| 1017 : | static long ignpoi,j,k,kflag,l,m; | ||
| 1018 : | static float b1,b2,c,c0,c1,c2,c3,d,del,difmuk,e,fk,fx,fy,g,omega,p,p0,px,py,q,s, | ||
| 1019 : | t,u,v,x,xx,pp[35]; | ||
| 1020 : | |||
| 1021 : | if(mu == muprev) goto S10; | ||
| 1022 : | if(mu < 10.0) goto S120; | ||
| 1023 : | /* | ||
| 1024 : | C A S E A. (RECALCULATION OF S,D,L IF MU HAS CHANGED) | ||
| 1025 : | */ | ||
| 1026 : | muprev = mu; | ||
| 1027 : | s = sqrt(mu); | ||
| 1028 : | d = 6.0*mu*mu; | ||
| 1029 : | /* | ||
| 1030 : | THE POISSON PROBABILITIES PK EXCEED THE DISCRETE NORMAL | ||
| 1031 : | PROBABILITIES FK WHENEVER K >= M(MU). L=IFIX(MU-1.1484) | ||
| 1032 : | IS AN UPPER BOUND TO M(MU) FOR ALL MU >= 10 . | ||
| 1033 : | */ | ||
| 1034 : | l = (long) (mu-1.1484); | ||
| 1035 : | S10: | ||
| 1036 : | /* | ||
| 1037 : | STEP N. NORMAL SAMPLE - SNORM(IR) FOR STANDARD NORMAL DEVIATE | ||
| 1038 : | */ | ||
| 1039 : | g = mu+s*snorm(); | ||
| 1040 : | if(g < 0.0) goto S20; | ||
| 1041 : | ignpoi = (long) (g); | ||
| 1042 : | /* | ||
| 1043 : | STEP I. IMMEDIATE ACCEPTANCE IF IGNPOI IS LARGE ENOUGH | ||
| 1044 : | */ | ||
| 1045 : | if(ignpoi >= l) return ignpoi; | ||
| 1046 : | /* | ||
| 1047 : | STEP S. SQUEEZE ACCEPTANCE - SUNIF(IR) FOR (0,1)-SAMPLE U | ||
| 1048 : | */ | ||
| 1049 : | fk = (float)ignpoi; | ||
| 1050 : | difmuk = mu-fk; | ||
| 1051 : | u = ranf(); | ||
| 1052 : | if(d*u >= difmuk*difmuk*difmuk) return ignpoi; | ||
| 1053 : | S20: | ||
| 1054 : | /* | ||
| 1055 : | STEP P. PREPARATIONS FOR STEPS Q AND H. | ||
| 1056 : | (RECALCULATIONS OF PARAMETERS IF NECESSARY) | ||
| 1057 : | .3989423=(2*PI)**(-.5) .416667E-1=1./24. .1428571=1./7. | ||
| 1058 : | THE QUANTITIES B1, B2, C3, C2, C1, C0 ARE FOR THE HERMITE | ||
| 1059 : | APPROXIMATIONS TO THE DISCRETE NORMAL PROBABILITIES FK. | ||
| 1060 : | C=.1069/MU GUARANTEES MAJORIZATION BY THE 'HAT'-FUNCTION. | ||
| 1061 : | */ | ||
| 1062 : | if(mu == muold) goto S30; | ||
| 1063 : | muold = mu; | ||
| 1064 : | omega = 0.3989423/s; | ||
| 1065 : | b1 = 4.166667E-2/mu; | ||
| 1066 : | b2 = 0.3*b1*b1; | ||
| 1067 : | c3 = 0.1428571*b1*b2; | ||
| 1068 : | c2 = b2-15.0*c3; | ||
| 1069 : | c1 = b1-6.0*b2+45.0*c3; | ||
| 1070 : | c0 = 1.0-b1+3.0*b2-15.0*c3; | ||
| 1071 : | c = 0.1069/mu; | ||
| 1072 : | S30: | ||
| 1073 : | if(g < 0.0) goto S50; | ||
| 1074 : | /* | ||
| 1075 : | 'SUBROUTINE' F IS CALLED (KFLAG=0 FOR CORRECT RETURN) | ||
| 1076 : | */ | ||
| 1077 : | kflag = 0; | ||
| 1078 : | goto S70; | ||
| 1079 : | S40: | ||
| 1080 : | /* | ||
| 1081 : | STEP Q. QUOTIENT ACCEPTANCE (RARE CASE) | ||
| 1082 : | */ | ||
| 1083 : | if(fy-u*fy <= py*exp(px-fx)) return ignpoi; | ||
| 1084 : | S50: | ||
| 1085 : | /* | ||
| 1086 : | STEP E. EXPONENTIAL SAMPLE - SEXPO(IR) FOR STANDARD EXPONENTIAL | ||
| 1087 : | DEVIATE E AND SAMPLE T FROM THE LAPLACE 'HAT' | ||
| 1088 : | (IF T <= -.6744 THEN PK < FK FOR ALL MU >= 10.) | ||
| 1089 : | */ | ||
| 1090 : | e = sexpo(); | ||
| 1091 : | u = ranf(); | ||
| 1092 : | u += (u-1.0); | ||
| 1093 : | t = 1.8+fsign(e,u); | ||
| 1094 : | if(t <= -0.6744) goto S50; | ||
| 1095 : | ignpoi = (long) (mu+s*t); | ||
| 1096 : | fk = (float)ignpoi; | ||
| 1097 : | difmuk = mu-fk; | ||
| 1098 : | /* | ||
| 1099 : | 'SUBROUTINE' F IS CALLED (KFLAG=1 FOR CORRECT RETURN) | ||
| 1100 : | */ | ||
| 1101 : | kflag = 1; | ||
| 1102 : | goto S70; | ||
| 1103 : | S60: | ||
| 1104 : | /* | ||
| 1105 : | STEP H. HAT ACCEPTANCE (E IS REPEATED ON REJECTION) | ||
| 1106 : | */ | ||
| 1107 : | if(c*fabs(u) > py*exp(px+e)-fy*exp(fx+e)) goto S50; | ||
| 1108 : | return ignpoi; | ||
| 1109 : | S70: | ||
| 1110 : | /* | ||
| 1111 : | STEP F. 'SUBROUTINE' F. CALCULATION OF PX,PY,FX,FY. | ||
| 1112 : | CASE IGNPOI .LT. 10 USES FACTORIALS FROM TABLE FACT | ||
| 1113 : | */ | ||
| 1114 : | if(ignpoi >= 10) goto S80; | ||
| 1115 : | px = -mu; | ||
| 1116 : | py = pow(mu,(double)ignpoi)/ *(fact+ignpoi); | ||
| 1117 : | goto S110; | ||
| 1118 : | S80: | ||
| 1119 : | /* | ||
| 1120 : | CASE IGNPOI .GE. 10 USES POLYNOMIAL APPROXIMATION | ||
| 1121 : | A0-A7 FOR ACCURACY WHEN ADVISABLE | ||
| 1122 : | .8333333E-1=1./12. .3989423=(2*PI)**(-.5) | ||
| 1123 : | */ | ||
| 1124 : | del = 8.333333E-2/fk; | ||
| 1125 : | del -= (4.8*del*del*del); | ||
| 1126 : | v = difmuk/fk; | ||
| 1127 : | if(fabs(v) <= 0.25) goto S90; | ||
| 1128 : | px = fk*log(1.0+v)-difmuk-del; | ||
| 1129 : | goto S100; | ||
| 1130 : | S90: | ||
| 1131 : | px = fk*v*v*(((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v+a0)-del; | ||
| 1132 : | S100: | ||
| 1133 : | py = 0.3989423/sqrt(fk); | ||
| 1134 : | S110: | ||
| 1135 : | x = (0.5-difmuk)/s; | ||
| 1136 : | xx = x*x; | ||
| 1137 : | fx = -0.5*xx; | ||
| 1138 : | fy = omega*(((c3*xx+c2)*xx+c1)*xx+c0); | ||
| 1139 : | if(kflag <= 0) goto S40; | ||
| 1140 : | goto S60; | ||
| 1141 : | S120: | ||
| 1142 : | /* | ||
| 1143 : | C A S E B. (START NEW TABLE AND CALCULATE P0 IF NECESSARY) | ||
| 1144 : | */ | ||
| 1145 : | muprev = 0.0; | ||
| 1146 : | if(mu == muold) goto S130; | ||
| 1147 : | muold = mu; | ||
| 1148 : | m = max(1L,(long) (mu)); | ||
| 1149 : | l = 0; | ||
| 1150 : | p = exp(-mu); | ||
| 1151 : | q = p0 = p; | ||
| 1152 : | S130: | ||
| 1153 : | /* | ||
| 1154 : | STEP U. UNIFORM SAMPLE FOR INVERSION METHOD | ||
| 1155 : | */ | ||
| 1156 : | u = ranf(); | ||
| 1157 : | ignpoi = 0; | ||
| 1158 : | if(u <= p0) return ignpoi; | ||
| 1159 : | /* | ||
| 1160 : | STEP T. TABLE COMPARISON UNTIL THE END PP(L) OF THE | ||
| 1161 : | PP-TABLE OF CUMULATIVE POISSON PROBABILITIES | ||
| 1162 : | (0.458=PP(9) FOR MU=10) | ||
| 1163 : | */ | ||
| 1164 : | if(l == 0) goto S150; | ||
| 1165 : | j = 1; | ||
| 1166 : | if(u > 0.458) j = min(l,m); | ||
| 1167 : | for(k=j; k<=l; k++) { | ||
| 1168 : | if(u <= *(pp+k-1)) goto S180; | ||
| 1169 : | } | ||
| 1170 : | if(l == 35) goto S130; | ||
| 1171 : | S150: | ||
| 1172 : | /* | ||
| 1173 : | STEP C. CREATION OF NEW POISSON PROBABILITIES P | ||
| 1174 : | AND THEIR CUMULATIVES Q=PP(K) | ||
| 1175 : | */ | ||
| 1176 : | l += 1; | ||
| 1177 : | for(k=l; k<=35; k++) { | ||
| 1178 : | p = p*mu/(float)k; | ||
| 1179 : | q += p; | ||
| 1180 : | *(pp+k-1) = q; | ||
| 1181 : | if(u <= q) goto S170; | ||
| 1182 : | } | ||
| 1183 : | l = 35; | ||
| 1184 : | goto S130; | ||
| 1185 : | S170: | ||
| 1186 : | l = k; | ||
| 1187 : | S180: | ||
| 1188 : | ignpoi = k; | ||
| 1189 : | return ignpoi; | ||
| 1190 : | } | ||
| 1191 : | long ignuin(long low,long high) | ||
| 1192 : | /* | ||
| 1193 : | ********************************************************************** | ||
| 1194 : | long ignuin(long low,long high) | ||
| 1195 : | GeNerate Uniform INteger | ||
| 1196 : | Function | ||
| 1197 : | Generates an integer uniformly distributed between LOW and HIGH. | ||
| 1198 : | Arguments | ||
| 1199 : | low --> Low bound (inclusive) on integer value to be generated | ||
| 1200 : | high --> High bound (inclusive) on integer value to be generated | ||
| 1201 : | Note | ||
| 1202 : | If (HIGH-LOW) > 2,147,483,561 prints error message on * unit and | ||
| 1203 : | stops the program. | ||
| 1204 : | ********************************************************************** | ||
| 1205 : | IGNLGI generates integers between 1 and 2147483562 | ||
| 1206 : | MAXNUM is 1 less than maximum generable value | ||
| 1207 : | */ | ||
| 1208 : | { | ||
| 1209 : | #define maxnum 2147483561L | ||
| 1210 : | static long ignuin,ign,maxnow,range,ranp1; | ||
| 1211 : | |||
| 1212 : | if(!(low > high)) goto S10; | ||
| 1213 : | fputs(" low > high in ignuin - ABORT\n",stderr); | ||
| 1214 : | PyErr_SetString (PyExc_ValueError, "Described above."); | ||
| 1215 : | return 0; | ||
| 1216 : | |||
| 1217 : | S10: | ||
| 1218 : | range = high-low; | ||
| 1219 : | if(!(range > maxnum)) goto S20; | ||
| 1220 : | fputs(" high - low too large in ignuin - ABORT\n",stderr); | ||
| 1221 : | PyErr_SetString (PyExc_ValueError, "Described above."); | ||
| 1222 : | return 0; | ||
| 1223 : | |||
| 1224 : | S20: | ||
| 1225 : | if(!(low == high)) goto S30; | ||
| 1226 : | ignuin = low; | ||
| 1227 : | return ignuin; | ||
| 1228 : | |||
| 1229 : | S30: | ||
| 1230 : | /* | ||
| 1231 : | Number to be generated should be in range 0..RANGE | ||
| 1232 : | Set MAXNOW so that the number of integers in 0..MAXNOW is an | ||
| 1233 : | integral multiple of the number in 0..RANGE | ||
| 1234 : | */ | ||
| 1235 : | ranp1 = range+1; | ||
| 1236 : | maxnow = maxnum/ranp1*ranp1; | ||
| 1237 : | S40: | ||
| 1238 : | ign = ignlgi()-1; | ||
| 1239 : | if(!(ign <= maxnow)) goto S50; | ||
| 1240 : | ignuin = low+ign%ranp1; | ||
| 1241 : | return ignuin; | ||
| 1242 : | S50: | ||
| 1243 : | goto S40; | ||
| 1244 : | #undef maxnum | ||
| 1245 : | #undef err1 | ||
| 1246 : | #undef err2 | ||
| 1247 : | } | ||
| 1248 : | long lennob( char *str ) | ||
| 1249 : | /* | ||
| 1250 : | Returns the length of str ignoring trailing blanks but not | ||
| 1251 : | other white space. | ||
| 1252 : | */ | ||
| 1253 : | { | ||
| 1254 : | long i, i_nb; | ||
| 1255 : | |||
| 1256 : | for (i=0, i_nb= -1L; *(str+i); i++) | ||
| 1257 : | if ( *(str+i) != ' ' ) i_nb = i; | ||
| 1258 : | return (i_nb+1); | ||
| 1259 : | } | ||
| 1260 : | long mltmod(long a,long s,long m) | ||
| 1261 : | /* | ||
| 1262 : | ********************************************************************** | ||
| 1263 : | long mltmod(long a,long s,long m) | ||
| 1264 : | Returns (A*S) MOD M | ||
| 1265 : | This is a transcription from Pascal to Fortran of routine | ||
| 1266 : | MULtMod_Decompos from the paper | ||
| 1267 : | L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package | ||
| 1268 : | with Splitting Facilities." ACM Transactions on Mathematical | ||
| 1269 : | Software, 17:98-111 (1991) | ||
| 1270 : | Arguments | ||
| 1271 : | a, s, m --> | ||
| 1272 : | ********************************************************************** | ||
| 1273 : | */ | ||
| 1274 : | { | ||
| 1275 : | #define h 32768L | ||
| 1276 : | static long mltmod,a0,a1,k,p,q,qh,rh; | ||
| 1277 : | /* | ||
| 1278 : | H = 2**((b-2)/2) where b = 32 because we are using a 32 bit | ||
| 1279 : | machine. On a different machine recompute H | ||
| 1280 : | */ | ||
| 1281 : | if(!(a <= 0 || a >= m || s <= 0 || s >= m)) goto S10; | ||
| 1282 : | fputs(" a, m, s out of order in mltmod - ABORT!\n",stderr); | ||
| 1283 : | fprintf(stderr," a = %12ld s = %12ld m = %12ld\n",a,s,m); | ||
| 1284 : | fputs(" mltmod requires: 0 < a < m; 0 < s < m\n",stderr); | ||
| 1285 : | PyErr_SetString (PyExc_ValueError, "Described above."); | ||
| 1286 : | return 0; | ||
| 1287 : | S10: | ||
| 1288 : | if(!(a < h)) goto S20; | ||
| 1289 : | a0 = a; | ||
| 1290 : | p = 0; | ||
| 1291 : | goto S120; | ||
| 1292 : | S20: | ||
| 1293 : | a1 = a/h; | ||
| 1294 : | a0 = a-h*a1; | ||
| 1295 : | qh = m/h; | ||
| 1296 : | rh = m-h*qh; | ||
| 1297 : | if(!(a1 >= h)) goto S50; | ||
| 1298 : | a1 -= h; | ||
| 1299 : | k = s/qh; | ||
| 1300 : | p = h*(s-k*qh)-k*rh; | ||
| 1301 : | S30: | ||
| 1302 : | if(!(p < 0)) goto S40; | ||
| 1303 : | p += m; | ||
| 1304 : | goto S30; | ||
| 1305 : | S40: | ||
| 1306 : | goto S60; | ||
| 1307 : | S50: | ||
| 1308 : | p = 0; | ||
| 1309 : | S60: | ||
| 1310 : | /* | ||
| 1311 : | P = (A2*S*H)MOD M | ||
| 1312 : | */ | ||
| 1313 : | if(!(a1 != 0)) goto S90; | ||
| 1314 : | q = m/a1; | ||
| 1315 : | k = s/q; | ||
| 1316 : | p -= (k*(m-a1*q)); | ||
| 1317 : | if(p > 0) p -= m; | ||
| 1318 : | p += (a1*(s-k*q)); | ||
| 1319 : | S70: | ||
| 1320 : | if(!(p < 0)) goto S80; | ||
| 1321 : | p += m; | ||
| 1322 : | goto S70; | ||
| 1323 : | S90: | ||
| 1324 : | S80: | ||
| 1325 : | k = p/qh; | ||
| 1326 : | /* | ||
| 1327 : | P = ((A2*H + A1)*S)MOD M | ||
| 1328 : | */ | ||
| 1329 : | p = h*(p-k*qh)-k*rh; | ||
| 1330 : | S100: | ||
| 1331 : | if(!(p < 0)) goto S110; | ||
| 1332 : | p += m; | ||
| 1333 : | goto S100; | ||
| 1334 : | S120: | ||
| 1335 : | S110: | ||
| 1336 : | if(!(a0 != 0)) goto S150; | ||
| 1337 : | /* | ||
| 1338 : | P = ((A2*H + A1)*H*S)MOD M | ||
| 1339 : | */ | ||
| 1340 : | q = m/a0; | ||
| 1341 : | k = s/q; | ||
| 1342 : | p -= (k*(m-a0*q)); | ||
| 1343 : | if(p > 0) p -= m; | ||
| 1344 : | p += (a0*(s-k*q)); | ||
| 1345 : | S130: | ||
| 1346 : | if(!(p < 0)) goto S140; | ||
| 1347 : | p += m; | ||
| 1348 : | goto S130; | ||
| 1349 : | S150: | ||
| 1350 : | S140: | ||
| 1351 : | mltmod = p; | ||
| 1352 : | return mltmod; | ||
| 1353 : | #undef h | ||
| 1354 : | } | ||
| 1355 : | void phrtsd(char* phrase,long *seed1,long *seed2) | ||
| 1356 : | /* | ||
| 1357 : | ********************************************************************** | ||
| 1358 : | void phrtsd(char* phrase,long *seed1,long *seed2) | ||
| 1359 : | PHRase To SeeDs | ||
| 1360 : | |||
| 1361 : | Function | ||
| 1362 : | |||
| 1363 : | Uses a phrase (character string) to generate two seeds for the RGN | ||
| 1364 : | random number generator. | ||
| 1365 : | Arguments | ||
| 1366 : | phrase --> Phrase to be used for random number generation | ||
| 1367 : | |||
| 1368 : | seed1 <-- First seed for generator | ||
| 1369 : | |||
| 1370 : | seed2 <-- Second seed for generator | ||
| 1371 : | |||
| 1372 : | Note | ||
| 1373 : | |||
| 1374 : | Trailing blanks are eliminated before the seeds are generated. | ||
| 1375 : | Generated seed values will fall in the range 1..2^30 | ||
| 1376 : | (1..1,073,741,824) | ||
| 1377 : | ********************************************************************** | ||
| 1378 : | */ | ||
| 1379 : | { | ||
| 1380 : | |||
| 1381 : | static char table[] = | ||
| 1382 : | "abcdefghijklmnopqrstuvwxyz\ | ||
| 1383 : | ABCDEFGHIJKLMNOPQRSTUVWXYZ\ | ||
| 1384 : | 0123456789\ | ||
| 1385 : | !@#$%^&*()_+[];:'\\\"<>?,./"; | ||
| 1386 : | |||
| 1387 : | long ix; | ||
| 1388 : | |||
| 1389 : | static long twop30 = 1073741824L; | ||
| 1390 : | static long shift[5] = { | ||
| 1391 : | 1L,64L,4096L,262144L,16777216L | ||
| 1392 : | }; | ||
| 1393 : | static long i,ichr,j,lphr,values[5]; | ||
| 1394 : | extern long lennob(char *str); | ||
| 1395 : | |||
| 1396 : | *seed1 = 1234567890L; | ||
| 1397 : | *seed2 = 123456789L; | ||
| 1398 : | lphr = lennob(phrase); | ||
| 1399 : | if(lphr < 1) return; | ||
| 1400 : | for(i=0; i<=(lphr-1); i++) { | ||
| 1401 : | for (ix=0; table[ix]; ix++) if (*(phrase+i) == table[ix]) break; | ||
| 1402 : | if (!table[ix]) ix = 0; | ||
| 1403 : | ichr = ix % 64; | ||
| 1404 : | if(ichr == 0) ichr = 63; | ||
| 1405 : | for(j=1; j<=5; j++) { | ||
| 1406 : | *(values+j-1) = ichr-j; | ||
| 1407 : | if(*(values+j-1) < 1) *(values+j-1) += 63; | ||
| 1408 : | } | ||
| 1409 : | for(j=1; j<=5; j++) { | ||
| 1410 : | *seed1 = ( *seed1+*(shift+j-1)**(values+j-1) ) % twop30; | ||
| 1411 : | *seed2 = ( *seed2+*(shift+j-1)**(values+6-j-1) ) % twop30; | ||
| 1412 : | } | ||
| 1413 : | } | ||
| 1414 : | #undef twop30 | ||
| 1415 : | } | ||
| 1416 : | double ranf(void) | ||
| 1417 : | /* | ||
| 1418 : | ********************************************************************** | ||
| 1419 : | double ranf(void) | ||
| 1420 : | RANDom number generator as a Function | ||
| 1421 : | Returns a random floating point number from a uniform distribution | ||
| 1422 : | over 0 - 1 (endpoints of this interval are not returned) using the | ||
| 1423 : | current generator | ||
| 1424 : | This is a transcription from Pascal to Fortran of routine | ||
| 1425 : | Uniform_01 from the paper | ||
| 1426 : | L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package | ||
| 1427 : | with Splitting Facilities." ACM Transactions on Mathematical | ||
| 1428 : | Software, 17:98-111 (1991) | ||
| 1429 : | ********************************************************************** | ||
| 1430 : | */ | ||
| 1431 : | { | ||
| 1432 : | static double ranf_; | ||
| 1433 : | /* | ||
| 1434 : | 4.656613057E-10 is 1/M1 M1 is set in a data statement in IGNLGI | ||
| 1435 : | and is currently 2147483563. If M1 changes, change this also. | ||
| 1436 : | */ | ||
| 1437 : | ranf_ = ignlgi()*4.656613057E-10; | ||
| 1438 : | return ranf_; | ||
| 1439 : | } | ||
| 1440 : | void setgmn(float *meanv,float *covm,long p,float *parm) | ||
| 1441 : | /* | ||
| 1442 : | ********************************************************************** | ||
| 1443 : | void setgmn(float *meanv,float *covm,long p,float *parm) | ||
| 1444 : | SET Generate Multivariate Normal random deviate | ||
| 1445 : | Function | ||
| 1446 : | Places P, MEANV, and the Cholesky factoriztion of COVM | ||
| 1447 : | in GENMN. | ||
| 1448 : | Arguments | ||
| 1449 : | meanv --> Mean vector of multivariate normal distribution. | ||
| 1450 : | covm <--> (Input) Covariance matrix of the multivariate | ||
| 1451 : | normal distribution | ||
| 1452 : | (Output) Destroyed on output | ||
| 1453 : | p --> Dimension of the normal, or length of MEANV. | ||
| 1454 : | parm <-- Array of parameters needed to generate multivariate norma | ||
| 1455 : | deviates (P, MEANV and Cholesky decomposition of | ||
| 1456 : | COVM). | ||
| 1457 : | 1 : 1 - P | ||
| 1458 : | 2 : P + 1 - MEANV | ||
| 1459 : | P+2 : P*(P+3)/2 + 1 - Cholesky decomposition of COVM | ||
| 1460 : | Needed dimension is (p*(p+3)/2 + 1) | ||
| 1461 : | ********************************************************************** | ||
| 1462 : | */ | ||
| 1463 : | { | ||
| 1464 : | extern void spofa(float *a,long lda,long n,long *info); | ||
| 1465 : | static long T1; | ||
| 1466 : | static long i,icount,info,j,D2,D3,D4,D5; | ||
| 1467 : | T1 = p*(p+3)/2+1; | ||
| 1468 : | /* | ||
| 1469 : | TEST THE INPUT | ||
| 1470 : | */ | ||
| 1471 : | if(!(p <= 0)) goto S10; | ||
| 1472 : | fputs("P nonpositive in SETGMN\n",stderr); | ||
| 1473 : | fprintf(stderr,"Value of P: %12ld\n",p); | ||
| 1474 : | PyErr_SetString (PyExc_ValueError, "Described above."); | ||
| 1475 : | return; | ||
| 1476 : | S10: | ||
| 1477 : | *parm = p; | ||
| 1478 : | /* | ||
| 1479 : | PUT P AND MEANV INTO PARM | ||
| 1480 : | */ | ||
| 1481 : | for(i=2,D2=1,D3=(p+1-i+D2)/D2; D3>0; D3--,i+=D2) *(parm+i-1) = *(meanv+i-2); | ||
| 1482 : | /* | ||
| 1483 : | Cholesky decomposition to find A s.t. trans(A)*(A) = COVM | ||
| 1484 : | */ | ||
| 1485 : | spofa(covm,p,p,&info); | ||
| 1486 : | if(!(info != 0)) goto S30; | ||
| 1487 : | fputs(" COVM not positive definite in SETGMN\n",stderr); | ||
| 1488 : | PyErr_SetString (PyExc_ValueError, "Described above."); | ||
| 1489 : | return; | ||
| 1490 : | S30: | ||
| 1491 : | icount = p+1; | ||
| 1492 : | /* | ||
| 1493 : | PUT UPPER HALF OF A, WHICH IS NOW THE CHOLESKY FACTOR, INTO PARM | ||
| 1494 : | COVM(1,1) = PARM(P+2) | ||
| 1495 : | COVM(1,2) = PARM(P+3) | ||
| 1496 : | : | ||
| 1497 : | COVM(1,P) = PARM(2P+1) | ||
| 1498 : | COVM(2,2) = PARM(2P+2) ... | ||
| 1499 : | */ | ||
| 1500 : | for(i=1,D4=1,D5=(p-i+D4)/D4; D5>0; D5--,i+=D4) { | ||
| 1501 : | for(j=i-1; j<p; j++) { | ||
| 1502 : | icount += 1; | ||
| 1503 : | *(parm+icount-1) = *(covm+i-1+j*p); | ||
| 1504 : | } | ||
| 1505 : | } | ||
| 1506 : | } | ||
| 1507 : | float sexpo(void) | ||
| 1508 : | /* | ||
| 1509 : | ********************************************************************** | ||
| 1510 : | |||
| 1511 : | |||
| 1512 : | (STANDARD-) E X P O N E N T I A L DISTRIBUTION | ||
| 1513 : | |||
| 1514 : | |||
| 1515 : | ********************************************************************** | ||
| 1516 : | ********************************************************************** | ||
| 1517 : | |||
| 1518 : | FOR DETAILS SEE: | ||
| 1519 : | |||
| 1520 : | AHRENS, J.H. AND DIETER, U. | ||
| 1521 : | COMPUTER METHODS FOR SAMPLING FROM THE | ||
| 1522 : | EXPONENTIAL AND NORMAL DISTRIBUTIONS. | ||
| 1523 : | COMM. ACM, 15,10 (OCT. 1972), 873 - 882. | ||
| 1524 : | |||
| 1525 : | ALL STATEMENT NUMBERS CORRESPOND TO THE STEPS OF ALGORITHM | ||
| 1526 : | 'SA' IN THE ABOVE PAPER (SLIGHTLY MODIFIED IMPLEMENTATION) | ||
| 1527 : | |||
| 1528 : | Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of | ||
| 1529 : | SUNIF. The argument IR thus goes away. | ||
| 1530 : | |||
| 1531 : | ********************************************************************** | ||
| 1532 : | Q(N) = SUM(ALOG(2.0)**K/K!) K=1,..,N , THE HIGHEST N | ||
| 1533 : | (HERE 8) IS DETERMINED BY Q(N)=1.0 WITHIN STANDARD PRECISION | ||
| 1534 : | */ | ||
| 1535 : | { | ||
| 1536 : | static float q[8] = { | ||
| 1537 : | 0.6931472,0.9333737,0.9888778,0.9984959,0.9998293,0.9999833,0.9999986,1.0 | ||
| 1538 : | }; | ||
| 1539 : | static long i; | ||
| 1540 : | static float sexpo,a,u,ustar,umin; | ||
| 1541 : | static float *q1 = q; | ||
| 1542 : | a = 0.0; | ||
| 1543 : | u = ranf(); | ||
| 1544 : | goto S30; | ||
| 1545 : | S20: | ||
| 1546 : | a += *q1; | ||
| 1547 : | S30: | ||
| 1548 : | u += u; | ||
| 1549 : | if(u <= 1.0) goto S20; | ||
| 1550 : | u -= 1.0; | ||
| 1551 : | if(u > *q1) goto S60; | ||
| 1552 : | sexpo = a+u; | ||
| 1553 : | return sexpo; | ||
| 1554 : | S60: | ||
| 1555 : | i = 1; | ||
| 1556 : | ustar = ranf(); | ||
| 1557 : | umin = ustar; | ||
| 1558 : | S70: | ||
| 1559 : | ustar = ranf(); | ||
| 1560 : | if(ustar < umin) umin = ustar; | ||
| 1561 : | i += 1; | ||
| 1562 : | if(u > *(q+i-1)) goto S70; | ||
| 1563 : | sexpo = a+umin**q1; | ||
| 1564 : | return sexpo; | ||
| 1565 : | } | ||
| 1566 : | float sgamma(float a) | ||
| 1567 : | /* | ||
| 1568 : | ********************************************************************** | ||
| 1569 : | |||
| 1570 : | |||
| 1571 : | (STANDARD-) G A M M A DISTRIBUTION | ||
| 1572 : | |||
| 1573 : | |||
| 1574 : | ********************************************************************** | ||
| 1575 : | ********************************************************************** | ||
| 1576 : | |||
| 1577 : | PARAMETER A >= 1.0 ! | ||
| 1578 : | |||
| 1579 : | ********************************************************************** | ||
| 1580 : | |||
| 1581 : | FOR DETAILS SEE: | ||
| 1582 : | |||
| 1583 : | AHRENS, J.H. AND DIETER, U. | ||
| 1584 : | GENERATING GAMMA VARIATES BY A | ||
| 1585 : | MODIFIED REJECTION TECHNIQUE. | ||
| 1586 : | COMM. ACM, 25,1 (JAN. 1982), 47 - 54. | ||
| 1587 : | |||
| 1588 : | STEP NUMBERS CORRESPOND TO ALGORITHM 'GD' IN THE ABOVE PAPER | ||
| 1589 : | (STRAIGHTFORWARD IMPLEMENTATION) | ||
| 1590 : | |||
| 1591 : | Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of | ||
| 1592 : | SUNIF. The argument IR thus goes away. | ||
| 1593 : | |||
| 1594 : | ********************************************************************** | ||
| 1595 : | |||
| 1596 : | PARAMETER 0.0 < A < 1.0 ! | ||
| 1597 : | |||
| 1598 : | ********************************************************************** | ||
| 1599 : | |||
| 1600 : | FOR DETAILS SEE: | ||
| 1601 : | |||
| 1602 : | AHRENS, J.H. AND DIETER, U. | ||
| 1603 : | COMPUTER METHODS FOR SAMPLING FROM GAMMA, | ||
| 1604 : | BETA, POISSON AND BINOMIAL DISTRIBUTIONS. | ||
| 1605 : | COMPUTING, 12 (1974), 223 - 246. | ||
| 1606 : | |||
| 1607 : | (ADAPTED IMPLEMENTATION OF ALGORITHM 'GS' IN THE ABOVE PAPER) | ||
| 1608 : | |||
| 1609 : | ********************************************************************** | ||
| 1610 : | INPUT: A =PARAMETER (MEAN) OF THE STANDARD GAMMA DISTRIBUTION | ||
| 1611 : | OUTPUT: SGAMMA = SAMPLE FROM THE GAMMA-(A)-DISTRIBUTION | ||
| 1612 : | COEFFICIENTS Q(K) - FOR Q0 = SUM(Q(K)*A**(-K)) | ||
| 1613 : | COEFFICIENTS A(K) - FOR Q = Q0+(T*T/2)*SUM(A(K)*V**K) | ||
| 1614 : | COEFFICIENTS E(K) - FOR EXP(Q)-1 = SUM(E(K)*Q**K) | ||
| 1615 : | PREVIOUS A PRE-SET TO ZERO - AA IS A', AAA IS A" | ||
| 1616 : | SQRT32 IS THE SQUAREROOT OF 32 = 5.656854249492380 | ||
| 1617 : | */ | ||
| 1618 : | { | ||
| 1619 : | extern float fsign( float num, float sign ); | ||
| 1620 : | static float q1 = 4.166669E-2; | ||
| 1621 : | static float q2 = 2.083148E-2; | ||
| 1622 : | static float q3 = 8.01191E-3; | ||
| 1623 : | static float q4 = 1.44121E-3; | ||
| 1624 : | static float q5 = -7.388E-5; | ||
| 1625 : | static float q6 = 2.4511E-4; | ||
| 1626 : | static float q7 = 2.424E-4; | ||
| 1627 : | static float a1 = 0.3333333; | ||
| 1628 : | static float a2 = -0.250003; | ||
| 1629 : | static float a3 = 0.2000062; | ||
| 1630 : | static float a4 = -0.1662921; | ||
| 1631 : | static float a5 = 0.1423657; | ||
| 1632 : | static float a6 = -0.1367177; | ||
| 1633 : | static float a7 = 0.1233795; | ||
| 1634 : | static float e1 = 1.0; | ||
| 1635 : | static float e2 = 0.4999897; | ||
| 1636 : | static float e3 = 0.166829; | ||
| 1637 : | static float e4 = 4.07753E-2; | ||
| 1638 : | static float e5 = 1.0293E-2; | ||
| 1639 : | static float aa = 0.0; | ||
| 1640 : | static float aaa = 0.0; | ||
| 1641 : | static float sqrt32 = 5.656854; | ||
| 1642 : | static float sgamma,s2,s,d,t,x,u,r,q0,b,si,c,v,q,e,w,p; | ||
| 1643 : | if(a == aa) goto S10; | ||
| 1644 : | if(a < 1.0) goto S120; | ||
| 1645 : | /* | ||
| 1646 : | STEP 1: RECALCULATIONS OF S2,S,D IF A HAS CHANGED | ||
| 1647 : | */ | ||
| 1648 : | aa = a; | ||
| 1649 : | s2 = a-0.5; | ||
| 1650 : | s = sqrt(s2); | ||
| 1651 : | d = sqrt32-12.0*s; | ||
| 1652 : | S10: | ||
| 1653 : | /* | ||
| 1654 : | STEP 2: T=STANDARD NORMAL DEVIATE, | ||
| 1655 : | X=(S,1/2)-NORMAL DEVIATE. | ||
| 1656 : | IMMEDIATE ACCEPTANCE (I) | ||
| 1657 : | */ | ||
| 1658 : | t = snorm(); | ||
| 1659 : | x = s+0.5*t; | ||
| 1660 : | sgamma = x*x; | ||
| 1661 : | if(t >= 0.0) return sgamma; | ||
| 1662 : | /* | ||
| 1663 : | STEP 3: U= 0,1 -UNIFORM SAMPLE. SQUEEZE ACCEPTANCE (S) | ||
| 1664 : | */ | ||
| 1665 : | u = ranf(); | ||
| 1666 : | if(d*u <= t*t*t) return sgamma; | ||
| 1667 : | /* | ||
| 1668 : | STEP 4: RECALCULATIONS OF Q0,B,SI,C IF NECESSARY | ||
| 1669 : | */ | ||
| 1670 : | if(a == aaa) goto S40; | ||
| 1671 : | aaa = a; | ||
| 1672 : | r = 1.0/ a; | ||
| 1673 : | q0 = ((((((q7*r+q6)*r+q5)*r+q4)*r+q3)*r+q2)*r+q1)*r; | ||
| 1674 : | /* | ||
| 1675 : | APPROXIMATION DEPENDING ON SIZE OF PARAMETER A | ||
| 1676 : | THE CONSTANTS IN THE EXPRESSIONS FOR B, SI AND | ||
| 1677 : | C WERE ESTABLISHED BY NUMERICAL EXPERIMENTS | ||
| 1678 : | */ | ||
| 1679 : | if(a <= 3.686) goto S30; | ||
| 1680 : | if(a <= 13.022) goto S20; | ||
| 1681 : | /* | ||
| 1682 : | CASE 3: A .GT. 13.022 | ||
| 1683 : | */ | ||
| 1684 : | b = 1.77; | ||
| 1685 : | si = 0.75; | ||
| 1686 : | c = 0.1515/s; | ||
| 1687 : | goto S40; | ||
| 1688 : | S20: | ||
| 1689 : | /* | ||
| 1690 : | CASE 2: 3.686 .LT. A .LE. 13.022 | ||
| 1691 : | */ | ||
| 1692 : | b = 1.654+7.6E-3*s2; | ||
| 1693 : | si = 1.68/s+0.275; | ||
| 1694 : | c = 6.2E-2/s+2.4E-2; | ||
| 1695 : | goto S40; | ||
| 1696 : | S30: | ||
| 1697 : | /* | ||
| 1698 : | CASE 1: A .LE. 3.686 | ||
| 1699 : | */ | ||
| 1700 : | b = 0.463+s+0.178*s2; | ||
| 1701 : | si = 1.235; | ||
| 1702 : | c = 0.195/s-7.9E-2+1.6E-1*s; | ||
| 1703 : | S40: | ||
| 1704 : | /* | ||
| 1705 : | STEP 5: NO QUOTIENT TEST IF X NOT POSITIVE | ||
| 1706 : | */ | ||
| 1707 : | if(x <= 0.0) goto S70; | ||
| 1708 : | /* | ||
| 1709 : | STEP 6: CALCULATION OF V AND QUOTIENT Q | ||
| 1710 : | */ | ||
| 1711 : | v = t/(s+s); | ||
| 1712 : | if(fabs(v) <= 0.25) goto S50; | ||
| 1713 : | q = q0-s*t+0.25*t*t+(s2+s2)*log(1.0+v); | ||
| 1714 : | goto S60; | ||
| 1715 : | S50: | ||
| 1716 : | q = q0+0.5*t*t*((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v; | ||
| 1717 : | S60: | ||
| 1718 : | /* | ||
| 1719 : | STEP 7: QUOTIENT ACCEPTANCE (Q) | ||
| 1720 : | */ | ||
| 1721 : | if(log(1.0-u) <= q) return sgamma; | ||
| 1722 : | S70: | ||
| 1723 : | /* | ||
| 1724 : | STEP 8: E=STANDARD EXPONENTIAL DEVIATE | ||
| 1725 : | U= 0,1 -UNIFORM DEVIATE | ||
| 1726 : | T=(B,SI)-DOUBLE EXPONENTIAL (LAPLACE) SAMPLE | ||
| 1727 : | */ | ||
| 1728 : | e = sexpo(); | ||
| 1729 : | u = ranf(); | ||
| 1730 : | u += (u-1.0); | ||
| 1731 : | t = b+fsign(si*e,u); | ||
| 1732 : | /* | ||
| 1733 : | STEP 9: REJECTION IF T .LT. TAU(1) = -.71874483771719 | ||
| 1734 : | */ | ||
| 1735 : | if(t < -0.7187449) goto S70; | ||
| 1736 : | /* | ||
| 1737 : | STEP 10: CALCULATION OF V AND QUOTIENT Q | ||
| 1738 : | */ | ||
| 1739 : | v = t/(s+s); | ||
| 1740 : | if(fabs(v) <= 0.25) goto S80; | ||
| 1741 : | q = q0-s*t+0.25*t*t+(s2+s2)*log(1.0+v); | ||
| 1742 : | goto S90; | ||
| 1743 : | S80: | ||
| 1744 : | q = q0+0.5*t*t*((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v; | ||
| 1745 : | S90: | ||
| 1746 : | /* | ||
| 1747 : | STEP 11: HAT ACCEPTANCE (H) (IF Q NOT POSITIVE GO TO STEP 8) | ||
| 1748 : | */ | ||
| 1749 : | if(q <= 0.0) goto S70; | ||
| 1750 : | if(q <= 0.5) goto S100; | ||
| 1751 : | w = exp(q)-1.0; | ||
| 1752 : | goto S110; | ||
| 1753 : | S100: | ||
| 1754 : | w = ((((e5*q+e4)*q+e3)*q+e2)*q+e1)*q; | ||
| 1755 : | S110: | ||
| 1756 : | /* | ||
| 1757 : | IF T IS REJECTED, SAMPLE AGAIN AT STEP 8 | ||
| 1758 : | */ | ||
| 1759 : | if(c*fabs(u) > w*exp(e-0.5*t*t)) goto S70; | ||
| 1760 : | x = s+0.5*t; | ||
| 1761 : | sgamma = x*x; | ||
| 1762 : | return sgamma; | ||
| 1763 : | S120: | ||
| 1764 : | /* | ||
| 1765 : | ALTERNATE METHOD FOR PARAMETERS A BELOW 1 (.3678794=EXP(-1.)) | ||
| 1766 : | */ | ||
| 1767 : | aa = 0.0; | ||
| 1768 : | b = 1.0+0.3678794*a; | ||
| 1769 : | S130: | ||
| 1770 : | p = b*ranf(); | ||
| 1771 : | if(p >= 1.0) goto S140; | ||
| 1772 : | sgamma = exp(log(p)/ a); | ||
| 1773 : | if(sexpo() < sgamma) goto S130; | ||
| 1774 : | return sgamma; | ||
| 1775 : | S140: | ||
| 1776 : | sgamma = -log((b-p)/ a); | ||
| 1777 : | if(sexpo() < (1.0-a)*log(sgamma)) goto S130; | ||
| 1778 : | return sgamma; | ||
| 1779 : | } | ||
| 1780 : | float snorm(void) | ||
| 1781 : | /* | ||
| 1782 : | ********************************************************************** | ||
| 1783 : | |||
| 1784 : | |||
| 1785 : | (STANDARD-) N O R M A L DISTRIBUTION | ||
| 1786 : | |||
| 1787 : | |||
| 1788 : | ********************************************************************** | ||
| 1789 : | ********************************************************************** | ||
| 1790 : | |||
| 1791 : | FOR DETAILS SEE: | ||
| 1792 : | |||
| 1793 : | AHRENS, J.H. AND DIETER, U. | ||
| 1794 : | EXTENSIONS OF FORSYTHE'S METHOD FOR RANDOM | ||
| 1795 : | SAMPLING FROM THE NORMAL DISTRIBUTION. | ||
| 1796 : | MATH. COMPUT., 27,124 (OCT. 1973), 927 - 937. | ||
| 1797 : | |||
| 1798 : | ALL STATEMENT NUMBERS CORRESPOND TO THE STEPS OF ALGORITHM 'FL' | ||
| 1799 : | (M=5) IN THE ABOVE PAPER (SLIGHTLY MODIFIED IMPLEMENTATION) | ||
| 1800 : | |||
| 1801 : | Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of | ||
| 1802 : | SUNIF. The argument IR thus goes away. | ||
| 1803 : | |||
| 1804 : | ********************************************************************** | ||
| 1805 : | THE DEFINITIONS OF THE CONSTANTS A(K), D(K), T(K) AND | ||
| 1806 : | H(K) ARE ACCORDING TO THE ABOVEMENTIONED ARTICLE | ||
| 1807 : | */ | ||
| 1808 : | { | ||
| 1809 : | static float a[32] = { | ||
| 1810 : | 0.0,3.917609E-2,7.841241E-2,0.11777,0.1573107,0.1970991,0.2372021,0.2776904, | ||
| 1811 : | 0.3186394,0.36013,0.4022501,0.4450965,0.4887764,0.5334097,0.5791322, | ||
| 1812 : | 0.626099,0.6744898,0.7245144,0.7764218,0.8305109,0.8871466,0.9467818, | ||
| 1813 : | 1.00999,1.077516,1.150349,1.229859,1.318011,1.417797,1.534121,1.67594, | ||
| 1814 : | 1.862732,2.153875 | ||
| 1815 : | }; | ||
| 1816 : | static float d[31] = { | ||
| 1817 : | 0.0,0.0,0.0,0.0,0.0,0.2636843,0.2425085,0.2255674,0.2116342,0.1999243, | ||
| 1818 : | 0.1899108,0.1812252,0.1736014,0.1668419,0.1607967,0.1553497,0.1504094, | ||
| 1819 : | 0.1459026,0.14177,0.1379632,0.1344418,0.1311722,0.128126,0.1252791, | ||
| 1820 : | 0.1226109,0.1201036,0.1177417,0.1155119,0.1134023,0.1114027,0.1095039 | ||
| 1821 : | }; | ||
| 1822 : | static float t[31] = { | ||
| 1823 : | 7.673828E-4,2.30687E-3,3.860618E-3,5.438454E-3,7.0507E-3,8.708396E-3, | ||
| 1824 : | 1.042357E-2,1.220953E-2,1.408125E-2,1.605579E-2,1.81529E-2,2.039573E-2, | ||
| 1825 : | 2.281177E-2,2.543407E-2,2.830296E-2,3.146822E-2,3.499233E-2,3.895483E-2, | ||
| 1826 : | 4.345878E-2,4.864035E-2,5.468334E-2,6.184222E-2,7.047983E-2,8.113195E-2, | ||
| 1827 : | 9.462444E-2,0.1123001,0.136498,0.1716886,0.2276241,0.330498,0.5847031 | ||
| 1828 : | }; | ||
| 1829 : | static float h[31] = { | ||
| 1830 : | 3.920617E-2,3.932705E-2,3.951E-2,3.975703E-2,4.007093E-2,4.045533E-2, | ||
| 1831 : | 4.091481E-2,4.145507E-2,4.208311E-2,4.280748E-2,4.363863E-2,4.458932E-2, | ||
| 1832 : | 4.567523E-2,4.691571E-2,4.833487E-2,4.996298E-2,5.183859E-2,5.401138E-2, | ||
| 1833 : | 5.654656E-2,5.95313E-2,6.308489E-2,6.737503E-2,7.264544E-2,7.926471E-2, | ||
| 1834 : | 8.781922E-2,9.930398E-2,0.11556,0.1404344,0.1836142,0.2790016,0.7010474 | ||
| 1835 : | }; | ||
| 1836 : | static long i; | ||
| 1837 : | static float snorm,u,s,ustar,aa,w,y,tt; | ||
| 1838 : | u = ranf(); | ||
| 1839 : | s = 0.0; | ||
| 1840 : | if(u > 0.5) s = 1.0; | ||
| 1841 : | u += (u-s); | ||
| 1842 : | u = 32.0*u; | ||
| 1843 : | i = (long) (u); | ||
| 1844 : | if(i == 32) i = 31; | ||
| 1845 : | if(i == 0) goto S100; | ||
| 1846 : | /* | ||
| 1847 : | START CENTER | ||
| 1848 : | */ | ||
| 1849 : | ustar = u-(float)i; | ||
| 1850 : | aa = *(a+i-1); | ||
| 1851 : | S40: | ||
| 1852 : | if(ustar <= *(t+i-1)) goto S60; | ||
| 1853 : | w = (ustar-*(t+i-1))**(h+i-1); | ||
| 1854 : | S50: | ||
| 1855 : | /* | ||
| 1856 : | EXIT (BOTH CASES) | ||
| 1857 : | */ | ||
| 1858 : | y = aa+w; | ||
| 1859 : | snorm = y; | ||
| 1860 : | if(s == 1.0) snorm = -y; | ||
| 1861 : | return snorm; | ||
| 1862 : | S60: | ||
| 1863 : | /* | ||
| 1864 : | CENTER CONTINUED | ||
| 1865 : | */ | ||
| 1866 : | u = ranf(); | ||
| 1867 : | w = u*(*(a+i)-aa); | ||
| 1868 : | tt = (0.5*w+aa)*w; | ||
| 1869 : | goto S80; | ||
| 1870 : | S70: | ||
| 1871 : | tt = u; | ||
| 1872 : | ustar = ranf(); | ||
| 1873 : | S80: | ||
| 1874 : | if(ustar > tt) goto S50; | ||
| 1875 : | u = ranf(); | ||
| 1876 : | if(ustar >= u) goto S70; | ||
| 1877 : | ustar = ranf(); | ||
| 1878 : | goto S40; | ||
| 1879 : | S100: | ||
| 1880 : | /* | ||
| 1881 : | START TAIL | ||
| 1882 : | */ | ||
| 1883 : | i = 6; | ||
| 1884 : | aa = *(a+31); | ||
| 1885 : | goto S120; | ||
| 1886 : | S110: | ||
| 1887 : | aa += *(d+i-1); | ||
| 1888 : | i += 1; | ||
| 1889 : | S120: | ||
| 1890 : | u += u; | ||
| 1891 : | if(u < 1.0) goto S110; | ||
| 1892 : | u -= 1.0; | ||
| 1893 : | S140: | ||
| 1894 : | w = u**(d+i-1); | ||
| 1895 : | tt = (0.5*w+aa)*w; | ||
| 1896 : | goto S160; | ||
| 1897 : | S150: | ||
| 1898 : | tt = u; | ||
| 1899 : | S160: | ||
| 1900 : | ustar = ranf(); | ||
| 1901 : | if(ustar > tt) goto S50; | ||
| 1902 : | u = ranf(); | ||
| 1903 : | if(ustar >= u) goto S150; | ||
| 1904 : | u = ranf(); | ||
| 1905 : | goto S140; | ||
| 1906 : | } | ||
| 1907 : | float fsign( float num, float sign ) | ||
| 1908 : | /* Transfers sign of argument sign to argument num */ | ||
| 1909 : | { | ||
| 1910 : | if ( ( sign>0.0f && num<0.0f ) || ( sign<0.0f && num>0.0f ) ) | ||
| 1911 : | return -num; | ||
| 1912 : | else return num; | ||
| 1913 : | } | ||
| 1914 : | /************************************************************************ | ||
| 1915 : | FTNSTOP: | ||
| 1916 : | Prints msg to standard error and then exits | ||
| 1917 : | ************************************************************************/ | ||
| 1918 : | void ftnstop(char* msg) | ||
| 1919 : | /* msg - error message */ | ||
| 1920 : | { | ||
| 1921 : | if (msg != NULL) fprintf(stderr,"%s\n",msg); | ||
| 1922 : | PyErr_SetString (PyExc_RuntimeError, "Described above."); | ||
| 1923 : | return; | ||
| 1924 : | } |
| MCS Webmaster | ViewVC Help |
| Powered by ViewVC 1.0.3 |