[2696] | 1 | MODULE KPP_ROOT_Random |
---|
| 2 | ! A module for random number generation from the following distributions: |
---|
| 3 | ! |
---|
| 4 | ! Distribution Function/subroutine name |
---|
| 5 | ! |
---|
| 6 | ! Normal (Gaussian) random_normal |
---|
| 7 | ! Gamma random_gamma |
---|
| 8 | ! Chi-squared random_chisq |
---|
| 9 | ! Exponential random_exponential |
---|
| 10 | ! Weibull random_Weibull |
---|
| 11 | ! Beta random_beta |
---|
| 12 | ! t random_t |
---|
| 13 | ! Multivariate normal random_mvnorm |
---|
| 14 | ! Generalized inverse Gaussian random_inv_gauss |
---|
| 15 | ! Poisson random_Poisson |
---|
| 16 | ! Binomial random_binomial1 * |
---|
| 17 | ! random_binomial2 * |
---|
| 18 | ! Negative binomial random_neg_binomial |
---|
| 19 | ! von Mises random_von_Mises |
---|
| 20 | ! Cauchy random_Cauchy |
---|
| 21 | ! |
---|
| 22 | ! Generate a random ordering of the integers 1 .. N |
---|
| 23 | ! random_order |
---|
| 24 | ! Initialize (seed) the uniform random number generator for ANY compiler |
---|
| 25 | ! seed_random_number |
---|
| 26 | |
---|
| 27 | ! Lognormal - see note below. |
---|
| 28 | |
---|
| 29 | ! ** Two functions are provided for the binomial distribution. |
---|
| 30 | ! If the parameter values remain constant, it is recommended that the |
---|
| 31 | ! first function is used (random_binomial1). If one or both of the |
---|
| 32 | ! parameters change, use the second function (random_binomial2). |
---|
| 33 | |
---|
| 34 | ! The compilers own random number generator, SUBROUTINE RANDOM_NUMBER(r), |
---|
| 35 | ! is used to provide a source of uniformly distributed random numbers. |
---|
| 36 | |
---|
| 37 | ! N.B. At this stage, only one random number is generated at each call to |
---|
| 38 | ! one of the functions above. |
---|
| 39 | |
---|
| 40 | ! The module uses the following functions which are included here: |
---|
| 41 | ! bin_prob to calculate a single binomial probability |
---|
| 42 | ! lngamma to calculate the logarithm to base e of the gamma function |
---|
| 43 | |
---|
| 44 | ! Some of the code is adapted from Dagpunar's book: |
---|
| 45 | ! Dagpunar, J. 'Principles of random variate generation' |
---|
| 46 | ! Clarendon Press, Oxford, 1988. ISBN 0-19-852202-9 |
---|
| 47 | ! |
---|
| 48 | ! In most of Dagpunar's routines, there is a test to see whether the value |
---|
| 49 | ! of one or two floating-point parameters has changed since the last call. |
---|
| 50 | ! These tests have been replaced by using a logical variable FIRST. |
---|
| 51 | ! This should be set to .TRUE. on the first call using new values of the |
---|
| 52 | ! parameters, and .FALSE. if the parameter values are the same as for the |
---|
| 53 | ! previous call. |
---|
| 54 | |
---|
| 55 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
| 56 | ! Lognormal distribution |
---|
| 57 | ! If X has a lognormal distribution, then log(X) is normally distributed. |
---|
| 58 | ! Here the logarithm is the natural logarithm, that is to base e, sometimes |
---|
| 59 | ! denoted as ln. To generate random variates from this distribution, generate |
---|
| 60 | ! a random deviate from the normal distribution with mean and variance equal |
---|
| 61 | ! to the mean and variance of the logarithms of X, then take its exponential. |
---|
| 62 | |
---|
| 63 | ! Relationship between the mean & variance of log(X) and the mean & variance |
---|
| 64 | ! of X, when X has a lognormal distribution. |
---|
| 65 | ! Let m = mean of log(X), and s^2 = variance of log(X) |
---|
| 66 | ! Then |
---|
| 67 | ! mean of X = exp(m + 0.5s^2) |
---|
| 68 | ! variance of X = (mean(X))^2.[exp(s^2) - 1] |
---|
| 69 | |
---|
| 70 | ! In the reverse direction (rarely used) |
---|
| 71 | ! variance of log(X) = log[1 + var(X)/(mean(X))^2] |
---|
| 72 | ! mean of log(X) = log(mean(X) - 0.5var(log(X)) |
---|
| 73 | |
---|
| 74 | ! N.B. The above formulae relate to population parameters; they will only be |
---|
| 75 | ! approximate if applied to sample values. |
---|
| 76 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
| 77 | |
---|
| 78 | ! Version 1.13, 2 October 2000 |
---|
| 79 | ! Changes from version 1.01 |
---|
| 80 | ! 1. The random_order, random_Poisson & random_binomial routines have been |
---|
| 81 | ! replaced with more efficient routines. |
---|
| 82 | ! 2. A routine, seed_random_number, has been added to seed the uniform random |
---|
| 83 | ! number generator. This requires input of the required number of seeds |
---|
| 84 | ! for the particular compiler from a specified I/O unit such as a keyboard. |
---|
| 85 | ! 3. Made compatible with Lahey's ELF90. |
---|
| 86 | ! 4. Marsaglia & Tsang algorithm used for random_gamma when shape parameter > 1. |
---|
| 87 | ! 5. INTENT for array f corrected in random_mvnorm. |
---|
| 88 | |
---|
| 89 | ! Author: Alan Miller |
---|
| 90 | ! CSIRO Division of Mathematical & Information Sciences |
---|
| 91 | ! Private Bag 10, Clayton South MDC |
---|
| 92 | ! Clayton 3169, Victoria, Australia |
---|
| 93 | ! Phone: (+61) 3 9545-8016 Fax: (+61) 3 9545-8080 |
---|
| 94 | ! e-mail: amiller @ bigpond.net.au |
---|
| 95 | |
---|
| 96 | IMPLICIT NONE |
---|
| 97 | REAL, PRIVATE :: zero = 0.0, half = 0.5, one = 1.0, two = 2.0, & |
---|
| 98 | vsmall = TINY(1.0), vlarge = HUGE(1.0) |
---|
| 99 | PRIVATE :: integral |
---|
| 100 | INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60) |
---|
| 101 | |
---|
| 102 | |
---|
| 103 | CONTAINS |
---|
| 104 | |
---|
| 105 | |
---|
| 106 | FUNCTION random_normal() RESULT(fn_val) |
---|
| 107 | |
---|
| 108 | ! Adapted from the following Fortran 77 code |
---|
| 109 | ! ALGORITHM 712, COLLECTED ALGORITHMS FROM ACM. |
---|
| 110 | ! THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, |
---|
| 111 | ! VOL. 18, NO. 4, DECEMBER, 1992, PP. 434-435. |
---|
| 112 | |
---|
| 113 | ! The function random_normal() returns a normally distributed pseudo-random |
---|
| 114 | ! number with zero mean and unit variance. |
---|
| 115 | |
---|
| 116 | ! The algorithm uses the ratio of uniforms method of A.J. Kinderman |
---|
| 117 | ! and J.F. Monahan augmented with quadratic bounding curves. |
---|
| 118 | |
---|
| 119 | REAL :: fn_val |
---|
| 120 | |
---|
| 121 | ! Local variables |
---|
| 122 | REAL :: s = 0.449871, t = -0.386595, a = 0.19600, b = 0.25472, & |
---|
| 123 | r1 = 0.27597, r2 = 0.27846, u, v, x, y, q |
---|
| 124 | |
---|
| 125 | ! Generate P = (u,v) uniform in rectangle enclosing acceptance region |
---|
| 126 | |
---|
| 127 | DO |
---|
| 128 | CALL RANDOM_NUMBER(u) |
---|
| 129 | CALL RANDOM_NUMBER(v) |
---|
| 130 | v = 1.7156 * (v - half) |
---|
| 131 | |
---|
| 132 | ! Evaluate the quadratic form |
---|
| 133 | x = u - s |
---|
| 134 | y = ABS(v) - t |
---|
| 135 | q = x**2 + y*(a*y - b*x) |
---|
| 136 | |
---|
| 137 | ! Accept P if inside inner ellipse |
---|
| 138 | IF (q < r1) EXIT |
---|
| 139 | ! Reject P if outside outer ellipse |
---|
| 140 | IF (q > r2) CYCLE |
---|
| 141 | ! Reject P if outside acceptance region |
---|
| 142 | IF (v**2 < -4.0*LOG(u)*u**2) EXIT |
---|
| 143 | END DO |
---|
| 144 | |
---|
| 145 | ! Return ratio of P's coordinates as the normal deviate |
---|
| 146 | fn_val = v/u |
---|
| 147 | RETURN |
---|
| 148 | |
---|
| 149 | END FUNCTION random_normal |
---|
| 150 | |
---|
| 151 | |
---|
| 152 | |
---|
| 153 | FUNCTION random_gamma(s, first) RESULT(fn_val) |
---|
| 154 | |
---|
| 155 | ! Adapted from Fortran 77 code from the book: |
---|
| 156 | ! Dagpunar, J. 'Principles of random variate generation' |
---|
| 157 | ! Clarendon Press, Oxford, 1988. ISBN 0-19-852202-9 |
---|
| 158 | |
---|
| 159 | ! FUNCTION GENERATES A RANDOM GAMMA VARIATE. |
---|
| 160 | ! CALLS EITHER random_gamma1 (S > 1.0) |
---|
| 161 | ! OR random_exponential (S = 1.0) |
---|
| 162 | ! OR random_gamma2 (S < 1.0). |
---|
| 163 | |
---|
| 164 | ! S = SHAPE PARAMETER OF DISTRIBUTION (0 < REAL). |
---|
| 165 | |
---|
| 166 | REAL, INTENT(IN) :: s |
---|
| 167 | LOGICAL, INTENT(IN) :: first |
---|
| 168 | REAL :: fn_val |
---|
| 169 | |
---|
| 170 | IF (s <= zero) THEN |
---|
| 171 | WRITE(*, *) 'SHAPE PARAMETER VALUE MUST BE POSITIVE' |
---|
| 172 | STOP |
---|
| 173 | END IF |
---|
| 174 | |
---|
| 175 | IF (s > one) THEN |
---|
| 176 | fn_val = random_gamma1(s, first) |
---|
| 177 | ELSE IF (s < one) THEN |
---|
| 178 | fn_val = random_gamma2(s, first) |
---|
| 179 | ELSE |
---|
| 180 | fn_val = random_exponential() |
---|
| 181 | END IF |
---|
| 182 | |
---|
| 183 | RETURN |
---|
| 184 | END FUNCTION random_gamma |
---|
| 185 | |
---|
| 186 | |
---|
| 187 | |
---|
| 188 | FUNCTION random_gamma1(s, first) RESULT(fn_val) |
---|
| 189 | |
---|
| 190 | ! Uses the algorithm in |
---|
| 191 | ! Marsaglia, G. and Tsang, W.W. (2000) `A simple method for generating |
---|
| 192 | ! gamma variables', Trans. om Math. Software (TOMS), vol.26(3), pp.363-372. |
---|
| 193 | |
---|
| 194 | ! Generates a random gamma deviate for shape parameter s >= 1. |
---|
| 195 | |
---|
| 196 | REAL, INTENT(IN) :: s |
---|
| 197 | LOGICAL, INTENT(IN) :: first |
---|
| 198 | REAL :: fn_val |
---|
| 199 | |
---|
| 200 | ! Local variables |
---|
| 201 | REAL, SAVE :: c, d |
---|
| 202 | REAL :: u, v, x |
---|
| 203 | |
---|
| 204 | IF (first) THEN |
---|
| 205 | d = s - one/3. |
---|
| 206 | c = one/SQRT(9.0*d) |
---|
| 207 | END IF |
---|
| 208 | |
---|
| 209 | ! Start of main loop |
---|
| 210 | DO |
---|
| 211 | |
---|
| 212 | ! Generate v = (1+cx)^3 where x is random normal; repeat if v <= 0. |
---|
| 213 | |
---|
| 214 | DO |
---|
| 215 | x = random_normal() |
---|
| 216 | v = (one + c*x)**3 |
---|
| 217 | IF (v > zero) EXIT |
---|
| 218 | END DO |
---|
| 219 | |
---|
| 220 | ! Generate uniform variable U |
---|
| 221 | |
---|
| 222 | CALL RANDOM_NUMBER(u) |
---|
| 223 | IF (u < one - 0.0331*x**4) THEN |
---|
| 224 | fn_val = d*v |
---|
| 225 | EXIT |
---|
| 226 | ELSE IF (LOG(u) < half*x**2 + d*(one - v + LOG(v))) THEN |
---|
| 227 | fn_val = d*v |
---|
| 228 | EXIT |
---|
| 229 | END IF |
---|
| 230 | END DO |
---|
| 231 | |
---|
| 232 | RETURN |
---|
| 233 | END FUNCTION random_gamma1 |
---|
| 234 | |
---|
| 235 | |
---|
| 236 | |
---|
| 237 | FUNCTION random_gamma2(s, first) RESULT(fn_val) |
---|
| 238 | |
---|
| 239 | ! Adapted from Fortran 77 code from the book: |
---|
| 240 | ! Dagpunar, J. 'Principles of random variate generation' |
---|
| 241 | ! Clarendon Press, Oxford, 1988. ISBN 0-19-852202-9 |
---|
| 242 | |
---|
| 243 | ! FUNCTION GENERATES A RANDOM VARIATE IN [0,INFINITY) FROM |
---|
| 244 | ! A GAMMA DISTRIBUTION WITH DENSITY PROPORTIONAL TO |
---|
| 245 | ! GAMMA2**(S-1) * EXP(-GAMMA2), |
---|
| 246 | ! USING A SWITCHING METHOD. |
---|
| 247 | |
---|
| 248 | ! S = SHAPE PARAMETER OF DISTRIBUTION |
---|
| 249 | ! (REAL < 1.0) |
---|
| 250 | |
---|
| 251 | REAL, INTENT(IN) :: s |
---|
| 252 | LOGICAL, INTENT(IN) :: first |
---|
| 253 | REAL :: fn_val |
---|
| 254 | |
---|
| 255 | ! Local variables |
---|
| 256 | REAL :: r, x, w |
---|
| 257 | REAL, SAVE :: a, p, c, uf, vr, d |
---|
| 258 | |
---|
| 259 | IF (s <= zero .OR. s >= one) THEN |
---|
| 260 | WRITE(*, *) 'SHAPE PARAMETER VALUE OUTSIDE PERMITTED RANGE' |
---|
| 261 | STOP |
---|
| 262 | END IF |
---|
| 263 | |
---|
| 264 | IF (first) THEN ! Initialization, if necessary |
---|
| 265 | a = one - s |
---|
| 266 | p = a/(a + s*EXP(-a)) |
---|
| 267 | IF (s < vsmall) THEN |
---|
| 268 | WRITE(*, *) 'SHAPE PARAMETER VALUE TOO SMALL' |
---|
| 269 | STOP |
---|
| 270 | END IF |
---|
| 271 | c = one/s |
---|
| 272 | uf = p*(vsmall/a)**s |
---|
| 273 | vr = one - vsmall |
---|
| 274 | d = a*LOG(a) |
---|
| 275 | END IF |
---|
| 276 | |
---|
| 277 | DO |
---|
| 278 | CALL RANDOM_NUMBER(r) |
---|
| 279 | IF (r >= vr) THEN |
---|
| 280 | CYCLE |
---|
| 281 | ELSE IF (r > p) THEN |
---|
| 282 | x = a - LOG((one - r)/(one - p)) |
---|
| 283 | w = a*LOG(x)-d |
---|
| 284 | ELSE IF (r > uf) THEN |
---|
| 285 | x = a*(r/p)**c |
---|
| 286 | w = x |
---|
| 287 | ELSE |
---|
| 288 | fn_val = zero |
---|
| 289 | RETURN |
---|
| 290 | END IF |
---|
| 291 | |
---|
| 292 | CALL RANDOM_NUMBER(r) |
---|
| 293 | IF (one-r <= w .AND. r > zero) THEN |
---|
| 294 | IF (r*(w + one) >= one) CYCLE |
---|
| 295 | IF (-LOG(r) <= w) CYCLE |
---|
| 296 | END IF |
---|
| 297 | EXIT |
---|
| 298 | END DO |
---|
| 299 | |
---|
| 300 | fn_val = x |
---|
| 301 | RETURN |
---|
| 302 | |
---|
| 303 | END FUNCTION random_gamma2 |
---|
| 304 | |
---|
| 305 | |
---|
| 306 | |
---|
| 307 | FUNCTION random_chisq(ndf, first) RESULT(fn_val) |
---|
| 308 | |
---|
| 309 | ! Generates a random variate from the chi-squared distribution with |
---|
| 310 | ! ndf degrees of freedom |
---|
| 311 | |
---|
| 312 | INTEGER, INTENT(IN) :: ndf |
---|
| 313 | LOGICAL, INTENT(IN) :: first |
---|
| 314 | REAL :: fn_val |
---|
| 315 | |
---|
| 316 | fn_val = two * random_gamma(half*ndf, first) |
---|
| 317 | RETURN |
---|
| 318 | |
---|
| 319 | END FUNCTION random_chisq |
---|
| 320 | |
---|
| 321 | |
---|
| 322 | |
---|
| 323 | FUNCTION random_exponential() RESULT(fn_val) |
---|
| 324 | |
---|
| 325 | ! Adapted from Fortran 77 code from the book: |
---|
| 326 | ! Dagpunar, J. 'Principles of random variate generation' |
---|
| 327 | ! Clarendon Press, Oxford, 1988. ISBN 0-19-852202-9 |
---|
| 328 | |
---|
| 329 | ! FUNCTION GENERATES A RANDOM VARIATE IN [0,INFINITY) FROM |
---|
| 330 | ! A NEGATIVE EXPONENTIAL DlSTRIBUTION WlTH DENSITY PROPORTIONAL |
---|
| 331 | ! TO EXP(-random_exponential), USING INVERSION. |
---|
| 332 | |
---|
| 333 | REAL :: fn_val |
---|
| 334 | |
---|
| 335 | ! Local variable |
---|
| 336 | REAL :: r |
---|
| 337 | |
---|
| 338 | DO |
---|
| 339 | CALL RANDOM_NUMBER(r) |
---|
| 340 | IF (r > zero) EXIT |
---|
| 341 | END DO |
---|
| 342 | |
---|
| 343 | fn_val = -LOG(r) |
---|
| 344 | RETURN |
---|
| 345 | |
---|
| 346 | END FUNCTION random_exponential |
---|
| 347 | |
---|
| 348 | |
---|
| 349 | |
---|
| 350 | FUNCTION random_Weibull(a) RESULT(fn_val) |
---|
| 351 | |
---|
| 352 | ! Generates a random variate from the Weibull distribution with |
---|
| 353 | ! probability density: |
---|
| 354 | ! a |
---|
| 355 | ! a-1 -x |
---|
| 356 | ! f(x) = a.x e |
---|
| 357 | |
---|
| 358 | REAL, INTENT(IN) :: a |
---|
| 359 | REAL :: fn_val |
---|
| 360 | |
---|
| 361 | ! For speed, there is no checking that a is not zero or very small. |
---|
| 362 | |
---|
| 363 | fn_val = random_exponential() ** (one/a) |
---|
| 364 | RETURN |
---|
| 365 | |
---|
| 366 | END FUNCTION random_Weibull |
---|
| 367 | |
---|
| 368 | |
---|
| 369 | |
---|
| 370 | FUNCTION random_beta(aa, bb, first) RESULT(fn_val) |
---|
| 371 | |
---|
| 372 | ! Adapted from Fortran 77 code from the book: |
---|
| 373 | ! Dagpunar, J. 'Principles of random variate generation' |
---|
| 374 | ! Clarendon Press, Oxford, 1988. ISBN 0-19-852202-9 |
---|
| 375 | |
---|
| 376 | ! FUNCTION GENERATES A RANDOM VARIATE IN [0,1] |
---|
| 377 | ! FROM A BETA DISTRIBUTION WITH DENSITY |
---|
| 378 | ! PROPORTIONAL TO BETA**(AA-1) * (1-BETA)**(BB-1). |
---|
| 379 | ! USING CHENG'S LOG LOGISTIC METHOD. |
---|
| 380 | |
---|
| 381 | ! AA = SHAPE PARAMETER FROM DISTRIBUTION (0 < REAL) |
---|
| 382 | ! BB = SHAPE PARAMETER FROM DISTRIBUTION (0 < REAL) |
---|
| 383 | |
---|
| 384 | REAL, INTENT(IN) :: aa, bb |
---|
| 385 | LOGICAL, INTENT(IN) :: first |
---|
| 386 | REAL :: fn_val |
---|
| 387 | |
---|
| 388 | ! Local variables |
---|
| 389 | REAL, PARAMETER :: aln4 = 1.3862944 |
---|
| 390 | REAL :: a, b, g, r, s, x, y, z |
---|
| 391 | REAL, SAVE :: d, f, h, t, c |
---|
| 392 | LOGICAL, SAVE :: swap |
---|
| 393 | |
---|
| 394 | IF (aa <= zero .OR. bb <= zero) THEN |
---|
| 395 | WRITE(*, *) 'IMPERMISSIBLE SHAPE PARAMETER VALUE(S)' |
---|
| 396 | STOP |
---|
| 397 | END IF |
---|
| 398 | |
---|
| 399 | IF (first) THEN ! Initialization, if necessary |
---|
| 400 | a = aa |
---|
| 401 | b = bb |
---|
| 402 | swap = b > a |
---|
| 403 | IF (swap) THEN |
---|
| 404 | g = b |
---|
| 405 | b = a |
---|
| 406 | a = g |
---|
| 407 | END IF |
---|
| 408 | d = a/b |
---|
| 409 | f = a+b |
---|
| 410 | IF (b > one) THEN |
---|
| 411 | h = SQRT((two*a*b - f)/(f - two)) |
---|
| 412 | t = one |
---|
| 413 | ELSE |
---|
| 414 | h = b |
---|
| 415 | t = one/(one + (a/(vlarge*b))**b) |
---|
| 416 | END IF |
---|
| 417 | c = a+h |
---|
| 418 | END IF |
---|
| 419 | |
---|
| 420 | DO |
---|
| 421 | CALL RANDOM_NUMBER(r) |
---|
| 422 | CALL RANDOM_NUMBER(x) |
---|
| 423 | s = r*r*x |
---|
| 424 | IF (r < vsmall .OR. s <= zero) CYCLE |
---|
| 425 | IF (r < t) THEN |
---|
| 426 | x = LOG(r/(one - r))/h |
---|
| 427 | y = d*EXP(x) |
---|
| 428 | z = c*x + f*LOG((one + d)/(one + y)) - aln4 |
---|
| 429 | IF (s - one > z) THEN |
---|
| 430 | IF (s - s*z > one) CYCLE |
---|
| 431 | IF (LOG(s) > z) CYCLE |
---|
| 432 | END IF |
---|
| 433 | fn_val = y/(one + y) |
---|
| 434 | ELSE |
---|
| 435 | IF (4.0*s > (one + one/d)**f) CYCLE |
---|
| 436 | fn_val = one |
---|
| 437 | END IF |
---|
| 438 | EXIT |
---|
| 439 | END DO |
---|
| 440 | |
---|
| 441 | IF (swap) fn_val = one - fn_val |
---|
| 442 | RETURN |
---|
| 443 | END FUNCTION random_beta |
---|
| 444 | |
---|
| 445 | |
---|
| 446 | |
---|
| 447 | FUNCTION random_t(m) RESULT(fn_val) |
---|
| 448 | |
---|
| 449 | ! Adapted from Fortran 77 code from the book: |
---|
| 450 | ! Dagpunar, J. 'Principles of random variate generation' |
---|
| 451 | ! Clarendon Press, Oxford, 1988. ISBN 0-19-852202-9 |
---|
| 452 | |
---|
| 453 | ! FUNCTION GENERATES A RANDOM VARIATE FROM A |
---|
| 454 | ! T DISTRIBUTION USING KINDERMAN AND MONAHAN'S RATIO METHOD. |
---|
| 455 | |
---|
| 456 | ! M = DEGREES OF FREEDOM OF DISTRIBUTION |
---|
| 457 | ! (1 <= 1NTEGER) |
---|
| 458 | |
---|
| 459 | INTEGER, INTENT(IN) :: m |
---|
| 460 | REAL :: fn_val |
---|
| 461 | |
---|
| 462 | ! Local variables |
---|
| 463 | REAL, SAVE :: s, c, a, f, g |
---|
| 464 | REAL :: r, x, v |
---|
| 465 | |
---|
| 466 | REAL, PARAMETER :: three = 3.0, four = 4.0, quart = 0.25, & |
---|
| 467 | five = 5.0, sixteen = 16.0 |
---|
| 468 | INTEGER :: mm = 0 |
---|
| 469 | |
---|
| 470 | IF (m < 1) THEN |
---|
| 471 | WRITE(*, *) 'IMPERMISSIBLE DEGREES OF FREEDOM' |
---|
| 472 | STOP |
---|
| 473 | END IF |
---|
| 474 | |
---|
| 475 | IF (m /= mm) THEN ! Initialization, if necessary |
---|
| 476 | s = m |
---|
| 477 | c = -quart*(s + one) |
---|
| 478 | a = four/(one + one/s)**c |
---|
| 479 | f = sixteen/a |
---|
| 480 | IF (m > 1) THEN |
---|
| 481 | g = s - one |
---|
| 482 | g = ((s + one)/g)**c*SQRT((s+s)/g) |
---|
| 483 | ELSE |
---|
| 484 | g = one |
---|
| 485 | END IF |
---|
| 486 | mm = m |
---|
| 487 | END IF |
---|
| 488 | |
---|
| 489 | DO |
---|
| 490 | CALL RANDOM_NUMBER(r) |
---|
| 491 | IF (r <= zero) CYCLE |
---|
| 492 | CALL RANDOM_NUMBER(v) |
---|
| 493 | x = (two*v - one)*g/r |
---|
| 494 | v = x*x |
---|
| 495 | IF (v > five - a*r) THEN |
---|
| 496 | IF (m >= 1 .AND. r*(v + three) > f) CYCLE |
---|
| 497 | IF (r > (one + v/s)**c) CYCLE |
---|
| 498 | END IF |
---|
| 499 | EXIT |
---|
| 500 | END DO |
---|
| 501 | |
---|
| 502 | fn_val = x |
---|
| 503 | RETURN |
---|
| 504 | END FUNCTION random_t |
---|
| 505 | |
---|
| 506 | |
---|
| 507 | |
---|
| 508 | SUBROUTINE random_mvnorm(n, h, d, f, first, x, ier) |
---|
| 509 | |
---|
| 510 | ! Adapted from Fortran 77 code from the book: |
---|
| 511 | ! Dagpunar, J. 'Principles of random variate generation' |
---|
| 512 | ! Clarendon Press, Oxford, 1988. ISBN 0-19-852202-9 |
---|
| 513 | |
---|
| 514 | ! N.B. An extra argument, ier, has been added to Dagpunar's routine |
---|
| 515 | |
---|
| 516 | ! SUBROUTINE GENERATES AN N VARIATE RANDOM NORMAL |
---|
| 517 | ! VECTOR USING A CHOLESKY DECOMPOSITION. |
---|
| 518 | |
---|
| 519 | ! ARGUMENTS: |
---|
| 520 | ! N = NUMBER OF VARIATES IN VECTOR |
---|
| 521 | ! (INPUT,INTEGER >= 1) |
---|
| 522 | ! H(J) = J'TH ELEMENT OF VECTOR OF MEANS |
---|
| 523 | ! (INPUT,REAL) |
---|
| 524 | ! X(J) = J'TH ELEMENT OF DELIVERED VECTOR |
---|
| 525 | ! (OUTPUT,REAL) |
---|
| 526 | ! |
---|
| 527 | ! D(J*(J-1)/2+I) = (I,J)'TH ELEMENT OF VARIANCE MATRIX (J> = I) |
---|
| 528 | ! (INPUT,REAL) |
---|
| 529 | ! F((J-1)*(2*N-J)/2+I) = (I,J)'TH ELEMENT OF LOWER TRIANGULAR |
---|
| 530 | ! DECOMPOSITION OF VARIANCE MATRIX (J <= I) |
---|
| 531 | ! (OUTPUT,REAL) |
---|
| 532 | |
---|
| 533 | ! FIRST = .TRUE. IF THIS IS THE FIRST CALL OF THE ROUTINE |
---|
| 534 | ! OR IF THE DISTRIBUTION HAS CHANGED SINCE THE LAST CALL OF THE ROUTINE. |
---|
| 535 | ! OTHERWISE SET TO .FALSE. |
---|
| 536 | ! (INPUT,LOGICAL) |
---|
| 537 | |
---|
| 538 | ! ier = 1 if the input covariance matrix is not +ve definite |
---|
| 539 | ! = 0 otherwise |
---|
| 540 | |
---|
| 541 | INTEGER, INTENT(IN) :: n |
---|
| 542 | REAL, INTENT(IN) :: h(:), d(:) ! d(n*(n+1)/2) |
---|
| 543 | REAL, INTENT(IN OUT) :: f(:) ! f(n*(n+1)/2) |
---|
| 544 | REAL, INTENT(OUT) :: x(:) |
---|
| 545 | LOGICAL, INTENT(IN) :: first |
---|
| 546 | INTEGER, INTENT(OUT) :: ier |
---|
| 547 | |
---|
| 548 | ! Local variables |
---|
| 549 | INTEGER :: j, i, m |
---|
| 550 | REAL :: y, v |
---|
| 551 | INTEGER, SAVE :: n2 |
---|
| 552 | |
---|
| 553 | IF (n < 1) THEN |
---|
| 554 | WRITE(*, *) 'SIZE OF VECTOR IS NON POSITIVE' |
---|
| 555 | STOP |
---|
| 556 | END IF |
---|
| 557 | |
---|
| 558 | ier = 0 |
---|
| 559 | IF (first) THEN ! Initialization, if necessary |
---|
| 560 | n2 = 2*n |
---|
| 561 | IF (d(1) < zero) THEN |
---|
| 562 | ier = 1 |
---|
| 563 | RETURN |
---|
| 564 | END IF |
---|
| 565 | |
---|
| 566 | f(1) = SQRT(d(1)) |
---|
| 567 | y = one/f(1) |
---|
| 568 | DO j = 2,n |
---|
| 569 | f(j) = d(1+j*(j-1)/2) * y |
---|
| 570 | END DO |
---|
| 571 | |
---|
| 572 | DO i = 2,n |
---|
| 573 | v = d(i*(i-1)/2+i) |
---|
| 574 | DO m = 1,i-1 |
---|
| 575 | v = v - f((m-1)*(n2-m)/2+i)**2 |
---|
| 576 | END DO |
---|
| 577 | |
---|
| 578 | IF (v < zero) THEN |
---|
| 579 | ier = 1 |
---|
| 580 | RETURN |
---|
| 581 | END IF |
---|
| 582 | |
---|
| 583 | v = SQRT(v) |
---|
| 584 | y = one/v |
---|
| 585 | f((i-1)*(n2-i)/2+i) = v |
---|
| 586 | DO j = i+1,n |
---|
| 587 | v = d(j*(j-1)/2+i) |
---|
| 588 | DO m = 1,i-1 |
---|
| 589 | v = v - f((m-1)*(n2-m)/2+i)*f((m-1)*(n2-m)/2 + j) |
---|
| 590 | END DO ! m = 1,i-1 |
---|
| 591 | f((i-1)*(n2-i)/2 + j) = v*y |
---|
| 592 | END DO ! j = i+1,n |
---|
| 593 | END DO ! i = 2,n |
---|
| 594 | END IF |
---|
| 595 | |
---|
| 596 | x(1:n) = h(1:n) |
---|
| 597 | DO j = 1,n |
---|
| 598 | y = random_normal() |
---|
| 599 | DO i = j,n |
---|
| 600 | x(i) = x(i) + f((j-1)*(n2-j)/2 + i) * y |
---|
| 601 | END DO ! i = j,n |
---|
| 602 | END DO ! j = 1,n |
---|
| 603 | |
---|
| 604 | RETURN |
---|
| 605 | END SUBROUTINE random_mvnorm |
---|
| 606 | |
---|
| 607 | |
---|
| 608 | |
---|
| 609 | FUNCTION random_inv_gauss(h, b, first) RESULT(fn_val) |
---|
| 610 | |
---|
| 611 | ! Adapted from Fortran 77 code from the book: |
---|
| 612 | ! Dagpunar, J. 'Principles of random variate generation' |
---|
| 613 | ! Clarendon Press, Oxford, 1988. ISBN 0-19-852202-9 |
---|
| 614 | |
---|
| 615 | ! FUNCTION GENERATES A RANDOM VARIATE IN [0,INFINITY] FROM |
---|
| 616 | ! A REPARAMETERISED GENERALISED INVERSE GAUSSIAN (GIG) DISTRIBUTION |
---|
| 617 | ! WITH DENSITY PROPORTIONAL TO GIG**(H-1) * EXP(-0.5*B*(GIG+1/GIG)) |
---|
| 618 | ! USING A RATIO METHOD. |
---|
| 619 | |
---|
| 620 | ! H = PARAMETER OF DISTRIBUTION (0 <= REAL) |
---|
| 621 | ! B = PARAMETER OF DISTRIBUTION (0 < REAL) |
---|
| 622 | |
---|
| 623 | REAL, INTENT(IN) :: h, b |
---|
| 624 | LOGICAL, INTENT(IN) :: first |
---|
| 625 | REAL :: fn_val |
---|
| 626 | |
---|
| 627 | ! Local variables |
---|
| 628 | REAL :: ym, xm, r, w, r1, r2, x |
---|
| 629 | REAL, SAVE :: a, c, d, e |
---|
| 630 | REAL, PARAMETER :: quart = 0.25 |
---|
| 631 | |
---|
| 632 | IF (h < zero .OR. b <= zero) THEN |
---|
| 633 | WRITE(*, *) 'IMPERMISSIBLE DISTRIBUTION PARAMETER VALUES' |
---|
| 634 | STOP |
---|
| 635 | END IF |
---|
| 636 | |
---|
| 637 | IF (first) THEN ! Initialization, if necessary |
---|
| 638 | IF (h > quart*b*SQRT(vlarge)) THEN |
---|
| 639 | WRITE(*, *) 'THE RATIO H:B IS TOO SMALL' |
---|
| 640 | STOP |
---|
| 641 | END IF |
---|
| 642 | e = b*b |
---|
| 643 | d = h + one |
---|
| 644 | ym = (-d + SQRT(d*d + e))/b |
---|
| 645 | IF (ym < vsmall) THEN |
---|
| 646 | WRITE(*, *) 'THE VALUE OF B IS TOO SMALL' |
---|
| 647 | STOP |
---|
| 648 | END IF |
---|
| 649 | |
---|
| 650 | d = h - one |
---|
| 651 | xm = (d + SQRT(d*d + e))/b |
---|
| 652 | d = half*d |
---|
| 653 | e = -quart*b |
---|
| 654 | r = xm + one/xm |
---|
| 655 | w = xm*ym |
---|
| 656 | a = w**(-half*h) * SQRT(xm/ym) * EXP(-e*(r - ym - one/ym)) |
---|
| 657 | IF (a < vsmall) THEN |
---|
| 658 | WRITE(*, *) 'THE VALUE OF H IS TOO LARGE' |
---|
| 659 | STOP |
---|
| 660 | END IF |
---|
| 661 | c = -d*LOG(xm) - e*r |
---|
| 662 | END IF |
---|
| 663 | |
---|
| 664 | DO |
---|
| 665 | CALL RANDOM_NUMBER(r1) |
---|
| 666 | IF (r1 <= zero) CYCLE |
---|
| 667 | CALL RANDOM_NUMBER(r2) |
---|
| 668 | x = a*r2/r1 |
---|
| 669 | IF (x <= zero) CYCLE |
---|
| 670 | IF (LOG(r1) < d*LOG(x) + e*(x + one/x) + c) EXIT |
---|
| 671 | END DO |
---|
| 672 | |
---|
| 673 | fn_val = x |
---|
| 674 | |
---|
| 675 | RETURN |
---|
| 676 | END FUNCTION random_inv_gauss |
---|
| 677 | |
---|
| 678 | |
---|
| 679 | |
---|
| 680 | FUNCTION random_Poisson(mu, first) RESULT(ival) |
---|
| 681 | !********************************************************************** |
---|
| 682 | ! Translated to Fortran 90 by Alan Miller from: |
---|
| 683 | ! RANLIB |
---|
| 684 | ! |
---|
| 685 | ! Library of Fortran Routines for Random Number Generation |
---|
| 686 | ! |
---|
| 687 | ! Compiled and Written by: |
---|
| 688 | ! |
---|
| 689 | ! Barry W. Brown |
---|
| 690 | ! James Lovato |
---|
| 691 | ! |
---|
| 692 | ! Department of Biomathematics, Box 237 |
---|
| 693 | ! The University of Texas, M.D. Anderson Cancer Center |
---|
| 694 | ! 1515 Holcombe Boulevard |
---|
| 695 | ! Houston, TX 77030 |
---|
| 696 | ! |
---|
| 697 | ! This work was supported by grant CA-16672 from the National Cancer Institute. |
---|
| 698 | |
---|
| 699 | ! GENerate POIsson random deviate |
---|
| 700 | |
---|
| 701 | ! Function |
---|
| 702 | |
---|
| 703 | ! Generates a single random deviate from a Poisson distribution with mean mu. |
---|
| 704 | |
---|
| 705 | ! Arguments |
---|
| 706 | |
---|
| 707 | ! mu --> The mean of the Poisson distribution from which |
---|
| 708 | ! a random deviate is to be generated. |
---|
| 709 | ! REAL mu |
---|
| 710 | |
---|
| 711 | ! Method |
---|
| 712 | |
---|
| 713 | ! For details see: |
---|
| 714 | |
---|
| 715 | ! Ahrens, J.H. and Dieter, U. |
---|
| 716 | ! Computer Generation of Poisson Deviates |
---|
| 717 | ! From Modified Normal Distributions. |
---|
| 718 | ! ACM Trans. Math. Software, 8, 2 |
---|
| 719 | ! (June 1982),163-179 |
---|
| 720 | |
---|
| 721 | ! TABLES: COEFFICIENTS A0-A7 FOR STEP F. FACTORIALS FACT |
---|
| 722 | ! COEFFICIENTS A(K) - FOR PX = FK*V*V*SUM(A(K)*V**K)-DEL |
---|
| 723 | |
---|
| 724 | ! SEPARATION OF CASES A AND B |
---|
| 725 | |
---|
| 726 | ! .. Scalar Arguments .. |
---|
| 727 | REAL, INTENT(IN) :: mu |
---|
| 728 | LOGICAL, INTENT(IN) :: first |
---|
| 729 | INTEGER :: ival |
---|
| 730 | ! .. |
---|
| 731 | ! .. Local Scalars .. |
---|
| 732 | REAL :: b1, b2, c, c0, c1, c2, c3, del, difmuk, e, fk, fx, fy, g, & |
---|
| 733 | omega, px, py, t, u, v, x, xx |
---|
| 734 | REAL, SAVE :: s, d, p, q, p0 |
---|
| 735 | INTEGER :: j, k, kflag |
---|
| 736 | LOGICAL, SAVE :: full_init |
---|
| 737 | INTEGER, SAVE :: l, m |
---|
| 738 | ! .. |
---|
| 739 | ! .. Local Arrays .. |
---|
| 740 | REAL, SAVE :: pp(35) |
---|
| 741 | ! .. |
---|
| 742 | ! .. Data statements .. |
---|
| 743 | REAL, PARAMETER :: a0 = -.5, a1 = .3333333, a2 = -.2500068, a3 = .2000118, & |
---|
| 744 | a4 = -.1661269, a5 = .1421878, a6 = -.1384794, & |
---|
| 745 | a7 = .1250060 |
---|
| 746 | |
---|
| 747 | REAL, PARAMETER :: fact(10) = (/ 1., 1., 2., 6., 24., 120., 720., 5040., & |
---|
| 748 | 40320., 362880. /) |
---|
| 749 | |
---|
| 750 | ! .. |
---|
| 751 | ! .. Executable Statements .. |
---|
| 752 | IF (mu > 10.0) THEN |
---|
| 753 | ! C A S E A. (RECALCULATION OF S, D, L IF MU HAS CHANGED) |
---|
| 754 | |
---|
| 755 | IF (first) THEN |
---|
| 756 | s = SQRT(mu) |
---|
| 757 | d = 6.0*mu*mu |
---|
| 758 | |
---|
| 759 | ! THE POISSON PROBABILITIES PK EXCEED THE DISCRETE NORMAL |
---|
| 760 | ! PROBABILITIES FK WHENEVER K >= M(MU). L=IFIX(MU-1.1484) |
---|
| 761 | ! IS AN UPPER BOUND TO M(MU) FOR ALL MU >= 10 . |
---|
| 762 | |
---|
| 763 | l = mu - 1.1484 |
---|
| 764 | full_init = .false. |
---|
| 765 | END IF |
---|
| 766 | |
---|
| 767 | |
---|
| 768 | ! STEP N. NORMAL SAMPLE - random_normal() FOR STANDARD NORMAL DEVIATE |
---|
| 769 | |
---|
| 770 | g = mu + s*random_normal() |
---|
| 771 | IF (g > 0.0) THEN |
---|
| 772 | ival = g |
---|
| 773 | |
---|
| 774 | ! STEP I. IMMEDIATE ACCEPTANCE IF ival IS LARGE ENOUGH |
---|
| 775 | |
---|
| 776 | IF (ival>=l) RETURN |
---|
| 777 | |
---|
| 778 | ! STEP S. SQUEEZE ACCEPTANCE - SAMPLE U |
---|
| 779 | |
---|
| 780 | fk = ival |
---|
| 781 | difmuk = mu - fk |
---|
| 782 | CALL RANDOM_NUMBER(u) |
---|
| 783 | IF (d*u >= difmuk*difmuk*difmuk) RETURN |
---|
| 784 | END IF |
---|
| 785 | |
---|
| 786 | ! STEP P. PREPARATIONS FOR STEPS Q AND H. |
---|
| 787 | ! (RECALCULATIONS OF PARAMETERS IF NECESSARY) |
---|
| 788 | ! .3989423=(2*PI)**(-.5) .416667E-1=1./24. .1428571=1./7. |
---|
| 789 | ! THE QUANTITIES B1, B2, C3, C2, C1, C0 ARE FOR THE HERMITE |
---|
| 790 | ! APPROXIMATIONS TO THE DISCRETE NORMAL PROBABILITIES FK. |
---|
| 791 | ! C=.1069/MU GUARANTEES MAJORIZATION BY THE 'HAT'-FUNCTION. |
---|
| 792 | |
---|
| 793 | IF (.NOT. full_init) THEN |
---|
| 794 | omega = .3989423/s |
---|
| 795 | b1 = .4166667E-1/mu |
---|
| 796 | b2 = .3*b1*b1 |
---|
| 797 | c3 = .1428571*b1*b2 |
---|
| 798 | c2 = b2 - 15.*c3 |
---|
| 799 | c1 = b1 - 6.*b2 + 45.*c3 |
---|
| 800 | c0 = 1. - b1 + 3.*b2 - 15.*c3 |
---|
| 801 | c = .1069/mu |
---|
| 802 | full_init = .true. |
---|
| 803 | END IF |
---|
| 804 | |
---|
| 805 | IF (g < 0.0) GO TO 50 |
---|
| 806 | |
---|
| 807 | ! 'SUBROUTINE' F IS CALLED (KFLAG=0 FOR CORRECT RETURN) |
---|
| 808 | |
---|
| 809 | kflag = 0 |
---|
| 810 | GO TO 70 |
---|
| 811 | |
---|
| 812 | ! STEP Q. QUOTIENT ACCEPTANCE (RARE CASE) |
---|
| 813 | |
---|
| 814 | 40 IF (fy-u*fy <= py*EXP(px-fx)) RETURN |
---|
| 815 | |
---|
| 816 | ! STEP E. EXPONENTIAL SAMPLE - random_exponential() FOR STANDARD EXPONENTIAL |
---|
| 817 | ! DEVIATE E AND SAMPLE T FROM THE LAPLACE 'HAT' |
---|
| 818 | ! (IF T <= -.6744 THEN PK < FK FOR ALL MU >= 10.) |
---|
| 819 | |
---|
| 820 | 50 e = random_exponential() |
---|
| 821 | CALL RANDOM_NUMBER(u) |
---|
| 822 | u = u + u - one |
---|
| 823 | t = 1.8 + SIGN(e, u) |
---|
| 824 | IF (t <= (-.6744)) GO TO 50 |
---|
| 825 | ival = mu + s*t |
---|
| 826 | fk = ival |
---|
| 827 | difmuk = mu - fk |
---|
| 828 | |
---|
| 829 | ! 'SUBROUTINE' F IS CALLED (KFLAG=1 FOR CORRECT RETURN) |
---|
| 830 | |
---|
| 831 | kflag = 1 |
---|
| 832 | GO TO 70 |
---|
| 833 | |
---|
| 834 | ! STEP H. HAT ACCEPTANCE (E IS REPEATED ON REJECTION) |
---|
| 835 | |
---|
| 836 | 60 IF (c*ABS(u) > py*EXP(px+e) - fy*EXP(fx+e)) GO TO 50 |
---|
| 837 | RETURN |
---|
| 838 | |
---|
| 839 | ! STEP F. 'SUBROUTINE' F. CALCULATION OF PX, PY, FX, FY. |
---|
| 840 | ! CASE ival < 10 USES FACTORIALS FROM TABLE FACT |
---|
| 841 | |
---|
| 842 | 70 IF (ival>=10) GO TO 80 |
---|
| 843 | px = -mu |
---|
| 844 | py = mu**ival/fact(ival+1) |
---|
| 845 | GO TO 110 |
---|
| 846 | |
---|
| 847 | ! CASE ival >= 10 USES POLYNOMIAL APPROXIMATION |
---|
| 848 | ! A0-A7 FOR ACCURACY WHEN ADVISABLE |
---|
| 849 | ! .8333333E-1=1./12. .3989423=(2*PI)**(-.5) |
---|
| 850 | |
---|
| 851 | 80 del = .8333333E-1/fk |
---|
| 852 | del = del - 4.8*del*del*del |
---|
| 853 | v = difmuk/fk |
---|
| 854 | IF (ABS(v)>0.25) THEN |
---|
| 855 | px = fk*LOG(one + v) - difmuk - del |
---|
| 856 | ELSE |
---|
| 857 | px = fk*v*v* (((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v+a0) - del |
---|
| 858 | END IF |
---|
| 859 | py = .3989423/SQRT(fk) |
---|
| 860 | 110 x = (half - difmuk)/s |
---|
| 861 | xx = x*x |
---|
| 862 | fx = -half*xx |
---|
| 863 | fy = omega* (((c3*xx + c2)*xx + c1)*xx + c0) |
---|
| 864 | IF (kflag <= 0) GO TO 40 |
---|
| 865 | GO TO 60 |
---|
| 866 | |
---|
| 867 | !--------------------------------------------------------------------------- |
---|
| 868 | ! C A S E B. mu < 10 |
---|
| 869 | ! START NEW TABLE AND CALCULATE P0 IF NECESSARY |
---|
| 870 | |
---|
| 871 | ELSE |
---|
| 872 | IF (first) THEN |
---|
| 873 | m = MAX(1, INT(mu)) |
---|
| 874 | l = 0 |
---|
| 875 | p = EXP(-mu) |
---|
| 876 | q = p |
---|
| 877 | p0 = p |
---|
| 878 | END IF |
---|
| 879 | |
---|
| 880 | ! STEP U. UNIFORM SAMPLE FOR INVERSION METHOD |
---|
| 881 | |
---|
| 882 | DO |
---|
| 883 | CALL RANDOM_NUMBER(u) |
---|
| 884 | ival = 0 |
---|
| 885 | IF (u <= p0) RETURN |
---|
| 886 | |
---|
| 887 | ! STEP T. TABLE COMPARISON UNTIL THE END PP(L) OF THE |
---|
| 888 | ! PP-TABLE OF CUMULATIVE POISSON PROBABILITIES |
---|
| 889 | ! (0.458=PP(9) FOR MU=10) |
---|
| 890 | |
---|
| 891 | IF (l == 0) GO TO 150 |
---|
| 892 | j = 1 |
---|
| 893 | IF (u > 0.458) j = MIN(l, m) |
---|
| 894 | DO k = j, l |
---|
| 895 | IF (u <= pp(k)) GO TO 180 |
---|
| 896 | END DO |
---|
| 897 | IF (l == 35) CYCLE |
---|
| 898 | |
---|
| 899 | ! STEP C. CREATION OF NEW POISSON PROBABILITIES P |
---|
| 900 | ! AND THEIR CUMULATIVES Q=PP(K) |
---|
| 901 | |
---|
| 902 | 150 l = l + 1 |
---|
| 903 | DO k = l, 35 |
---|
| 904 | p = p*mu / k |
---|
| 905 | q = q + p |
---|
| 906 | pp(k) = q |
---|
| 907 | IF (u <= q) GO TO 170 |
---|
| 908 | END DO |
---|
| 909 | l = 35 |
---|
| 910 | END DO |
---|
| 911 | |
---|
| 912 | 170 l = k |
---|
| 913 | 180 ival = k |
---|
| 914 | RETURN |
---|
| 915 | END IF |
---|
| 916 | |
---|
| 917 | RETURN |
---|
| 918 | END FUNCTION random_Poisson |
---|
| 919 | |
---|
| 920 | |
---|
| 921 | |
---|
| 922 | FUNCTION random_binomial1(n, p, first) RESULT(ival) |
---|
| 923 | |
---|
| 924 | ! FUNCTION GENERATES A RANDOM BINOMIAL VARIATE USING C.D.Kemp's method. |
---|
| 925 | ! This algorithm is suitable when many random variates are required |
---|
| 926 | ! with the SAME parameter values for n & p. |
---|
| 927 | |
---|
| 928 | ! P = BERNOULLI SUCCESS PROBABILITY |
---|
| 929 | ! (0 <= REAL <= 1) |
---|
| 930 | ! N = NUMBER OF BERNOULLI TRIALS |
---|
| 931 | ! (1 <= INTEGER) |
---|
| 932 | ! FIRST = .TRUE. for the first call using the current parameter values |
---|
| 933 | ! = .FALSE. if the values of (n,p) are unchanged from last call |
---|
| 934 | |
---|
| 935 | ! Reference: Kemp, C.D. (1986). `A modal method for generating binomial |
---|
| 936 | ! variables', Commun. Statist. - Theor. Meth. 15(3), 805-813. |
---|
| 937 | |
---|
| 938 | INTEGER, INTENT(IN) :: n |
---|
| 939 | REAL, INTENT(IN) :: p |
---|
| 940 | LOGICAL, INTENT(IN) :: first |
---|
| 941 | INTEGER :: ival |
---|
| 942 | |
---|
| 943 | ! Local variables |
---|
| 944 | |
---|
| 945 | INTEGER :: ru, rd |
---|
| 946 | INTEGER, SAVE :: r0 |
---|
| 947 | REAL :: u, pd, pu |
---|
| 948 | REAL, SAVE :: odds_ratio, p_r |
---|
| 949 | REAL, PARAMETER :: zero = 0.0, one = 1.0 |
---|
| 950 | |
---|
| 951 | IF (first) THEN |
---|
| 952 | r0 = (n+1)*p |
---|
| 953 | p_r = bin_prob(n, p, r0) |
---|
| 954 | odds_ratio = p / (one - p) |
---|
| 955 | END IF |
---|
| 956 | |
---|
| 957 | CALL RANDOM_NUMBER(u) |
---|
| 958 | u = u - p_r |
---|
| 959 | IF (u < zero) THEN |
---|
| 960 | ival = r0 |
---|
| 961 | RETURN |
---|
| 962 | END IF |
---|
| 963 | |
---|
| 964 | pu = p_r |
---|
| 965 | ru = r0 |
---|
| 966 | pd = p_r |
---|
| 967 | rd = r0 |
---|
| 968 | DO |
---|
| 969 | rd = rd - 1 |
---|
| 970 | IF (rd >= 0) THEN |
---|
| 971 | pd = pd * (rd+1) / (odds_ratio * (n-rd)) |
---|
| 972 | u = u - pd |
---|
| 973 | IF (u < zero) THEN |
---|
| 974 | ival = rd |
---|
| 975 | RETURN |
---|
| 976 | END IF |
---|
| 977 | END IF |
---|
| 978 | |
---|
| 979 | ru = ru + 1 |
---|
| 980 | IF (ru <= n) THEN |
---|
| 981 | pu = pu * (n-ru+1) * odds_ratio / ru |
---|
| 982 | u = u - pu |
---|
| 983 | IF (u < zero) THEN |
---|
| 984 | ival = ru |
---|
| 985 | RETURN |
---|
| 986 | END IF |
---|
| 987 | END IF |
---|
| 988 | END DO |
---|
| 989 | |
---|
| 990 | ! This point should not be reached, but just in case: |
---|
| 991 | |
---|
| 992 | ival = r0 |
---|
| 993 | RETURN |
---|
| 994 | |
---|
| 995 | END FUNCTION random_binomial1 |
---|
| 996 | |
---|
| 997 | |
---|
| 998 | |
---|
| 999 | FUNCTION bin_prob(n, p, r) RESULT(fn_val) |
---|
| 1000 | ! Calculate a binomial probability |
---|
| 1001 | |
---|
| 1002 | INTEGER, INTENT(IN) :: n, r |
---|
| 1003 | REAL, INTENT(IN) :: p |
---|
| 1004 | REAL :: fn_val |
---|
| 1005 | |
---|
| 1006 | ! Local variable |
---|
| 1007 | REAL :: one = 1.0 |
---|
| 1008 | |
---|
| 1009 | fn_val = EXP( lngamma(DBLE(n+1)) - lngamma(DBLE(r+1)) - lngamma(DBLE(n-r+1)) & |
---|
| 1010 | + r*LOG(p) + (n-r)*LOG(one - p) ) |
---|
| 1011 | RETURN |
---|
| 1012 | |
---|
| 1013 | END FUNCTION bin_prob |
---|
| 1014 | |
---|
| 1015 | |
---|
| 1016 | |
---|
| 1017 | FUNCTION lngamma(x) RESULT(fn_val) |
---|
| 1018 | |
---|
| 1019 | ! Logarithm to base e of the gamma function. |
---|
| 1020 | ! |
---|
| 1021 | ! Accurate to about 1.e-14. |
---|
| 1022 | ! Programmer: Alan Miller |
---|
| 1023 | |
---|
| 1024 | ! Latest revision of Fortran 77 version - 28 February 1988 |
---|
| 1025 | |
---|
| 1026 | REAL (dp), INTENT(IN) :: x |
---|
| 1027 | REAL (dp) :: fn_val |
---|
| 1028 | |
---|
| 1029 | ! Local variables |
---|
| 1030 | |
---|
| 1031 | REAL (dp) :: a1 = -4.166666666554424D-02, a2 = 2.430554511376954D-03, & |
---|
| 1032 | a3 = -7.685928044064347D-04, a4 = 5.660478426014386D-04, & |
---|
| 1033 | temp, arg, product, lnrt2pi = 9.189385332046727D-1, & |
---|
| 1034 | pi = 3.141592653589793D0 |
---|
| 1035 | LOGICAL :: reflect |
---|
| 1036 | |
---|
| 1037 | ! lngamma is not defined if x = 0 or a negative integer. |
---|
| 1038 | |
---|
| 1039 | IF (x > 0.d0) GO TO 10 |
---|
| 1040 | IF (x /= INT(x)) GO TO 10 |
---|
| 1041 | fn_val = 0.d0 |
---|
| 1042 | RETURN |
---|
| 1043 | |
---|
| 1044 | ! If x < 0, use the reflection formula: |
---|
| 1045 | ! gamma(x) * gamma(1-x) = pi * cosec(pi.x) |
---|
| 1046 | |
---|
| 1047 | 10 reflect = (x < 0.d0) |
---|
| 1048 | IF (reflect) THEN |
---|
| 1049 | arg = 1.d0 - x |
---|
| 1050 | ELSE |
---|
| 1051 | arg = x |
---|
| 1052 | END IF |
---|
| 1053 | |
---|
| 1054 | ! Increase the argument, if necessary, to make it > 10. |
---|
| 1055 | |
---|
| 1056 | product = 1.d0 |
---|
| 1057 | 20 IF (arg <= 10.d0) THEN |
---|
| 1058 | product = product * arg |
---|
| 1059 | arg = arg + 1.d0 |
---|
| 1060 | GO TO 20 |
---|
| 1061 | END IF |
---|
| 1062 | |
---|
| 1063 | ! Use a polynomial approximation to Stirling's formula. |
---|
| 1064 | ! N.B. The real Stirling's formula is used here, not the simpler, but less |
---|
| 1065 | ! accurate formula given by De Moivre in a letter to Stirling, which |
---|
| 1066 | ! is the one usually quoted. |
---|
| 1067 | |
---|
| 1068 | arg = arg - 0.5D0 |
---|
| 1069 | temp = 1.d0/arg**2 |
---|
| 1070 | fn_val = lnrt2pi + arg * (LOG(arg) - 1.d0 + & |
---|
| 1071 | (((a4*temp + a3)*temp + a2)*temp + a1)*temp) - LOG(product) |
---|
| 1072 | IF (reflect) THEN |
---|
| 1073 | temp = SIN(pi * x) |
---|
| 1074 | fn_val = LOG(pi/temp) - fn_val |
---|
| 1075 | END IF |
---|
| 1076 | RETURN |
---|
| 1077 | END FUNCTION lngamma |
---|
| 1078 | |
---|
| 1079 | |
---|
| 1080 | |
---|
| 1081 | FUNCTION random_binomial2(n, pp, first) RESULT(ival) |
---|
| 1082 | !********************************************************************** |
---|
| 1083 | ! Translated to Fortran 90 by Alan Miller from: |
---|
| 1084 | ! RANLIB |
---|
| 1085 | ! |
---|
| 1086 | ! Library of Fortran Routines for Random Number Generation |
---|
| 1087 | ! |
---|
| 1088 | ! Compiled and Written by: |
---|
| 1089 | ! |
---|
| 1090 | ! Barry W. Brown |
---|
| 1091 | ! James Lovato |
---|
| 1092 | ! |
---|
| 1093 | ! Department of Biomathematics, Box 237 |
---|
| 1094 | ! The University of Texas, M.D. Anderson Cancer Center |
---|
| 1095 | ! 1515 Holcombe Boulevard |
---|
| 1096 | ! Houston, TX 77030 |
---|
| 1097 | ! |
---|
| 1098 | ! This work was supported by grant CA-16672 from the National Cancer Institute. |
---|
| 1099 | |
---|
| 1100 | ! GENerate BINomial random deviate |
---|
| 1101 | |
---|
| 1102 | ! Function |
---|
| 1103 | |
---|
| 1104 | ! Generates a single random deviate from a binomial |
---|
| 1105 | ! distribution whose number of trials is N and whose |
---|
| 1106 | ! probability of an event in each trial is P. |
---|
| 1107 | |
---|
| 1108 | ! Arguments |
---|
| 1109 | |
---|
| 1110 | ! N --> The number of trials in the binomial distribution |
---|
| 1111 | ! from which a random deviate is to be generated. |
---|
| 1112 | ! INTEGER N |
---|
| 1113 | |
---|
| 1114 | ! P --> The probability of an event in each trial of the |
---|
| 1115 | ! binomial distribution from which a random deviate |
---|
| 1116 | ! is to be generated. |
---|
| 1117 | ! REAL P |
---|
| 1118 | |
---|
| 1119 | ! FIRST --> Set FIRST = .TRUE. for the first call to perform initialization |
---|
| 1120 | ! the set FIRST = .FALSE. for further calls using the same pair |
---|
| 1121 | ! of parameter values (N, P). |
---|
| 1122 | ! LOGICAL FIRST |
---|
| 1123 | |
---|
| 1124 | ! random_binomial2 <-- A random deviate yielding the number of events |
---|
| 1125 | ! from N independent trials, each of which has |
---|
| 1126 | ! a probability of event P. |
---|
| 1127 | ! INTEGER random_binomial |
---|
| 1128 | |
---|
| 1129 | ! Method |
---|
| 1130 | |
---|
| 1131 | ! This is algorithm BTPE from: |
---|
| 1132 | |
---|
| 1133 | ! Kachitvichyanukul, V. and Schmeiser, B. W. |
---|
| 1134 | ! Binomial Random Variate Generation. |
---|
| 1135 | ! Communications of the ACM, 31, 2 (February, 1988) 216. |
---|
| 1136 | |
---|
| 1137 | !********************************************************************** |
---|
| 1138 | |
---|
| 1139 | !*****DETERMINE APPROPRIATE ALGORITHM AND WHETHER SETUP IS NECESSARY |
---|
| 1140 | |
---|
| 1141 | ! .. |
---|
| 1142 | ! .. Scalar Arguments .. |
---|
| 1143 | REAL, INTENT(IN) :: pp |
---|
| 1144 | INTEGER, INTENT(IN) :: n |
---|
| 1145 | LOGICAL, INTENT(IN) :: first |
---|
| 1146 | INTEGER :: ival |
---|
| 1147 | ! .. |
---|
| 1148 | ! .. Local Scalars .. |
---|
| 1149 | REAL :: alv, amaxp, f, f1, f2, u, v, w, w2, x, x1, x2, ynorm, z, z2 |
---|
| 1150 | REAL, PARAMETER :: zero = 0.0, half = 0.5, one = 1.0 |
---|
| 1151 | INTEGER :: i, ix, ix1, k, mp |
---|
| 1152 | INTEGER, SAVE :: m |
---|
| 1153 | REAL, SAVE :: p, q, xnp, ffm, fm, xnpq, p1, xm, xl, xr, c, al, xll, & |
---|
| 1154 | xlr, p2, p3, p4, qn, r, g |
---|
| 1155 | |
---|
| 1156 | ! .. |
---|
| 1157 | ! .. Executable Statements .. |
---|
| 1158 | |
---|
| 1159 | !*****SETUP, PERFORM ONLY WHEN PARAMETERS CHANGE |
---|
| 1160 | |
---|
| 1161 | IF (first) THEN |
---|
| 1162 | p = MIN(pp, one-pp) |
---|
| 1163 | q = one - p |
---|
| 1164 | xnp = n * p |
---|
| 1165 | END IF |
---|
| 1166 | |
---|
| 1167 | IF (xnp > 30.) THEN |
---|
| 1168 | IF (first) THEN |
---|
| 1169 | ffm = xnp + p |
---|
| 1170 | m = ffm |
---|
| 1171 | fm = m |
---|
| 1172 | xnpq = xnp * q |
---|
| 1173 | p1 = INT(2.195*SQRT(xnpq) - 4.6*q) + half |
---|
| 1174 | xm = fm + half |
---|
| 1175 | xl = xm - p1 |
---|
| 1176 | xr = xm + p1 |
---|
| 1177 | c = 0.134 + 20.5 / (15.3 + fm) |
---|
| 1178 | al = (ffm-xl) / (ffm - xl*p) |
---|
| 1179 | xll = al * (one + half*al) |
---|
| 1180 | al = (xr - ffm) / (xr*q) |
---|
| 1181 | xlr = al * (one + half*al) |
---|
| 1182 | p2 = p1 * (one + c + c) |
---|
| 1183 | p3 = p2 + c / xll |
---|
| 1184 | p4 = p3 + c / xlr |
---|
| 1185 | END IF |
---|
| 1186 | |
---|
| 1187 | !*****GENERATE VARIATE, Binomial mean at least 30. |
---|
| 1188 | |
---|
| 1189 | 20 CALL RANDOM_NUMBER(u) |
---|
| 1190 | u = u * p4 |
---|
| 1191 | CALL RANDOM_NUMBER(v) |
---|
| 1192 | |
---|
| 1193 | ! TRIANGULAR REGION |
---|
| 1194 | |
---|
| 1195 | IF (u <= p1) THEN |
---|
| 1196 | ix = xm - p1 * v + u |
---|
| 1197 | GO TO 110 |
---|
| 1198 | END IF |
---|
| 1199 | |
---|
| 1200 | ! PARALLELOGRAM REGION |
---|
| 1201 | |
---|
| 1202 | IF (u <= p2) THEN |
---|
| 1203 | x = xl + (u-p1) / c |
---|
| 1204 | v = v * c + one - ABS(xm-x) / p1 |
---|
| 1205 | IF (v > one .OR. v <= zero) GO TO 20 |
---|
| 1206 | ix = x |
---|
| 1207 | ELSE |
---|
| 1208 | |
---|
| 1209 | ! LEFT TAIL |
---|
| 1210 | |
---|
| 1211 | IF (u <= p3) THEN |
---|
| 1212 | ix = xl + LOG(v) / xll |
---|
| 1213 | IF (ix < 0) GO TO 20 |
---|
| 1214 | v = v * (u-p2) * xll |
---|
| 1215 | ELSE |
---|
| 1216 | |
---|
| 1217 | ! RIGHT TAIL |
---|
| 1218 | |
---|
| 1219 | ix = xr - LOG(v) / xlr |
---|
| 1220 | IF (ix > n) GO TO 20 |
---|
| 1221 | v = v * (u-p3) * xlr |
---|
| 1222 | END IF |
---|
| 1223 | END IF |
---|
| 1224 | |
---|
| 1225 | !*****DETERMINE APPROPRIATE WAY TO PERFORM ACCEPT/REJECT TEST |
---|
| 1226 | |
---|
| 1227 | k = ABS(ix-m) |
---|
| 1228 | IF (k <= 20 .OR. k >= xnpq/2-1) THEN |
---|
| 1229 | |
---|
| 1230 | ! EXPLICIT EVALUATION |
---|
| 1231 | |
---|
| 1232 | f = one |
---|
| 1233 | r = p / q |
---|
| 1234 | g = (n+1) * r |
---|
| 1235 | IF (m < ix) THEN |
---|
| 1236 | mp = m + 1 |
---|
| 1237 | DO i = mp, ix |
---|
| 1238 | f = f * (g/i-r) |
---|
| 1239 | END DO |
---|
| 1240 | |
---|
| 1241 | ELSE IF (m > ix) THEN |
---|
| 1242 | ix1 = ix + 1 |
---|
| 1243 | DO i = ix1, m |
---|
| 1244 | f = f / (g/i-r) |
---|
| 1245 | END DO |
---|
| 1246 | END IF |
---|
| 1247 | |
---|
| 1248 | IF (v > f) THEN |
---|
| 1249 | GO TO 20 |
---|
| 1250 | ELSE |
---|
| 1251 | GO TO 110 |
---|
| 1252 | END IF |
---|
| 1253 | END IF |
---|
| 1254 | |
---|
| 1255 | ! SQUEEZING USING UPPER AND LOWER BOUNDS ON LOG(F(X)) |
---|
| 1256 | |
---|
| 1257 | amaxp = (k/xnpq) * ((k*(k/3. + .625) + .1666666666666)/xnpq + half) |
---|
| 1258 | ynorm = -k * k / (2.*xnpq) |
---|
| 1259 | alv = LOG(v) |
---|
| 1260 | IF (alv<ynorm - amaxp) GO TO 110 |
---|
| 1261 | IF (alv>ynorm + amaxp) GO TO 20 |
---|
| 1262 | |
---|
| 1263 | ! STIRLING'S (actually de Moivre's) FORMULA TO MACHINE ACCURACY FOR |
---|
| 1264 | ! THE FINAL ACCEPTANCE/REJECTION TEST |
---|
| 1265 | |
---|
| 1266 | x1 = ix + 1 |
---|
| 1267 | f1 = fm + one |
---|
| 1268 | z = n + 1 - fm |
---|
| 1269 | w = n - ix + one |
---|
| 1270 | z2 = z * z |
---|
| 1271 | x2 = x1 * x1 |
---|
| 1272 | f2 = f1 * f1 |
---|
| 1273 | w2 = w * w |
---|
| 1274 | IF (alv - (xm*LOG(f1/x1) + (n-m+half)*LOG(z/w) + (ix-m)*LOG(w*p/(x1*q)) + & |
---|
| 1275 | (13860.-(462.-(132.-(99.-140./f2)/f2)/f2)/f2)/f1/166320. + & |
---|
| 1276 | (13860.-(462.-(132.-(99.-140./z2)/z2)/z2)/z2)/z/166320. + & |
---|
| 1277 | (13860.-(462.-(132.-(99.-140./x2)/x2)/x2)/x2)/x1/166320. + & |
---|
| 1278 | (13860.-(462.-(132.-(99.-140./w2)/w2)/w2)/w2)/w/166320.) > zero) THEN |
---|
| 1279 | GO TO 20 |
---|
| 1280 | ELSE |
---|
| 1281 | GO TO 110 |
---|
| 1282 | END IF |
---|
| 1283 | |
---|
| 1284 | ELSE |
---|
| 1285 | ! INVERSE CDF LOGIC FOR MEAN LESS THAN 30 |
---|
| 1286 | IF (first) THEN |
---|
| 1287 | qn = q ** n |
---|
| 1288 | r = p / q |
---|
| 1289 | g = r * (n+1) |
---|
| 1290 | END IF |
---|
| 1291 | |
---|
| 1292 | 90 ix = 0 |
---|
| 1293 | f = qn |
---|
| 1294 | CALL RANDOM_NUMBER(u) |
---|
| 1295 | 100 IF (u >= f) THEN |
---|
| 1296 | IF (ix > 110) GO TO 90 |
---|
| 1297 | u = u - f |
---|
| 1298 | ix = ix + 1 |
---|
| 1299 | f = f * (g/ix - r) |
---|
| 1300 | GO TO 100 |
---|
| 1301 | END IF |
---|
| 1302 | END IF |
---|
| 1303 | |
---|
| 1304 | 110 IF (pp > half) ix = n - ix |
---|
| 1305 | ival = ix |
---|
| 1306 | RETURN |
---|
| 1307 | |
---|
| 1308 | END FUNCTION random_binomial2 |
---|
| 1309 | |
---|
| 1310 | |
---|
| 1311 | |
---|
| 1312 | |
---|
| 1313 | FUNCTION random_neg_binomial(sk, p) RESULT(ival) |
---|
| 1314 | |
---|
| 1315 | ! Adapted from Fortran 77 code from the book: |
---|
| 1316 | ! Dagpunar, J. 'Principles of random variate generation' |
---|
| 1317 | ! Clarendon Press, Oxford, 1988. ISBN 0-19-852202-9 |
---|
| 1318 | |
---|
| 1319 | ! FUNCTION GENERATES A RANDOM NEGATIVE BINOMIAL VARIATE USING UNSTORED |
---|
| 1320 | ! INVERSION AND/OR THE REPRODUCTIVE PROPERTY. |
---|
| 1321 | |
---|
| 1322 | ! SK = NUMBER OF FAILURES REQUIRED (Dagpunar's words!) |
---|
| 1323 | ! = the `power' parameter of the negative binomial |
---|
| 1324 | ! (0 < REAL) |
---|
| 1325 | ! P = BERNOULLI SUCCESS PROBABILITY |
---|
| 1326 | ! (0 < REAL < 1) |
---|
| 1327 | |
---|
| 1328 | ! THE PARAMETER H IS SET SO THAT UNSTORED INVERSION ONLY IS USED WHEN P <= H, |
---|
| 1329 | ! OTHERWISE A COMBINATION OF UNSTORED INVERSION AND |
---|
| 1330 | ! THE REPRODUCTIVE PROPERTY IS USED. |
---|
| 1331 | |
---|
| 1332 | REAL, INTENT(IN) :: sk, p |
---|
| 1333 | INTEGER :: ival |
---|
| 1334 | |
---|
| 1335 | ! Local variables |
---|
| 1336 | ! THE PARAMETER ULN = -LOG(MACHINE'S SMALLEST REAL NUMBER). |
---|
| 1337 | |
---|
| 1338 | REAL, PARAMETER :: h = 0.7 |
---|
| 1339 | REAL :: q, x, st, uln, v, r, s, y, g |
---|
| 1340 | INTEGER :: k, i, n |
---|
| 1341 | |
---|
| 1342 | IF (sk <= zero .OR. p <= zero .OR. p >= one) THEN |
---|
| 1343 | WRITE(*, *) 'IMPERMISSIBLE DISTRIBUTION PARAMETER VALUES' |
---|
| 1344 | STOP |
---|
| 1345 | END IF |
---|
| 1346 | |
---|
| 1347 | q = one - p |
---|
| 1348 | x = zero |
---|
| 1349 | st = sk |
---|
| 1350 | IF (p > h) THEN |
---|
| 1351 | v = one/LOG(p) |
---|
| 1352 | k = st |
---|
| 1353 | DO i = 1,k |
---|
| 1354 | DO |
---|
| 1355 | CALL RANDOM_NUMBER(r) |
---|
| 1356 | IF (r > zero) EXIT |
---|
| 1357 | END DO |
---|
| 1358 | n = v*LOG(r) |
---|
| 1359 | x = x + n |
---|
| 1360 | END DO |
---|
| 1361 | st = st - k |
---|
| 1362 | END IF |
---|
| 1363 | |
---|
| 1364 | s = zero |
---|
| 1365 | uln = -LOG(vsmall) |
---|
| 1366 | IF (st > -uln/LOG(q)) THEN |
---|
| 1367 | WRITE(*, *) ' P IS TOO LARGE FOR THIS VALUE OF SK' |
---|
| 1368 | STOP |
---|
| 1369 | END IF |
---|
| 1370 | |
---|
| 1371 | y = q**st |
---|
| 1372 | g = st |
---|
| 1373 | CALL RANDOM_NUMBER(r) |
---|
| 1374 | DO |
---|
| 1375 | IF (y > r) EXIT |
---|
| 1376 | r = r - y |
---|
| 1377 | s = s + one |
---|
| 1378 | y = y*p*g/s |
---|
| 1379 | g = g + one |
---|
| 1380 | END DO |
---|
| 1381 | |
---|
| 1382 | ival = x + s + half |
---|
| 1383 | RETURN |
---|
| 1384 | END FUNCTION random_neg_binomial |
---|
| 1385 | |
---|
| 1386 | |
---|
| 1387 | |
---|
| 1388 | FUNCTION random_von_Mises(k, first) RESULT(fn_val) |
---|
| 1389 | |
---|
| 1390 | ! Algorithm VMD from: |
---|
| 1391 | ! Dagpunar, J.S. (1990) `Sampling from the von Mises distribution via a |
---|
| 1392 | ! comparison of random numbers', J. of Appl. Statist., 17, 165-168. |
---|
| 1393 | |
---|
| 1394 | ! Fortran 90 code by Alan Miller |
---|
| 1395 | ! CSIRO Division of Mathematical & Information Sciences |
---|
| 1396 | |
---|
| 1397 | ! Arguments: |
---|
| 1398 | ! k (real) parameter of the von Mises distribution. |
---|
| 1399 | ! first (logical) set to .TRUE. the first time that the function |
---|
| 1400 | ! is called, or the first time with a new value |
---|
| 1401 | ! for k. When first = .TRUE., the function sets |
---|
| 1402 | ! up starting values and may be very much slower. |
---|
| 1403 | |
---|
| 1404 | REAL, INTENT(IN) :: k |
---|
| 1405 | LOGICAL, INTENT(IN) :: first |
---|
| 1406 | REAL :: fn_val |
---|
| 1407 | |
---|
| 1408 | ! Local variables |
---|
| 1409 | |
---|
| 1410 | INTEGER :: j, n |
---|
| 1411 | INTEGER, SAVE :: nk |
---|
| 1412 | REAL, PARAMETER :: pi = 3.14159265 |
---|
| 1413 | REAL, SAVE :: p(20), theta(0:20) |
---|
| 1414 | REAL :: sump, r, th, lambda, rlast |
---|
| 1415 | REAL (dp) :: dk |
---|
| 1416 | |
---|
| 1417 | IF (first) THEN ! Initialization, if necessary |
---|
| 1418 | IF (k < zero) THEN |
---|
| 1419 | WRITE(*, *) '** Error: argument k for random_von_Mises = ', k |
---|
| 1420 | RETURN |
---|
| 1421 | END IF |
---|
| 1422 | |
---|
| 1423 | nk = k + k + one |
---|
| 1424 | IF (nk > 20) THEN |
---|
| 1425 | WRITE(*, *) '** Error: argument k for random_von_Mises = ', k |
---|
| 1426 | RETURN |
---|
| 1427 | END IF |
---|
| 1428 | |
---|
| 1429 | dk = k |
---|
| 1430 | theta(0) = zero |
---|
| 1431 | IF (k > half) THEN |
---|
| 1432 | |
---|
| 1433 | ! Set up array p of probabilities. |
---|
| 1434 | |
---|
| 1435 | sump = zero |
---|
| 1436 | DO j = 1, nk |
---|
| 1437 | IF (j < nk) THEN |
---|
| 1438 | theta(j) = ACOS(one - j/k) |
---|
| 1439 | ELSE |
---|
| 1440 | theta(nk) = pi |
---|
| 1441 | END IF |
---|
| 1442 | |
---|
| 1443 | ! Numerical integration of e^[k.cos(x)] from theta(j-1) to theta(j) |
---|
| 1444 | |
---|
| 1445 | CALL integral(theta(j-1), theta(j), p(j), dk) |
---|
| 1446 | sump = sump + p(j) |
---|
| 1447 | END DO |
---|
| 1448 | p(1:nk) = p(1:nk) / sump |
---|
| 1449 | ELSE |
---|
| 1450 | p(1) = one |
---|
| 1451 | theta(1) = pi |
---|
| 1452 | END IF ! if k > 0.5 |
---|
| 1453 | END IF ! if first |
---|
| 1454 | |
---|
| 1455 | CALL RANDOM_NUMBER(r) |
---|
| 1456 | DO j = 1, nk |
---|
| 1457 | r = r - p(j) |
---|
| 1458 | IF (r < zero) EXIT |
---|
| 1459 | END DO |
---|
| 1460 | r = -r/p(j) |
---|
| 1461 | |
---|
| 1462 | DO |
---|
| 1463 | th = theta(j-1) + r*(theta(j) - theta(j-1)) |
---|
| 1464 | lambda = k - j + one - k*COS(th) |
---|
| 1465 | n = 1 |
---|
| 1466 | rlast = lambda |
---|
| 1467 | |
---|
| 1468 | DO |
---|
| 1469 | CALL RANDOM_NUMBER(r) |
---|
| 1470 | IF (r > rlast) EXIT |
---|
| 1471 | n = n + 1 |
---|
| 1472 | rlast = r |
---|
| 1473 | END DO |
---|
| 1474 | |
---|
| 1475 | IF (n .NE. 2*(n/2)) EXIT ! is n even? |
---|
| 1476 | CALL RANDOM_NUMBER(r) |
---|
| 1477 | END DO |
---|
| 1478 | |
---|
| 1479 | fn_val = SIGN(th, (r - rlast)/(one - rlast) - half) |
---|
| 1480 | RETURN |
---|
| 1481 | END FUNCTION random_von_Mises |
---|
| 1482 | |
---|
| 1483 | |
---|
| 1484 | |
---|
| 1485 | SUBROUTINE integral(a, b, result, dk) |
---|
| 1486 | |
---|
| 1487 | ! Gaussian integration of exp(k.cosx) from a to b. |
---|
| 1488 | |
---|
| 1489 | REAL (dp), INTENT(IN) :: dk |
---|
| 1490 | REAL, INTENT(IN) :: a, b |
---|
| 1491 | REAL, INTENT(OUT) :: result |
---|
| 1492 | |
---|
| 1493 | ! Local variables |
---|
| 1494 | |
---|
| 1495 | REAL (dp) :: xmid, range, x1, x2, & |
---|
| 1496 | x(3) = (/0.238619186083197_dp, 0.661209386466265_dp, 0.932469514203152_dp/), & |
---|
| 1497 | w(3) = (/0.467913934572691_dp, 0.360761573048139_dp, 0.171324492379170_dp/) |
---|
| 1498 | INTEGER :: i |
---|
| 1499 | |
---|
| 1500 | xmid = (a + b)/2._dp |
---|
| 1501 | range = (b - a)/2._dp |
---|
| 1502 | |
---|
| 1503 | result = 0._dp |
---|
| 1504 | DO i = 1, 3 |
---|
| 1505 | x1 = xmid + x(i)*range |
---|
| 1506 | x2 = xmid - x(i)*range |
---|
| 1507 | result = result + w(i)*(EXP(dk*COS(x1)) + EXP(dk*COS(x2))) |
---|
| 1508 | END DO |
---|
| 1509 | |
---|
| 1510 | result = result * range |
---|
| 1511 | RETURN |
---|
| 1512 | END SUBROUTINE integral |
---|
| 1513 | |
---|
| 1514 | |
---|
| 1515 | |
---|
| 1516 | FUNCTION random_Cauchy() RESULT(fn_val) |
---|
| 1517 | |
---|
| 1518 | ! Generate a random deviate from the standard Cauchy distribution |
---|
| 1519 | |
---|
| 1520 | REAL :: fn_val |
---|
| 1521 | |
---|
| 1522 | ! Local variables |
---|
| 1523 | REAL :: v(2) |
---|
| 1524 | |
---|
| 1525 | DO |
---|
| 1526 | CALL RANDOM_NUMBER(v) |
---|
| 1527 | v = two*(v - half) |
---|
| 1528 | IF (ABS(v(2)) < vsmall) CYCLE ! Test for zero |
---|
| 1529 | IF (v(1)**2 + v(2)**2 < one) EXIT |
---|
| 1530 | END DO |
---|
| 1531 | fn_val = v(1) / v(2) |
---|
| 1532 | |
---|
| 1533 | RETURN |
---|
| 1534 | END FUNCTION random_Cauchy |
---|
| 1535 | |
---|
| 1536 | |
---|
| 1537 | |
---|
| 1538 | SUBROUTINE random_order(order, n) |
---|
| 1539 | |
---|
| 1540 | ! Generate a random ordering of the integers 1 ... n. |
---|
| 1541 | |
---|
| 1542 | INTEGER, INTENT(IN) :: n |
---|
| 1543 | INTEGER, INTENT(OUT) :: order(n) |
---|
| 1544 | |
---|
| 1545 | ! Local variables |
---|
| 1546 | |
---|
| 1547 | INTEGER :: i, j, k |
---|
| 1548 | REAL :: wk |
---|
| 1549 | |
---|
| 1550 | DO i = 1, n |
---|
| 1551 | order(i) = i |
---|
| 1552 | END DO |
---|
| 1553 | |
---|
| 1554 | ! Starting at the end, swap the current last indicator with one |
---|
| 1555 | ! randomly chosen from those preceeding it. |
---|
| 1556 | |
---|
| 1557 | DO i = n, 2, -1 |
---|
| 1558 | CALL RANDOM_NUMBER(wk) |
---|
| 1559 | j = 1 + i * wk |
---|
| 1560 | IF (j < i) THEN |
---|
| 1561 | k = order(i) |
---|
| 1562 | order(i) = order(j) |
---|
| 1563 | order(j) = k |
---|
| 1564 | END IF |
---|
| 1565 | END DO |
---|
| 1566 | |
---|
| 1567 | RETURN |
---|
| 1568 | END SUBROUTINE random_order |
---|
| 1569 | |
---|
| 1570 | |
---|
| 1571 | |
---|
| 1572 | SUBROUTINE seed_random_number(iounit) |
---|
| 1573 | |
---|
| 1574 | INTEGER, INTENT(IN) :: iounit |
---|
| 1575 | |
---|
| 1576 | ! Local variables |
---|
| 1577 | |
---|
| 1578 | INTEGER :: k |
---|
| 1579 | INTEGER, ALLOCATABLE :: seed(:) |
---|
| 1580 | |
---|
| 1581 | CALL RANDOM_SEED(SIZE=k) |
---|
| 1582 | ALLOCATE( seed(k) ) |
---|
| 1583 | |
---|
| 1584 | WRITE(*, '(a, i2, a)')' Enter ', k, ' integers for random no. seeds: ' |
---|
| 1585 | READ(*, *) seed |
---|
| 1586 | WRITE(iounit, '(a, (7i10))') ' Random no. seeds: ', seed |
---|
| 1587 | CALL RANDOM_SEED(PUT=seed) |
---|
| 1588 | |
---|
| 1589 | DEALLOCATE( seed ) |
---|
| 1590 | |
---|
| 1591 | RETURN |
---|
| 1592 | END SUBROUTINE seed_random_number |
---|
| 1593 | |
---|
| 1594 | |
---|
| 1595 | END MODULE KPP_ROOT_Random |
---|
| 1596 | |
---|
| 1597 | MODULE KPP_ROOT_Integrator |
---|
| 1598 | USE KPP_ROOT_Random |
---|
| 1599 | USE KPP_ROOT_Parameters, ONLY : NVAR, NFIX, NREACT |
---|
| 1600 | USE KPP_ROOT_Global, ONLY : TIME, RCONST, Volume |
---|
| 1601 | USE KPP_ROOT_Stoichiom |
---|
| 1602 | USE KPP_ROOT_Stochastic |
---|
| 1603 | USE KPP_ROOT_Rates |
---|
| 1604 | USE KPP_ROOT_Random, ddp => dp |
---|
| 1605 | IMPLICIT NONE |
---|
| 1606 | |
---|
| 1607 | CONTAINS |
---|
| 1608 | |
---|
| 1609 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 1610 | SUBROUTINE TauLeap(Nsteps, Tau, T, SCT, NmlcV, NmlcF) |
---|
| 1611 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 1612 | ! |
---|
| 1613 | ! Tau-leap stochastic integration |
---|
| 1614 | ! INPUT: |
---|
| 1615 | ! Nsteps = no. of tau-leap steps to be simulated |
---|
| 1616 | ! Tau = time step length |
---|
| 1617 | ! T = time |
---|
| 1618 | ! SCT = stochastic rate constants |
---|
| 1619 | ! NmlcV, NmlcF = no. of molecules for variable and fixed species |
---|
| 1620 | ! OUTPUT: |
---|
| 1621 | ! T = updated time (after Nsteps) |
---|
| 1622 | ! NmlcV = updated no. of molecules for variable species |
---|
| 1623 | ! |
---|
| 1624 | ! |
---|
| 1625 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 1626 | IMPLICIT NONE |
---|
| 1627 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 1628 | |
---|
| 1629 | KPP_REAL:: T, Tau |
---|
| 1630 | INTEGER :: Nsteps |
---|
| 1631 | INTEGER :: NmlcV(NVAR), NmlcF(NFIX) |
---|
| 1632 | INTEGER :: i, j, irct, id, istep, Nfirings(NREACT) |
---|
| 1633 | REAL :: mu |
---|
| 1634 | KPP_REAL :: A(NREACT), SCT(NREACT), x |
---|
| 1635 | LOGICAL, SAVE :: First = .TRUE. |
---|
| 1636 | |
---|
| 1637 | DO istep = 1, Nsteps |
---|
| 1638 | |
---|
| 1639 | ! Propensity vector |
---|
| 1640 | CALL Propensity ( NmlcV, NmlcF, SCT, A ) |
---|
| 1641 | |
---|
| 1642 | ! Index of next reaction |
---|
| 1643 | DO irct = 1, NREACT |
---|
| 1644 | mu = A(irct)*Tau |
---|
| 1645 | Nfirings(irct) = random_Poisson(mu, First) |
---|
| 1646 | First = .TRUE. |
---|
| 1647 | END DO |
---|
| 1648 | |
---|
| 1649 | ! Update time with the leap interval |
---|
| 1650 | T = T + Tau; |
---|
| 1651 | |
---|
| 1652 | ! Directly update state vector |
---|
| 1653 | DO irct = 1, NREACT |
---|
| 1654 | DO i = CCOL_STOICM(irct), CCOL_STOICM(irct+1)-1 |
---|
| 1655 | id = IROW_STOICM(i) |
---|
| 1656 | NmlcV(id) = MAX(0, NmlcV(id) + Nfirings(irct)*INT(STOICM(i))) |
---|
| 1657 | END DO |
---|
| 1658 | END DO |
---|
| 1659 | |
---|
| 1660 | ! Update state vector |
---|
| 1661 | ! DO irct = 1, NREACT |
---|
| 1662 | ! DO j = 1, Nfirings(irct) |
---|
| 1663 | ! CALL MoleculeChange( irct, NmlcV ) |
---|
| 1664 | ! END DO |
---|
| 1665 | ! END DO |
---|
| 1666 | |
---|
| 1667 | END DO |
---|
| 1668 | |
---|
| 1669 | CONTAINS |
---|
| 1670 | |
---|
| 1671 | SUBROUTINE PropensityTemplate( T, NmlcV, NmlcF, Prop ) |
---|
| 1672 | KPP_REAL, INTENT(IN) :: T |
---|
| 1673 | INTEGER, INTENT(IN) :: NmlcV(NVAR), NmlcF(NFIX) |
---|
| 1674 | KPP_REAL, INTENT(OUT) :: Prop(NREACT) |
---|
| 1675 | KPP_REAL :: Tsave |
---|
| 1676 | ! Update the stochastic reaction rates, which may be time dependent |
---|
| 1677 | Tsave = TIME |
---|
| 1678 | TIME = T |
---|
| 1679 | CALL Update_RCONST() |
---|
| 1680 | CALL StochasticRates( RCONST, Volume, SCT ) |
---|
| 1681 | CALL Propensity ( NmlcV, NmlcF, SCT, Prop ) |
---|
| 1682 | TIME = Tsave |
---|
| 1683 | END SUBROUTINE PropensityTemplate |
---|
| 1684 | |
---|
| 1685 | END SUBROUTINE TauLeap |
---|
| 1686 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 1687 | |
---|
| 1688 | END MODULE KPP_ROOT_Integrator |
---|