source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/int/tau_leap.f90 @ 3997

Last change on this file since 3997 was 2696, checked in by kanani, 6 years ago

Merge of branch palm4u into trunk

File size: 41.1 KB
Line 
1MODULE 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
96IMPLICIT NONE
97REAL, PRIVATE      :: zero = 0.0, half = 0.5, one = 1.0, two = 2.0,   &
98                      vsmall = TINY(1.0), vlarge = HUGE(1.0)
99PRIVATE            :: integral
100INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60)
101
102
103CONTAINS
104
105
106FUNCTION 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
119REAL :: fn_val
120
121!     Local variables
122REAL     :: 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
127DO
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
143END DO
144
145!     Return ratio of P's coordinates as the normal deviate
146fn_val = v/u
147RETURN
148
149END FUNCTION random_normal
150
151
152
153FUNCTION 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
166REAL, INTENT(IN)    :: s
167LOGICAL, INTENT(IN) :: first
168REAL                :: fn_val
169
170IF (s <= zero) THEN
171  WRITE(*, *) 'SHAPE PARAMETER VALUE MUST BE POSITIVE'
172  STOP
173END IF
174
175IF (s > one) THEN
176  fn_val = random_gamma1(s, first)
177ELSE IF (s < one) THEN
178  fn_val = random_gamma2(s, first)
179ELSE
180  fn_val = random_exponential()
181END IF
182
183RETURN
184END FUNCTION random_gamma
185
186
187
188FUNCTION 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
196REAL, INTENT(IN)    :: s
197LOGICAL, INTENT(IN) :: first
198REAL                :: fn_val
199
200! Local variables
201REAL, SAVE  :: c, d
202REAL        :: u, v, x
203
204IF (first) THEN
205  d = s - one/3.
206  c = one/SQRT(9.0*d)
207END IF
208
209! Start of main loop
210DO
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
230END DO
231
232RETURN
233END FUNCTION random_gamma1
234
235
236
237FUNCTION 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
251REAL, INTENT(IN)    :: s
252LOGICAL, INTENT(IN) :: first
253REAL                :: fn_val
254
255!     Local variables
256REAL       :: r, x, w
257REAL, SAVE :: a, p, c, uf, vr, d
258
259IF (s <= zero .OR. s >= one) THEN
260  WRITE(*, *) 'SHAPE PARAMETER VALUE OUTSIDE PERMITTED RANGE'
261  STOP
262END IF
263
264IF (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)
275END IF
276
277DO
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
298END DO
299
300fn_val = x
301RETURN
302
303END FUNCTION random_gamma2
304
305
306
307FUNCTION 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
312INTEGER, INTENT(IN) :: ndf
313LOGICAL, INTENT(IN) :: first
314REAL                :: fn_val
315
316fn_val = two * random_gamma(half*ndf, first)
317RETURN
318
319END FUNCTION random_chisq
320
321
322
323FUNCTION 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
333REAL  :: fn_val
334
335!     Local variable
336REAL  :: r
337
338DO
339  CALL RANDOM_NUMBER(r)
340  IF (r > zero) EXIT
341END DO
342
343fn_val = -LOG(r)
344RETURN
345
346END FUNCTION random_exponential
347
348
349
350FUNCTION 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
358REAL, INTENT(IN) :: a
359REAL             :: fn_val
360
361!     For speed, there is no checking that a is not zero or very small.
362
363fn_val = random_exponential() ** (one/a)
364RETURN
365
366END FUNCTION random_Weibull
367
368
369
370FUNCTION 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
384REAL, INTENT(IN)    :: aa, bb
385LOGICAL, INTENT(IN) :: first
386REAL                :: fn_val
387
388!     Local variables
389REAL, PARAMETER  :: aln4 = 1.3862944
390REAL             :: a, b, g, r, s, x, y, z
391REAL, SAVE       :: d, f, h, t, c
392LOGICAL, SAVE    :: swap
393
394IF (aa <= zero .OR. bb <= zero) THEN
395  WRITE(*, *) 'IMPERMISSIBLE SHAPE PARAMETER VALUE(S)'
396  STOP
397END IF
398
399IF (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
418END IF
419
420DO
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
439END DO
440
441IF (swap) fn_val = one - fn_val
442RETURN
443END FUNCTION random_beta
444
445
446
447FUNCTION 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
459INTEGER, INTENT(IN) :: m
460REAL                :: fn_val
461
462!     Local variables
463REAL, SAVE      :: s, c, a, f, g
464REAL            :: r, x, v
465
466REAL, PARAMETER :: three = 3.0, four = 4.0, quart = 0.25,   &
467                   five = 5.0, sixteen = 16.0
468INTEGER         :: mm = 0
469
470IF (m < 1) THEN
471  WRITE(*, *) 'IMPERMISSIBLE DEGREES OF FREEDOM'
472  STOP
473END IF
474
475IF (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
487END IF
488
489DO
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
500END DO
501
502fn_val = x
503RETURN
504END FUNCTION random_t
505
506
507
508SUBROUTINE 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
541INTEGER, INTENT(IN)   :: n
542REAL, INTENT(IN)      :: h(:), d(:)   ! d(n*(n+1)/2)
543REAL, INTENT(IN OUT)  :: f(:)         ! f(n*(n+1)/2)
544REAL, INTENT(OUT)     :: x(:)
545LOGICAL, INTENT(IN)   :: first
546INTEGER, INTENT(OUT)  :: ier
547
548!     Local variables
549INTEGER       :: j, i, m
550REAL          :: y, v
551INTEGER, SAVE :: n2
552
553IF (n < 1) THEN
554  WRITE(*, *) 'SIZE OF VECTOR IS NON POSITIVE'
555  STOP
556END IF
557
558ier = 0
559IF (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
594END IF
595
596x(1:n) = h(1:n)
597DO 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
602END DO ! j = 1,n
603
604RETURN
605END SUBROUTINE random_mvnorm
606
607
608
609FUNCTION 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
623REAL, INTENT(IN)    :: h, b
624LOGICAL, INTENT(IN) :: first
625REAL                :: fn_val
626
627!     Local variables
628REAL            :: ym, xm, r, w, r1, r2, x
629REAL, SAVE      :: a, c, d, e
630REAL, PARAMETER :: quart = 0.25
631
632IF (h < zero .OR. b <= zero) THEN
633  WRITE(*, *) 'IMPERMISSIBLE DISTRIBUTION PARAMETER VALUES'
634  STOP
635END IF
636
637IF (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
662END IF
663
664DO
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
671END DO
672
673fn_val = x
674
675RETURN
676END FUNCTION random_inv_gauss
677
678
679
680FUNCTION 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 ..
727REAL, INTENT(IN)    :: mu
728LOGICAL, INTENT(IN) :: first
729INTEGER             :: ival
730!     ..
731!     .. Local Scalars ..
732REAL          :: b1, b2, c, c0, c1, c2, c3, del, difmuk, e, fk, fx, fy, g,  &
733                 omega, px, py, t, u, v, x, xx
734REAL, SAVE    :: s, d, p, q, p0
735INTEGER       :: j, k, kflag
736LOGICAL, SAVE :: full_init
737INTEGER, SAVE :: l, m
738!     ..
739!     .. Local Arrays ..
740REAL, SAVE    :: pp(35)
741!     ..
742!     .. Data statements ..
743REAL, PARAMETER :: a0 = -.5, a1 = .3333333, a2 = -.2500068, a3 = .2000118,  &
744                   a4 = -.1661269, a5 = .1421878, a6 = -.1384794,   &
745                   a7 = .1250060
746
747REAL, PARAMETER :: fact(10) = (/ 1., 1., 2., 6., 24., 120., 720., 5040.,  &
748                                 40320., 362880. /)
749
750!     ..
751!     .. Executable Statements ..
752IF (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
871ELSE
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
915END IF
916
917RETURN
918END FUNCTION random_Poisson
919
920
921
922FUNCTION 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
938INTEGER, INTENT(IN) :: n
939REAL, INTENT(IN)    :: p
940LOGICAL, INTENT(IN) :: first
941INTEGER             :: ival
942
943!     Local variables
944
945INTEGER         :: ru, rd
946INTEGER, SAVE   :: r0
947REAL            :: u, pd, pu
948REAL, SAVE      :: odds_ratio, p_r
949REAL, PARAMETER :: zero = 0.0, one = 1.0
950
951IF (first) THEN
952  r0 = (n+1)*p
953  p_r = bin_prob(n, p, r0)
954  odds_ratio = p / (one - p)
955END IF
956
957CALL RANDOM_NUMBER(u)
958u = u - p_r
959IF (u < zero) THEN
960  ival = r0
961  RETURN
962END IF
963
964pu = p_r
965ru = r0
966pd = p_r
967rd = r0
968DO
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
988END DO
989
990!     This point should not be reached, but just in case:
991
992ival = r0
993RETURN
994
995END FUNCTION random_binomial1
996
997
998
999FUNCTION bin_prob(n, p, r) RESULT(fn_val)
1000!     Calculate a binomial probability
1001
1002INTEGER, INTENT(IN) :: n, r
1003REAL, INTENT(IN)    :: p
1004REAL                :: fn_val
1005
1006!     Local variable
1007REAL                :: one = 1.0
1008
1009fn_val = EXP( lngamma(DBLE(n+1)) - lngamma(DBLE(r+1)) - lngamma(DBLE(n-r+1)) &
1010              + r*LOG(p) + (n-r)*LOG(one - p) )
1011RETURN
1012
1013END FUNCTION bin_prob
1014
1015
1016
1017FUNCTION 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
1026REAL (dp), INTENT(IN) :: x
1027REAL (dp)             :: fn_val
1028
1029!       Local variables
1030
1031REAL (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
1035LOGICAL   :: reflect
1036
1037!       lngamma is not defined if x = 0 or a negative integer.
1038
1039IF (x > 0.d0) GO TO 10
1040IF (x /= INT(x)) GO TO 10
1041fn_val = 0.d0
1042RETURN
1043
1044!       If x < 0, use the reflection formula:
1045!               gamma(x) * gamma(1-x) = pi * cosec(pi.x)
1046
104710 reflect = (x < 0.d0)
1048IF (reflect) THEN
1049  arg = 1.d0 - x
1050ELSE
1051  arg = x
1052END IF
1053
1054!       Increase the argument, if necessary, to make it > 10.
1055
1056product = 1.d0
105720 IF (arg <= 10.d0) THEN
1058  product = product * arg
1059  arg = arg + 1.d0
1060  GO TO 20
1061END 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
1068arg = arg - 0.5D0
1069temp = 1.d0/arg**2
1070fn_val = lnrt2pi + arg * (LOG(arg) - 1.d0 + &
1071                  (((a4*temp + a3)*temp + a2)*temp + a1)*temp) - LOG(product)
1072IF (reflect) THEN
1073  temp = SIN(pi * x)
1074  fn_val = LOG(pi/temp) - fn_val
1075END IF
1076RETURN
1077END FUNCTION lngamma
1078
1079
1080
1081FUNCTION 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 ..
1143REAL, INTENT(IN)    :: pp
1144INTEGER, INTENT(IN) :: n
1145LOGICAL, INTENT(IN) :: first
1146INTEGER             :: ival
1147!     ..
1148!     .. Local Scalars ..
1149REAL            :: alv, amaxp, f, f1, f2, u, v, w, w2, x, x1, x2, ynorm, z, z2
1150REAL, PARAMETER :: zero = 0.0, half = 0.5, one = 1.0
1151INTEGER         :: i, ix, ix1, k, mp
1152INTEGER, SAVE   :: m
1153REAL, 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
1161IF (first) THEN
1162  p = MIN(pp, one-pp)
1163  q = one - p
1164  xnp = n * p
1165END IF
1166
1167IF (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
1284ELSE
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
1302END IF
1303
1304110 IF (pp > half) ix = n - ix
1305ival = ix
1306RETURN
1307
1308END FUNCTION random_binomial2
1309
1310
1311
1312
1313FUNCTION 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
1332REAL, INTENT(IN)   :: sk, p
1333INTEGER            :: ival
1334
1335!     Local variables
1336! THE PARAMETER ULN = -LOG(MACHINE'S SMALLEST REAL NUMBER).
1337
1338REAL, PARAMETER    :: h = 0.7
1339REAL               :: q, x, st, uln, v, r, s, y, g
1340INTEGER            :: k, i, n
1341
1342IF (sk <= zero .OR. p <= zero .OR. p >= one) THEN
1343  WRITE(*, *) 'IMPERMISSIBLE DISTRIBUTION PARAMETER VALUES'
1344  STOP
1345END IF
1346
1347q = one - p
1348x = zero
1349st = sk
1350IF (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
1362END IF
1363
1364s = zero
1365uln = -LOG(vsmall)
1366IF (st > -uln/LOG(q)) THEN
1367  WRITE(*, *) ' P IS TOO LARGE FOR THIS VALUE OF SK'
1368  STOP
1369END IF
1370
1371y = q**st
1372g = st
1373CALL RANDOM_NUMBER(r)
1374DO
1375  IF (y > r) EXIT
1376  r = r - y
1377  s = s + one
1378  y = y*p*g/s
1379  g = g + one
1380END DO
1381
1382ival = x + s + half
1383RETURN
1384END FUNCTION random_neg_binomial
1385
1386
1387
1388FUNCTION 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
1404REAL, INTENT(IN)     :: k
1405LOGICAL, INTENT(IN)  :: first
1406REAL                 :: fn_val
1407
1408!     Local variables
1409
1410INTEGER          :: j, n
1411INTEGER, SAVE    :: nk
1412REAL, PARAMETER  :: pi = 3.14159265
1413REAL, SAVE       :: p(20), theta(0:20)
1414REAL             :: sump, r, th, lambda, rlast
1415REAL (dp)        :: dk
1416
1417IF (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
1453END IF                           ! if first
1454
1455CALL RANDOM_NUMBER(r)
1456DO j = 1, nk
1457  r = r - p(j)
1458  IF (r < zero) EXIT
1459END DO
1460r = -r/p(j)
1461
1462DO
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)
1477END DO
1478
1479fn_val = SIGN(th, (r - rlast)/(one - rlast) - half)
1480RETURN
1481END FUNCTION random_von_Mises
1482
1483
1484
1485SUBROUTINE integral(a, b, result, dk)
1486
1487!     Gaussian integration of exp(k.cosx) from a to b.
1488
1489REAL (dp), INTENT(IN) :: dk
1490REAL, INTENT(IN)      :: a, b
1491REAL, INTENT(OUT)     :: result
1492
1493!     Local variables
1494
1495REAL (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/)
1498INTEGER    :: i
1499
1500xmid = (a + b)/2._dp
1501range = (b - a)/2._dp
1502
1503result = 0._dp
1504DO 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)))
1508END DO
1509
1510result = result * range
1511RETURN
1512END SUBROUTINE integral
1513
1514
1515
1516FUNCTION random_Cauchy() RESULT(fn_val)
1517
1518!     Generate a random deviate from the standard Cauchy distribution
1519
1520REAL     :: fn_val
1521
1522!     Local variables
1523REAL     :: v(2)
1524
1525DO
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
1530END DO
1531fn_val = v(1) / v(2)
1532
1533RETURN
1534END FUNCTION random_Cauchy
1535
1536
1537
1538SUBROUTINE random_order(order, n)
1539
1540!     Generate a random ordering of the integers 1 ... n.
1541
1542INTEGER, INTENT(IN)  :: n
1543INTEGER, INTENT(OUT) :: order(n)
1544
1545!     Local variables
1546
1547INTEGER :: i, j, k
1548REAL    :: wk
1549
1550DO i = 1, n
1551  order(i) = i
1552END DO
1553
1554!     Starting at the end, swap the current last indicator with one
1555!     randomly chosen from those preceeding it.
1556
1557DO 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
1565END DO
1566
1567RETURN
1568END SUBROUTINE random_order
1569
1570
1571
1572SUBROUTINE seed_random_number(iounit)
1573
1574INTEGER, INTENT(IN)  :: iounit
1575
1576! Local variables
1577
1578INTEGER              :: k
1579INTEGER, ALLOCATABLE :: seed(:)
1580
1581CALL RANDOM_SEED(SIZE=k)
1582ALLOCATE( seed(k) )
1583
1584WRITE(*, '(a, i2, a)')' Enter ', k, ' integers for random no. seeds: '
1585READ(*, *) seed
1586WRITE(iounit, '(a, (7i10))') ' Random no. seeds: ', seed
1587CALL RANDOM_SEED(PUT=seed)
1588
1589DEALLOCATE( seed )
1590
1591RETURN
1592END SUBROUTINE seed_random_number
1593
1594
1595END MODULE KPP_ROOT_Random
1596
1597MODULE 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
1607CONTAINS
1608
1609!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1610SUBROUTINE 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     
1669CONTAINS
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   
1685END SUBROUTINE TauLeap
1686!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1687
1688END MODULE KPP_ROOT_Integrator
Note: See TracBrowser for help on using the repository browser.