 # module for random number generation

## Found at: ftp.icm.edu.pl:70/packages/netlib/random/random.f90

```MODULE random
```
```! A module for random number generation from the following distributions:
```
```!
```
```!     Distribution                    Function/subroutine name
```
```!
```
```!     Normal (Gaussian)               random_normal
```
```!     Gamma                           random_gamma
```
```!     Chi-squared                     random_chisq
```
```!     Exponential                     random_exponential
```
```!     Weibull                         random_Weibull
```
```!     Beta                            random_beta
```
```!     t                               random_t
```
```!     Multivariate normal             random_mvnorm
```
```!     Generalized inverse Gaussian    random_inv_gauss
```
```!     Poisson                         random_Poisson
```
```!     Binomial                        random_binomial1   *
```
```!                                     random_binomial2   *
```
```!     Negative binomial               random_neg_binomial
```
```!     von Mises                       random_von_Mises
```
```!     Cauchy                          random_Cauchy
```
```!
```
```!  Generate a random ordering of the integers 1 .. N
```
```!                                     random_order
```
```!     Initialize (seed) the uniform random number generator for ANY compiler
```
```!                                     seed_random_number
```

```!     Lognormal - see note below.
```

```!  ** Two functions are provided for the binomial distribution.
```
```!  If the parameter values remain constant, it is recommended that the
```
```!  first function is used (random_binomial1).   If one or both of the
```
```!  parameters change, use the second function (random_binomial2).
```

```! The compilers own random number generator, SUBROUTINE RANDOM_NUMBER(r),
```
```! is used to provide a source of uniformly distributed random numbers.
```

```! N.B. At this stage, only one random number is generated at each call to
```
```!      one of the functions above.
```

```! The module uses the following functions which are included here:
```
```! bin_prob to calculate a single binomial probability
```
```! lngamma  to calculate the logarithm to base e of the gamma function
```

```! Some of the code is adapted from Dagpunar's book:
```
```!     Dagpunar, J. 'Principles of random variate generation'
```
```!     Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9
```
```!
```
```! In most of Dagpunar's routines, there is a test to see whether the value
```
```! of one or two floating-point parameters has changed since the last call.
```
```! These tests have been replaced by using a logical variable FIRST.
```
```! This should be set to .TRUE. on the first call using new values of the
```
```! parameters, and .FALSE. if the parameter values are the same as for the
```
```! previous call.
```

```!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
```
```! Lognormal distribution
```
```! If X has a lognormal distribution, then log(X) is normally distributed.
```
```! Here the logarithm is the natural logarithm, that is to base e, sometimes
```
```! denoted as ln.  To generate random variates from this distribution, generate
```
```! a random deviate from the normal distribution with mean and variance equal
```
```! to the mean and variance of the logarithms of X, then take its exponential.
```

```! Relationship between the mean & variance of log(X) and the mean & variance
```
```! of X, when X has a lognormal distribution.
```
```! Let m = mean of log(X), and s^2 = variance of log(X)
```
```! Then
```
```! mean of X     = exp(m + 0.5s^2)
```
```! variance of X = (mean(X))^2.[exp(s^2) - 1]
```

```! In the reverse direction (rarely used)
```
```! variance of log(X) = log[1 + var(X)/(mean(X))^2]
```
```! mean of log(X)     = log(mean(X) - 0.5var(log(X))
```

```! N.B. The above formulae relate to population parameters; they will only be
```
```!      approximate if applied to sample values.
```
```!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
```

```! Version 1.13, 2 October 2000
```
```! Changes from version 1.01
```
```! 1. The random_order, random_Poisson & random_binomial routines have been
```
```!    replaced with more efficient routines.
```
```! 2. A routine, seed_random_number, has been added to seed the uniform random
```
```!    number generator.   This requires input of the required number of seeds
```
```!    for the particular compiler from a specified I/O unit such as a keyboard.
```
```! 3. Made compatible with Lahey's ELF90.
```
```! 4. Marsaglia & Tsang algorithm used for random_gamma when shape parameter > 1.
```
```! 5. INTENT for array f corrected in random_mvnorm.
```

```!     Author: Alan Miller
```
```!             CSIRO Division of Mathematical & Information Sciences
```
```!             Private Bag 10, Clayton South MDC
```
```!             Clayton 3169, Victoria, Australia
```
```!     Phone: (+61) 3 9545-8016      Fax: (+61) 3 9545-8080
```
```!     e-mail: amiller @ bigpond.net.au
```

```REAL, PRIVATE      :: zero = 0.0, half = 0.5, one = 1.0, two = 2.0,   &
```
```                      vsmall = TINY(1.0), vlarge = HUGE(1.0)
```

```CONTAINS
```

```FUNCTION random_normal() RESULT(fn_val)
```

```! Adapted from the following Fortran 77 code
```
```!      ALGORITHM 712, COLLECTED ALGORITHMS FROM ACM.
```
```!      THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
```
```!      VOL. 18, NO. 4, DECEMBER, 1992, PP. 434-435.
```

```!  The function random_normal() returns a normally distributed pseudo-random
```
```!  number with zero mean and unit variance.
```

```!  The algorithm uses the ratio of uniforms method of A.J. Kinderman
```
```!  and J.F. Monahan augmented with quadratic bounding curves.
```

```REAL :: fn_val
```

```!     Local variables
```
```REAL     :: s = 0.449871, t = -0.386595, a = 0.19600, b = 0.25472,           &
```
```            r1 = 0.27597, r2 = 0.27846, u, v, x, y, q
```

```!     Generate P = (u,v) uniform in rectangle enclosing acceptance region
```

```DO
```
```  CALL RANDOM_NUMBER(u)
```
```  CALL RANDOM_NUMBER(v)
```
```  v = 1.7156 * (v - half)
```

```!     Evaluate the quadratic form
```
```  x = u - s
```
```  y = ABS(v) - t
```
```  q = x**2 + y*(a*y - b*x)
```

```!     Accept P if inside inner ellipse
```
```  IF (q < r1) EXIT
```
```!     Reject P if outside outer ellipse
```
```  IF (q > r2) CYCLE
```
```!     Reject P if outside acceptance region
```
```  IF (v**2 < -4.0*LOG(u)*u**2) EXIT
```
```END DO
```

```!     Return ratio of P's coordinates as the normal deviate
```
```fn_val = v/u
```
```RETURN
```

```END FUNCTION random_normal
```

```FUNCTION random_gamma(s, first) RESULT(fn_val)
```

```! Adapted from Fortran 77 code from the book:
```
```!     Dagpunar, J. 'Principles of random variate generation'
```
```!     Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9
```

```!     FUNCTION GENERATES A RANDOM GAMMA VARIATE.
```
```!     CALLS EITHER random_gamma1 (S > 1.0)
```
```!     OR random_exponential (S = 1.0)
```
```!     OR random_gamma2 (S < 1.0).
```

```!     S = SHAPE PARAMETER OF DISTRIBUTION (0 < REAL).
```

```REAL, INTENT(IN)    :: s
```
```LOGICAL, INTENT(IN) :: first
```
```REAL                :: fn_val
```

```  WRITE(*, *) 'SHAPE PARAMETER VALUE MUST BE POSITIVE'
```
```  STOP
```
```END IF
```

```  fn_val = random_gamma1(s, first)
```
```ELSE IF (s < one) THEN
```
```  fn_val = random_gamma2(s, first)
```
```ELSE
```
```  fn_val = random_exponential()
```
```END IF
```

```RETURN
```
```END FUNCTION random_gamma
```

```FUNCTION random_gamma1(s, first) RESULT(fn_val)
```

```! Uses the algorithm in
```
```! Marsaglia, G. and Tsang, W.W. (2000) `A simple method for generating
```
```! gamma variables', Trans. om Math. Software (TOMS), vol.26(3), pp.363-372.
```

```! Generates a random gamma deviate for shape parameter s >= 1.
```

```REAL, INTENT(IN)    :: s
```
```LOGICAL, INTENT(IN) :: first
```
```REAL                :: fn_val
```

```! Local variables
```
```REAL, SAVE  :: c, d
```
```REAL        :: u, v, x
```

```  d = s - one/3.
```
```  c = one/SQRT(9.0*d)
```
```END IF
```

```! Start of main loop
```
```DO
```

```! Generate v = (1+cx)^3 where x is random normal; repeat if v <= 0.
```

```  DO
```
```    x = random_normal()
```
```    v = (one + c*x)**3
```
```    IF (v > zero) EXIT
```
```  END DO
```

```! Generate uniform variable U
```

```  CALL RANDOM_NUMBER(u)
```
```  IF (u < one - 0.0331*x**4) THEN
```
```    fn_val = d*v
```
```    EXIT
```
```  ELSE IF (LOG(u) < half*x**2 + d*(one - v + LOG(v))) THEN
```
```    fn_val = d*v
```
```    EXIT
```
```  END IF
```
```END DO
```

```RETURN
```
```END FUNCTION random_gamma1
```

```FUNCTION random_gamma2(s, first) RESULT(fn_val)
```

```! Adapted from Fortran 77 code from the book:
```
```!     Dagpunar, J. 'Principles of random variate generation'
```
```!     Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9
```

```! FUNCTION GENERATES A RANDOM VARIATE IN [0,INFINITY) FROM
```
```! A GAMMA DISTRIBUTION WITH DENSITY PROPORTIONAL TO
```
```! GAMMA2**(S-1) * EXP(-GAMMA2),
```
```! USING A SWITCHING METHOD.
```

```!    S = SHAPE PARAMETER OF DISTRIBUTION
```
```!          (REAL < 1.0)
```

```REAL, INTENT(IN)    :: s
```
```LOGICAL, INTENT(IN) :: first
```
```REAL                :: fn_val
```

```!     Local variables
```
```REAL       :: r, x, w
```
```REAL, SAVE :: a, p, c, uf, vr, d
```

```  WRITE(*, *) 'SHAPE PARAMETER VALUE OUTSIDE PERMITTED RANGE'
```
```  STOP
```
```END IF
```

```  a = one - s
```
```  p = a/(a + s*EXP(-a))
```
```  IF (s < vsmall) THEN
```
```    WRITE(*, *) 'SHAPE PARAMETER VALUE TOO SMALL'
```
```    STOP
```
```  END IF
```
```  c = one/s
```
```  uf = p*(vsmall/a)**s
```
```  vr = one - vsmall
```
```  d = a*LOG(a)
```
```END IF
```

```DO
```
```  CALL RANDOM_NUMBER(r)
```
```  IF (r >= vr) THEN
```
```    CYCLE
```
```  ELSE IF (r > p) THEN
```
```    x = a - LOG((one - r)/(one - p))
```
```    w = a*LOG(x)-d
```
```  ELSE IF (r > uf) THEN
```
```    x = a*(r/p)**c
```
```    w = x
```
```  ELSE
```
```    fn_val = zero
```
```    RETURN
```
```  END IF
```

```  CALL RANDOM_NUMBER(r)
```
```  IF (one-r <= w .AND. r > zero) THEN
```
```    IF (r*(w + one) >= one) CYCLE
```
```    IF (-LOG(r) <= w) CYCLE
```
```  END IF
```
```  EXIT
```
```END DO
```

```fn_val = x
```
```RETURN
```

```END FUNCTION random_gamma2
```

```FUNCTION random_chisq(ndf, first) RESULT(fn_val)
```

```!     Generates a random variate from the chi-squared distribution with
```
```!     ndf degrees of freedom
```

```LOGICAL, INTENT(IN) :: first
```
```REAL                :: fn_val
```

```fn_val = two * random_gamma(half*ndf, first)
```
```RETURN
```

```END FUNCTION random_chisq
```

```FUNCTION random_exponential() RESULT(fn_val)
```

```! Adapted from Fortran 77 code from the book:
```
```!     Dagpunar, J. 'Principles of random variate generation'
```
```!     Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9
```

```! FUNCTION GENERATES A RANDOM VARIATE IN [0,INFINITY) FROM
```
```! A NEGATIVE EXPONENTIAL DlSTRIBUTION WlTH DENSITY PROPORTIONAL
```
```! TO EXP(-random_exponential), USING INVERSION.
```

```REAL  :: fn_val
```

```!     Local variable
```
```REAL  :: r
```

```DO
```
```  CALL RANDOM_NUMBER(r)
```
```  IF (r > zero) EXIT
```
```END DO
```

```fn_val = -LOG(r)
```
```RETURN
```

```END FUNCTION random_exponential
```

```FUNCTION random_Weibull(a) RESULT(fn_val)
```

```!     Generates a random variate from the Weibull distribution with
```
```!     probability density:
```
```!                      a
```
```!               a-1  -x
```
```!     f(x) = a.x    e
```

```REAL, INTENT(IN) :: a
```
```REAL             :: fn_val
```

```!     For speed, there is no checking that a is not zero or very small.
```

```fn_val = random_exponential() ** (one/a)
```
```RETURN
```

```END FUNCTION random_Weibull
```

```FUNCTION random_beta(aa, bb, first) RESULT(fn_val)
```

```! Adapted from Fortran 77 code from the book:
```
```!     Dagpunar, J. 'Principles of random variate generation'
```
```!     Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9
```

```! FUNCTION GENERATES A RANDOM VARIATE IN [0,1]
```
```! FROM A BETA DISTRIBUTION WITH DENSITY
```
```! PROPORTIONAL TO BETA**(AA-1) * (1-BETA)**(BB-1).
```
```! USING CHENG'S LOG LOGISTIC METHOD.
```

```!     AA = SHAPE PARAMETER FROM DISTRIBUTION (0 < REAL)
```
```!     BB = SHAPE PARAMETER FROM DISTRIBUTION (0 < REAL)
```

```REAL, INTENT(IN)    :: aa, bb
```
```LOGICAL, INTENT(IN) :: first
```
```REAL                :: fn_val
```

```!     Local variables
```
```REAL, PARAMETER  :: aln4 = 1.3862944
```
```REAL             :: a, b, g, r, s, x, y, z
```
```REAL, SAVE       :: d, f, h, t, c
```
```LOGICAL, SAVE    :: swap
```

```  WRITE(*, *) 'IMPERMISSIBLE SHAPE PARAMETER VALUE(S)'
```
```  STOP
```
```END IF
```

```  a = aa
```
```  b = bb
```
```  swap = b > a
```
```  IF (swap) THEN
```
```    g = b
```
```    b = a
```
```    a = g
```
```  END IF
```
```  d = a/b
```
```  f = a+b
```
```  IF (b > one) THEN
```
```    h = SQRT((two*a*b - f)/(f - two))
```
```    t = one
```
```  ELSE
```
```    h = b
```
```    t = one/(one + (a/(vlarge*b))**b)
```
```  END IF
```
```  c = a+h
```
```END IF
```

```DO
```
```  CALL RANDOM_NUMBER(r)
```
```  CALL RANDOM_NUMBER(x)
```
```  s = r*r*x
```
```  IF (r < vsmall .OR. s <= zero) CYCLE
```
```  IF (r < t) THEN
```
```    x = LOG(r/(one - r))/h
```
```    y = d*EXP(x)
```
```    z = c*x + f*LOG((one + d)/(one + y)) - aln4
```
```    IF (s - one > z) THEN
```
```      IF (s - s*z > one) CYCLE
```
```      IF (LOG(s) > z) CYCLE
```
```    END IF
```
```    fn_val = y/(one + y)
```
```  ELSE
```
```    IF (4.0*s > (one + one/d)**f) CYCLE
```
```    fn_val = one
```
```  END IF
```
```  EXIT
```
```END DO
```

```RETURN
```
```END FUNCTION random_beta
```

```FUNCTION random_t(m) RESULT(fn_val)
```

```! Adapted from Fortran 77 code from the book:
```
```!     Dagpunar, J. 'Principles of random variate generation'
```
```!     Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9
```

```! FUNCTION GENERATES A RANDOM VARIATE FROM A
```
```! T DISTRIBUTION USING KINDERMAN AND MONAHAN'S RATIO METHOD.
```

```!     M = DEGREES OF FREEDOM OF DISTRIBUTION
```
```!           (1 <= 1NTEGER)
```

```REAL                :: fn_val
```

```!     Local variables
```
```REAL, SAVE      :: s, c, a, f, g
```
```REAL            :: r, x, v
```

```REAL, PARAMETER :: three = 3.0, four = 4.0, quart = 0.25,   &
```
```                   five = 5.0, sixteen = 16.0
```

```  WRITE(*, *) 'IMPERMISSIBLE DEGREES OF FREEDOM'
```
```  STOP
```
```END IF
```

```  s = m
```
```  c = -quart*(s + one)
```
```  a = four/(one + one/s)**c
```
```  f = sixteen/a
```
```  IF (m > 1) THEN
```
```    g = s - one
```
```    g = ((s + one)/g)**c*SQRT((s+s)/g)
```
```  ELSE
```
```    g = one
```
```  END IF
```
```  mm = m
```
```END IF
```

```DO
```
```  CALL RANDOM_NUMBER(r)
```
```  IF (r <= zero) CYCLE
```
```  CALL RANDOM_NUMBER(v)
```
```  x = (two*v - one)*g/r
```
```  v = x*x
```
```  IF (v > five - a*r) THEN
```
```    IF (m >= 1 .AND. r*(v + three) > f) CYCLE
```
```    IF (r > (one + v/s)**c) CYCLE
```
```  END IF
```
```  EXIT
```
```END DO
```

```fn_val = x
```
```RETURN
```
```END FUNCTION random_t
```

```SUBROUTINE random_mvnorm(n, h, d, f, first, x, ier)
```

```! Adapted from Fortran 77 code from the book:
```
```!     Dagpunar, J. 'Principles of random variate generation'
```
```!     Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9
```

```! N.B. An extra argument, ier, has been added to Dagpunar's routine
```

```!     SUBROUTINE GENERATES AN N VARIATE RANDOM NORMAL
```
```!     VECTOR USING A CHOLESKY DECOMPOSITION.
```

```! ARGUMENTS:
```
```!        N = NUMBER OF VARIATES IN VECTOR
```
```!           (INPUT,INTEGER >= 1)
```
```!     H(J) = J'TH ELEMENT OF VECTOR OF MEANS
```
```!           (INPUT,REAL)
```
```!     X(J) = J'TH ELEMENT OF DELIVERED VECTOR
```
```!           (OUTPUT,REAL)
```
```!
```
```!    D(J*(J-1)/2+I) = (I,J)'TH ELEMENT OF VARIANCE MATRIX (J> = I)
```
```!            (INPUT,REAL)
```
```!    F((J-1)*(2*N-J)/2+I) = (I,J)'TH ELEMENT OF LOWER TRIANGULAR
```
```!           DECOMPOSITION OF VARIANCE MATRIX (J <= I)
```
```!            (OUTPUT,REAL)
```

```!    FIRST = .TRUE. IF THIS IS THE FIRST CALL OF THE ROUTINE
```
```!    OR IF THE DISTRIBUTION HAS CHANGED SINCE THE LAST CALL OF THE ROUTINE.
```
```!    OTHERWISE SET TO .FALSE.
```
```!            (INPUT,LOGICAL)
```

```!    ier = 1 if the input covariance matrix is not +ve definite
```
```!        = 0 otherwise
```

```REAL, INTENT(IN)      :: h(:), d(:)   ! d(n*(n+1)/2)
```
```REAL, INTENT(IN OUT)  :: f(:)         ! f(n*(n+1)/2)
```
```REAL, INTENT(OUT)     :: x(:)
```
```LOGICAL, INTENT(IN)   :: first
```

```!     Local variables
```
```REAL          :: y, v
```

```  WRITE(*, *) 'SIZE OF VECTOR IS NON POSITIVE'
```
```  STOP
```
```END IF
```

`er = 0`
```  n2 = 2*n
```
```  IF (d(1) < zero) THEN
```
```    ier = 1
```
```    RETURN
```
```  END IF
```

```  f(1) = SQRT(d(1))
```
```  y = one/f(1)
```
```  DO j = 2,n
```
```    f(j) = d(1+j*(j-1)/2) * y
```
```  END DO
```

```  DO i = 2,n
```
```    v = d(i*(i-1)/2+i)
```
```    DO m = 1,i-1
```
```      v = v - f((m-1)*(n2-m)/2+i)**2
```
```    END DO
```

```    IF (v < zero) THEN
```
```      ier = 1
```
```      RETURN
```
```    END IF
```

```    v = SQRT(v)
```
```    y = one/v
```
```    f((i-1)*(n2-i)/2+i) = v
```
```    DO j = i+1,n
```
```      v = d(j*(j-1)/2+i)
```
```      DO m = 1,i-1
```
```        v = v - f((m-1)*(n2-m)/2+i)*f((m-1)*(n2-m)/2 + j)
```
```      END DO ! m = 1,i-1
```
```      f((i-1)*(n2-i)/2 + j) = v*y
```
```    END DO ! j = i+1,n
```
```  END DO ! i = 2,n
```
```END IF
```

```x(1:n) = h(1:n)
```
```DO j = 1,n
```
```  y = random_normal()
```
```  DO i = j,n
```
```    x(i) = x(i) + f((j-1)*(n2-j)/2 + i) * y
```
```  END DO ! i = j,n
```
```END DO ! j = 1,n
```

```RETURN
```
```END SUBROUTINE random_mvnorm
```

```FUNCTION random_inv_gauss(h, b, first) RESULT(fn_val)
```

```! Adapted from Fortran 77 code from the book:
```
```!     Dagpunar, J. 'Principles of random variate generation'
```
```!     Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9
```

```! FUNCTION GENERATES A RANDOM VARIATE IN [0,INFINITY] FROM
```
```! A REPARAMETERISED GENERALISED INVERSE GAUSSIAN (GIG) DISTRIBUTION
```
```! WITH DENSITY PROPORTIONAL TO  GIG**(H-1) * EXP(-0.5*B*(GIG+1/GIG))
```
```! USING A RATIO METHOD.
```

```!     H = PARAMETER OF DISTRIBUTION (0 <= REAL)
```
```!     B = PARAMETER OF DISTRIBUTION (0 < REAL)
```

```REAL, INTENT(IN)    :: h, b
```
```LOGICAL, INTENT(IN) :: first
```
```REAL                :: fn_val
```

```!     Local variables
```
```REAL            :: ym, xm, r, w, r1, r2, x
```
```REAL, SAVE      :: a, c, d, e
```
```REAL, PARAMETER :: quart = 0.25
```

```  WRITE(*, *) 'IMPERMISSIBLE DISTRIBUTION PARAMETER VALUES'
```
```  STOP
```
```END IF
```

```  IF (h > quart*b*SQRT(vlarge)) THEN
```
```    WRITE(*, *) 'THE RATIO H:B IS TOO SMALL'
```
```    STOP
```
```  END IF
```
```  e = b*b
```
```  d = h + one
```
```  ym = (-d + SQRT(d*d + e))/b
```
```  IF (ym < vsmall) THEN
```
```    WRITE(*, *) 'THE VALUE OF B IS TOO SMALL'
```
```    STOP
```
```  END IF
```

```  d = h - one
```
```  xm = (d + SQRT(d*d + e))/b
```
```  d = half*d
```
```  e = -quart*b
```
```  r = xm + one/xm
```
```  w = xm*ym
```
```  a = w**(-half*h) * SQRT(xm/ym) * EXP(-e*(r - ym - one/ym))
```
```  IF (a < vsmall) THEN
```
```    WRITE(*, *) 'THE VALUE OF H IS TOO LARGE'
```
```    STOP
```
```  END IF
```
```  c = -d*LOG(xm) - e*r
```
```END IF
```

```DO
```
```  CALL RANDOM_NUMBER(r1)
```
```  IF (r1 <= zero) CYCLE
```
```  CALL RANDOM_NUMBER(r2)
```
```  x = a*r2/r1
```
```  IF (x <= zero) CYCLE
```
```  IF (LOG(r1) < d*LOG(x) + e*(x + one/x) + c) EXIT
```
```END DO
```

```fn_val = x
```

```RETURN
```
```END FUNCTION random_inv_gauss
```

```FUNCTION random_Poisson(mu, first) RESULT(ival)
```
```!**********************************************************************
```
```!     Translated to Fortran 90 by Alan Miller from:
```
```!                           RANLIB
```
```!
```
```!     Library of Fortran Routines for Random Number Generation
```
```!
```
```!                    Compiled and Written by:
```
```!
```
```!                         Barry W. Brown
```
```!                          James Lovato
```
```!
```
```!             Department of Biomathematics, Box 237
```
```!             The University of Texas, M.D. Anderson Cancer Center
```
```!             1515 Holcombe Boulevard
```
```!             Houston, TX      77030
```
```!
```
```! This work was supported by grant CA-16672 from the National Cancer Institute.
```

```!                    GENerate POIsson random deviate
```

```!                            Function
```

```! Generates a single random deviate from a Poisson distribution with mean mu.
```

```!                            Arguments
```

```!     mu --> The mean of the Poisson distribution from which
```
```!            a random deviate is to be generated.
```
```!                              REAL mu
```

```!                              Method
```

```!     For details see:
```

```!               Ahrens, J.H. and Dieter, U.
```
```!               Computer Generation of Poisson Deviates
```
```!               From Modified Normal Distributions.
```
```!               ACM Trans. Math. Software, 8, 2
```
```!               (June 1982),163-179
```

```!     TABLES: COEFFICIENTS A0-A7 FOR STEP F. FACTORIALS FACT
```
```!     COEFFICIENTS A(K) - FOR PX = FK*V*V*SUM(A(K)*V**K)-DEL
```

```!     SEPARATION OF CASES A AND B
```

```!     .. Scalar Arguments ..
```
```REAL, INTENT(IN)    :: mu
```
```LOGICAL, INTENT(IN) :: first
```
```!     ..
```
```!     .. Local Scalars ..
```
```REAL          :: b1, b2, c, c0, c1, c2, c3, del, difmuk, e, fk, fx, fy, g,  &
```
```                 omega, px, py, t, u, v, x, xx
```
```REAL, SAVE    :: s, d, p, q, p0
```
```LOGICAL, SAVE :: full_init
```
```!     ..
```
```!     .. Local Arrays ..
```
```REAL, SAVE    :: pp(35)
```
```!     ..
```
```!     .. Data statements ..
```
```REAL, PARAMETER :: a0 = -.5, a1 = .3333333, a2 = -.2500068, a3 = .2000118,  &
```
```                   a4 = -.1661269, a5 = .1421878, a6 = -.1384794,   &
```
```                   a7 = .1250060
```

```REAL, PARAMETER :: fact(10) = (/ 1., 1., 2., 6., 24., 120., 720., 5040.,  &
```
```                                 40320., 362880. /)
```

```!     ..
```
```!     .. Executable Statements ..
```
```!     C A S E  A. (RECALCULATION OF S, D, L IF MU HAS CHANGED)
```

```  IF (first) THEN
```
```    s = SQRT(mu)
```
```    d = 6.0*mu*mu
```

```!             THE POISSON PROBABILITIES PK EXCEED THE DISCRETE NORMAL
```
```!             PROBABILITIES FK WHENEVER K >= M(MU). L=IFIX(MU-1.1484)
```
```!             IS AN UPPER BOUND TO M(MU) FOR ALL MU >= 10 .
```

```    l = mu - 1.1484
```
```    full_init = .false.
```
```  END IF
```

```!     STEP N. NORMAL SAMPLE - random_normal() FOR STANDARD NORMAL DEVIATE
```

```  g = mu + s*random_normal()
```
```  IF (g > 0.0) THEN
```
```    ival = g
```

```!     STEP I. IMMEDIATE ACCEPTANCE IF ival IS LARGE ENOUGH
```

```    IF (ival>=l) RETURN
```

```!     STEP S. SQUEEZE ACCEPTANCE - SAMPLE U
```

```    fk = ival
```
```    difmuk = mu - fk
```
```    CALL RANDOM_NUMBER(u)
```
```    IF (d*u >= difmuk*difmuk*difmuk) RETURN
```
```  END IF
```

```!     STEP P. PREPARATIONS FOR STEPS Q AND H.
```
```!             (RECALCULATIONS OF PARAMETERS IF NECESSARY)
```
```!             .3989423=(2*PI)**(-.5)  .416667E-1=1./24.  .1428571=1./7.
```
```!             THE QUANTITIES B1, B2, C3, C2, C1, C0 ARE FOR THE HERMITE
```
```!             APPROXIMATIONS TO THE DISCRETE NORMAL PROBABILITIES FK.
```
```!             C=.1069/MU GUARANTEES MAJORIZATION BY THE 'HAT'-FUNCTION.
```

```  IF (.NOT. full_init) THEN
```
```    omega = .3989423/s
```
```    b1 = .4166667E-1/mu
```
```    b2 = .3*b1*b1
```
```    c3 = .1428571*b1*b2
```
```    c2 = b2 - 15.*c3
```
```    c1 = b1 - 6.*b2 + 45.*c3
```
```    c0 = 1. - b1 + 3.*b2 - 15.*c3
```
```    c = .1069/mu
```
```    full_init = .true.
```
```  END IF
```

```  IF (g < 0.0) GO TO 50
```

```!             'SUBROUTINE' F IS CALLED (KFLAG=0 FOR CORRECT RETURN)
```

```  kflag = 0
```
```  GO TO 70
```

```!     STEP Q. QUOTIENT ACCEPTANCE (RARE CASE)
```

```  40 IF (fy-u*fy <= py*EXP(px-fx)) RETURN
```

```!     STEP E. EXPONENTIAL SAMPLE - random_exponential() FOR STANDARD EXPONENTIAL
```
```!             DEVIATE E AND SAMPLE T FROM THE LAPLACE 'HAT'
```
```!             (IF T <= -.6744 THEN PK < FK FOR ALL MU >= 10.)
```

```  50 e = random_exponential()
```
```  CALL RANDOM_NUMBER(u)
```
```  u = u + u - one
```
```  t = 1.8 + SIGN(e, u)
```
```  IF (t <= (-.6744)) GO TO 50
```
```  ival = mu + s*t
```
```  fk = ival
```
```  difmuk = mu - fk
```

```!             'SUBROUTINE' F IS CALLED (KFLAG=1 FOR CORRECT RETURN)
```

```  kflag = 1
```
```  GO TO 70
```

```!     STEP H. HAT ACCEPTANCE (E IS REPEATED ON REJECTION)
```

```  60 IF (c*ABS(u) > py*EXP(px+e) - fy*EXP(fx+e)) GO TO 50
```
```  RETURN
```

```!     STEP F. 'SUBROUTINE' F. CALCULATION OF PX, PY, FX, FY.
```
```!             CASE ival < 10 USES FACTORIALS FROM TABLE FACT
```

```  70 IF (ival>=10) GO TO 80
```
```  px = -mu
```
```  py = mu**ival/fact(ival+1)
```
```  GO TO 110
```

```!             CASE ival >= 10 USES POLYNOMIAL APPROXIMATION
```
```!             A0-A7 FOR ACCURACY WHEN ADVISABLE
```
```!             .8333333E-1=1./12.  .3989423=(2*PI)**(-.5)
```

```  80 del = .8333333E-1/fk
```
```  del = del - 4.8*del*del*del
```
```  v = difmuk/fk
```
```  IF (ABS(v)>0.25) THEN
```
```    px = fk*LOG(one + v) - difmuk - del
```
```  ELSE
```
```    px = fk*v*v* (((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v+a0) - del
```
```  END IF
```
```  py = .3989423/SQRT(fk)
```
```  110 x = (half - difmuk)/s
```
```  xx = x*x
```
```  fx = -half*xx
```
```  fy = omega* (((c3*xx + c2)*xx + c1)*xx + c0)
```
```  IF (kflag <= 0) GO TO 40
```
```  GO TO 60
```

```!---------------------------------------------------------------------------
```
```!     C A S E  B.    mu < 10
```
```!     START NEW TABLE AND CALCULATE P0 IF NECESSARY
```

```ELSE
```
```  IF (first) THEN
```
```    m = MAX(1, INT(mu))
```
```    l = 0
```
```    p = EXP(-mu)
```
```    q = p
```
```    p0 = p
```
```  END IF
```

```!     STEP U. UNIFORM SAMPLE FOR INVERSION METHOD
```

```  DO
```
```    CALL RANDOM_NUMBER(u)
```
```    ival = 0
```
```    IF (u <= p0) RETURN
```

```!     STEP T. TABLE COMPARISON UNTIL THE END PP(L) OF THE
```
```!             PP-TABLE OF CUMULATIVE POISSON PROBABILITIES
```
```!             (0.458=PP(9) FOR MU=10)
```

```    IF (l == 0) GO TO 150
```
```    j = 1
```
```    IF (u > 0.458) j = MIN(l, m)
```
```    DO k = j, l
```
```      IF (u <= pp(k)) GO TO 180
```
```    END DO
```
```    IF (l == 35) CYCLE
```

```!     STEP C. CREATION OF NEW POISSON PROBABILITIES P
```
```!             AND THEIR CUMULATIVES Q=PP(K)
```

```    150 l = l + 1
```
```    DO k = l, 35
```
```      p = p*mu / k
```
```      q = q + p
```
```      pp(k) = q
```
```      IF (u <= q) GO TO 170
```
```    END DO
```
```    l = 35
```
```  END DO
```

```  170 l = k
```
```  180 ival = k
```
```  RETURN
```
```END IF
```

```RETURN
```
```END FUNCTION random_Poisson
```

```FUNCTION random_binomial1(n, p, first) RESULT(ival)
```

```! FUNCTION GENERATES A RANDOM BINOMIAL VARIATE USING C.D.Kemp's method.
```
```! This algorithm is suitable when many random variates are required
```
```! with the SAME parameter values for n & p.
```

```!    P = BERNOULLI SUCCESS PROBABILITY
```
```!           (0 <= REAL <= 1)
```
```!    N = NUMBER OF BERNOULLI TRIALS
```
```!           (1 <= INTEGER)
```
```!    FIRST = .TRUE. for the first call using the current parameter values
```
```!          = .FALSE. if the values of (n,p) are unchanged from last call
```

```! Reference: Kemp, C.D. (1986). `A modal method for generating binomial
```
```!            variables', Commun. Statist. - Theor. Meth. 15(3), 805-813.
```

```REAL, INTENT(IN)    :: p
```
```LOGICAL, INTENT(IN) :: first
```

```!     Local variables
```

```REAL            :: u, pd, pu
```
```REAL, SAVE      :: odds_ratio, p_r
```
```REAL, PARAMETER :: zero = 0.0, one = 1.0
```

```  r0 = (n+1)*p
```
```  p_r = bin_prob(n, p, r0)
```
```  odds_ratio = p / (one - p)
```
```END IF
```

```CALL RANDOM_NUMBER(u)
```
```u = u - p_r
```
```  ival = r0
```
```  RETURN
```
```END IF
```

```DO
```
```  rd = rd - 1
```
```  IF (rd >= 0) THEN
```
```    pd = pd * (rd+1) / (odds_ratio * (n-rd))
```
```    u = u - pd
```
```    IF (u < zero) THEN
```
```      ival = rd
```
```      RETURN
```
```    END IF
```
```  END IF
```

```  ru = ru + 1
```
```  IF (ru <= n) THEN
```
```    pu = pu * (n-ru+1) * odds_ratio / ru
```
```    u = u - pu
```
```    IF (u < zero) THEN
```
```      ival = ru
```
```      RETURN
```
```    END IF
```
```  END IF
```
```END DO
```

```!     This point should not be reached, but just in case:
```

`val = r0`
```RETURN
```

```END FUNCTION random_binomial1
```

```FUNCTION bin_prob(n, p, r) RESULT(fn_val)
```
```!     Calculate a binomial probability
```

```REAL, INTENT(IN)    :: p
```
```REAL                :: fn_val
```

```!     Local variable
```
```REAL                :: one = 1.0
```

```fn_val = EXP( lngamma(DBLE(n+1)) - lngamma(DBLE(r+1)) - lngamma(DBLE(n-r+1)) &
```
```              + r*LOG(p) + (n-r)*LOG(one - p) )
```
```RETURN
```

```END FUNCTION bin_prob
```

```FUNCTION lngamma(x) RESULT(fn_val)
```

```! Logarithm to base e of the gamma function.
```
```!
```
```! Accurate to about 1.e-14.
```
```! Programmer: Alan Miller
```

```! Latest revision of Fortran 77 version - 28 February 1988
```

```REAL (dp), INTENT(IN) :: x
```
```REAL (dp)             :: fn_val
```

```!       Local variables
```

```REAL (dp) :: a1 = -4.166666666554424D-02, a2 = 2.430554511376954D-03,  &
```
```             a3 = -7.685928044064347D-04, a4 = 5.660478426014386D-04,  &
```
```             temp, arg, product, lnrt2pi = 9.189385332046727D-1,       &
```
```             pi = 3.141592653589793D0
```
```LOGICAL   :: reflect
```

```!       lngamma is not defined if x = 0 or a negative integer.
```

```fn_val = 0.d0
```
```RETURN
```

```!       If x < 0, use the reflection formula:
```
```!               gamma(x) * gamma(1-x) = pi * cosec(pi.x)
```

```  arg = 1.d0 - x
```
```ELSE
```
```  arg = x
```
```END IF
```

```!       Increase the argument, if necessary, to make it > 10.
```

```  product = product * arg
```
```  arg = arg + 1.d0
```
```  GO TO 20
```
```END IF
```

```!  Use a polynomial approximation to Stirling's formula.
```
```!  N.B. The real Stirling's formula is used here, not the simpler, but less
```
```!       accurate formula given by De Moivre in a letter to Stirling, which
```
```!       is the one usually quoted.
```

```arg = arg - 0.5D0
```
```temp = 1.d0/arg**2
```
```fn_val = lnrt2pi + arg * (LOG(arg) - 1.d0 + &
```
```                  (((a4*temp + a3)*temp + a2)*temp + a1)*temp) - LOG(product)
```
```  temp = SIN(pi * x)
```
```  fn_val = LOG(pi/temp) - fn_val
```
```END IF
```
```RETURN
```
```END FUNCTION lngamma
```

```FUNCTION random_binomial2(n, pp, first) RESULT(ival)
```
```!**********************************************************************
```
```!     Translated to Fortran 90 by Alan Miller from:
```
```!                              RANLIB
```
```!
```
```!     Library of Fortran Routines for Random Number Generation
```
```!
```
```!                      Compiled and Written by:
```
```!
```
```!                           Barry W. Brown
```
```!                            James Lovato
```
```!
```
```!               Department of Biomathematics, Box 237
```
```!               The University of Texas, M.D. Anderson Cancer Center
```
```!               1515 Holcombe Boulevard
```
```!               Houston, TX      77030
```
```!
```
```! This work was supported by grant CA-16672 from the National Cancer Institute.
```

```!                    GENerate BINomial random deviate
```

```!                              Function
```

```!     Generates a single random deviate from a binomial
```
```!     distribution whose number of trials is N and whose
```
```!     probability of an event in each trial is P.
```

```!                              Arguments
```

```!     N  --> The number of trials in the binomial distribution
```
```!            from which a random deviate is to be generated.
```
```!                              INTEGER N
```

```!     P  --> The probability of an event in each trial of the
```
```!            binomial distribution from which a random deviate
```
```!            is to be generated.
```
```!                              REAL P
```

```!     FIRST --> Set FIRST = .TRUE. for the first call to perform initialization
```
```!               the set FIRST = .FALSE. for further calls using the same pair
```
```!               of parameter values (N, P).
```
```!                              LOGICAL FIRST
```

```!     random_binomial2 <-- A random deviate yielding the number of events
```
```!                from N independent trials, each of which has
```
```!                a probability of event P.
```
```!                              INTEGER random_binomial
```

```!                              Method
```

```!     This is algorithm BTPE from:
```

```!         Kachitvichyanukul, V. and Schmeiser, B. W.
```
```!         Binomial Random Variate Generation.
```
```!         Communications of the ACM, 31, 2 (February, 1988) 216.
```

```!**********************************************************************
```

```!*****DETERMINE APPROPRIATE ALGORITHM AND WHETHER SETUP IS NECESSARY
```

```!     ..
```
```!     .. Scalar Arguments ..
```
```REAL, INTENT(IN)    :: pp
```
```LOGICAL, INTENT(IN) :: first
```
```!     ..
```
```!     .. Local Scalars ..
```
```REAL            :: alv, amaxp, f, f1, f2, u, v, w, w2, x, x1, x2, ynorm, z, z2
```
```REAL, PARAMETER :: zero = 0.0, half = 0.5, one = 1.0
```
```REAL, SAVE      :: p, q, xnp, ffm, fm, xnpq, p1, xm, xl, xr, c, al, xll,  &
```
```                   xlr, p2, p3, p4, qn, r, g
```

```!     ..
```
```!     .. Executable Statements ..
```

```!*****SETUP, PERFORM ONLY WHEN PARAMETERS CHANGE
```

```  p = MIN(pp, one-pp)
```
```  q = one - p
```
```  xnp = n * p
```
```END IF
```

```  IF (first) THEN
```
```    ffm = xnp + p
```
```    m = ffm
```
```    fm = m
```
```    xnpq = xnp * q
```
```    p1 = INT(2.195*SQRT(xnpq) - 4.6*q) + half
```
```    xm = fm + half
```
```    xl = xm - p1
```
```    xr = xm + p1
```
```    c = 0.134 + 20.5 / (15.3 + fm)
```
```    al = (ffm-xl) / (ffm - xl*p)
```
```    xll = al * (one + half*al)
```
```    al = (xr - ffm) / (xr*q)
```
```    xlr = al * (one + half*al)
```
```    p2 = p1 * (one + c + c)
```
```    p3 = p2 + c / xll
```
```    p4 = p3 + c / xlr
```
```  END IF
```

```!*****GENERATE VARIATE, Binomial mean at least 30.
```

```  20 CALL RANDOM_NUMBER(u)
```
```  u = u * p4
```
```  CALL RANDOM_NUMBER(v)
```

```!     TRIANGULAR REGION
```

```  IF (u <= p1) THEN
```
```    ix = xm - p1 * v + u
```
```    GO TO 110
```
```  END IF
```

```!     PARALLELOGRAM REGION
```

```  IF (u <= p2) THEN
```
```    x = xl + (u-p1) / c
```
```    v = v * c + one - ABS(xm-x) / p1
```
```    IF (v > one .OR. v <= zero) GO TO 20
```
```    ix = x
```
```  ELSE
```

```!     LEFT TAIL
```

```    IF (u <= p3) THEN
```
```      ix = xl + LOG(v) / xll
```
```      IF (ix < 0) GO TO 20
```
```      v = v * (u-p2) * xll
```
```    ELSE
```

```!     RIGHT TAIL
```

```      ix = xr - LOG(v) / xlr
```
```      IF (ix > n) GO TO 20
```
```      v = v * (u-p3) * xlr
```
```    END IF
```
```  END IF
```

```!*****DETERMINE APPROPRIATE WAY TO PERFORM ACCEPT/REJECT TEST
```

```  k = ABS(ix-m)
```
```  IF (k <= 20 .OR. k >= xnpq/2-1) THEN
```

```!     EXPLICIT EVALUATION
```

```    f = one
```
```    r = p / q
```
```    g = (n+1) * r
```
```    IF (m < ix) THEN
```
```      mp = m + 1
```
```      DO i = mp, ix
```
```        f = f * (g/i-r)
```
```      END DO
```

```    ELSE IF (m > ix) THEN
```
```      ix1 = ix + 1
```
```      DO i = ix1, m
```
```        f = f / (g/i-r)
```
```      END DO
```
```    END IF
```

```    IF (v > f) THEN
```
```      GO TO 20
```
```    ELSE
```
```      GO TO 110
```
```    END IF
```
```  END IF
```

```!     SQUEEZING USING UPPER AND LOWER BOUNDS ON LOG(F(X))
```

```  amaxp = (k/xnpq) * ((k*(k/3. + .625) + .1666666666666)/xnpq + half)
```
```  ynorm = -k * k / (2.*xnpq)
```
```  alv = LOG(v)
```
```  IF (alv<ynorm - amaxp) GO TO 110
```
```  IF (alv>ynorm + amaxp) GO TO 20
```

```!     STIRLING'S (actually de Moivre's) FORMULA TO MACHINE ACCURACY FOR
```
```!     THE FINAL ACCEPTANCE/REJECTION TEST
```

```  x1 = ix + 1
```
```  f1 = fm + one
```
```  z = n + 1 - fm
```
```  w = n - ix + one
```
```  z2 = z * z
```
```  x2 = x1 * x1
```
```  f2 = f1 * f1
```
```  w2 = w * w
```
```  IF (alv - (xm*LOG(f1/x1) + (n-m+half)*LOG(z/w) + (ix-m)*LOG(w*p/(x1*q)) +    &
```
```      (13860.-(462.-(132.-(99.-140./f2)/f2)/f2)/f2)/f1/166320. +               &
```
```      (13860.-(462.-(132.-(99.-140./z2)/z2)/z2)/z2)/z/166320. +                &
```
```      (13860.-(462.-(132.-(99.-140./x2)/x2)/x2)/x2)/x1/166320. +               &
```
```      (13860.-(462.-(132.-(99.-140./w2)/w2)/w2)/w2)/w/166320.) > zero) THEN
```
```    GO TO 20
```
```  ELSE
```
```    GO TO 110
```
```  END IF
```

```ELSE
```
```!     INVERSE CDF LOGIC FOR MEAN LESS THAN 30
```
```  IF (first) THEN
```
```    qn = q ** n
```
```    r = p / q
```
```    g = r * (n+1)
```
```  END IF
```

```  90 ix = 0
```
```  f = qn
```
```  CALL RANDOM_NUMBER(u)
```
```  100 IF (u >= f) THEN
```
```    IF (ix > 110) GO TO 90
```
```    u = u - f
```
```    ix = ix + 1
```
```    f = f * (g/ix - r)
```
```    GO TO 100
```
```  END IF
```
```END IF
```

`val = ix`
```RETURN
```

```END FUNCTION random_binomial2
```

```FUNCTION random_neg_binomial(sk, p) RESULT(ival)
```

```! Adapted from Fortran 77 code from the book:
```
```!     Dagpunar, J. 'Principles of random variate generation'
```
```!     Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9
```

```! FUNCTION GENERATES A RANDOM NEGATIVE BINOMIAL VARIATE USING UNSTORED
```
```! INVERSION AND/OR THE REPRODUCTIVE PROPERTY.
```

```!    SK = NUMBER OF FAILURES REQUIRED (Dagpunar's words!)
```
```!       = the `power' parameter of the negative binomial
```
```!           (0 < REAL)
```
```!    P = BERNOULLI SUCCESS PROBABILITY
```
```!           (0 < REAL < 1)
```

```! THE PARAMETER H IS SET SO THAT UNSTORED INVERSION ONLY IS USED WHEN P <= H,
```
```! OTHERWISE A COMBINATION OF UNSTORED INVERSION AND
```
```! THE REPRODUCTIVE PROPERTY IS USED.
```

```REAL, INTENT(IN)   :: sk, p
```

```!     Local variables
```
```! THE PARAMETER ULN = -LOG(MACHINE'S SMALLEST REAL NUMBER).
```

```REAL, PARAMETER    :: h = 0.7
```
```REAL               :: q, x, st, uln, v, r, s, y, g
```

```  WRITE(*, *) 'IMPERMISSIBLE DISTRIBUTION PARAMETER VALUES'
```
```  STOP
```
```END IF
```

```q = one - p
```
```x = zero
```
```  v = one/LOG(p)
```
```  k = st
```
```  DO i = 1,k
```
```    DO
```
```      CALL RANDOM_NUMBER(r)
```
```      IF (r > zero) EXIT
```
```    END DO
```
```    n = v*LOG(r)
```
```    x = x + n
```
```  END DO
```
```  st = st - k
```
```END IF
```

```uln = -LOG(vsmall)
```
```  WRITE(*, *) ' P IS TOO LARGE FOR THIS VALUE OF SK'
```
```  STOP
```
```END IF
```

```y = q**st
```
```CALL RANDOM_NUMBER(r)
```
```DO
```
```  IF (y > r) EXIT
```
```  r = r - y
```
```  s = s + one
```
```  y = y*p*g/s
```
```  g = g + one
```
```END DO
```

`val = x + s + half`
```RETURN
```
```END FUNCTION random_neg_binomial
```

```FUNCTION random_von_Mises(k, first) RESULT(fn_val)
```

```!     Algorithm VMD from:
```
```!     Dagpunar, J.S. (1990) `Sampling from the von Mises distribution via a
```
```!     comparison of random numbers', J. of Appl. Statist., 17, 165-168.
```

```!     Fortran 90 code by Alan Miller
```
```!     CSIRO Division of Mathematical & Information Sciences
```

```!     Arguments:
```
```!     k (real)        parameter of the von Mises distribution.
```
```!     first (logical) set to .TRUE. the first time that the function
```
```!                     is called, or the first time with a new value
```
```!                     for k.   When first = .TRUE., the function sets
```
```!                     up starting values and may be very much slower.
```

```REAL, INTENT(IN)     :: k
```
```LOGICAL, INTENT(IN)  :: first
```
```REAL                 :: fn_val
```

```!     Local variables
```

```REAL, PARAMETER  :: pi = 3.14159265
```
```REAL, SAVE       :: p(20), theta(0:20)
```
```REAL             :: sump, r, th, lambda, rlast
```
```REAL (dp)        :: dk
```

```  IF (k < zero) THEN
```
```    WRITE(*, *) '** Error: argument k for random_von_Mises = ', k
```
```    RETURN
```
```  END IF
```

```  nk = k + k + one
```
```  IF (nk > 20) THEN
```
```    WRITE(*, *) '** Error: argument k for random_von_Mises = ', k
```
```    RETURN
```
```  END IF
```

```  dk = k
```
```  theta(0) = zero
```
```  IF (k > half) THEN
```

```!     Set up array p of probabilities.
```

```    sump = zero
```
```    DO j = 1, nk
```
```      IF (j < nk) THEN
```
```        theta(j) = ACOS(one - j/k)
```
```      ELSE
```
```        theta(nk) = pi
```
```      END IF
```

```!     Numerical integration of e^[k.cos(x)] from theta(j-1) to theta(j)
```

```      CALL integral(theta(j-1), theta(j), p(j), dk)
```
```      sump = sump + p(j)
```
```    END DO
```
```    p(1:nk) = p(1:nk) / sump
```
```  ELSE
```
```    p(1) = one
```
```    theta(1) = pi
```
```  END IF                         ! if k > 0.5
```
```END IF                           ! if first
```

```CALL RANDOM_NUMBER(r)
```
```DO j = 1, nk
```
```  r = r - p(j)
```
```  IF (r < zero) EXIT
```
```END DO
```

```DO
```
```  th = theta(j-1) + r*(theta(j) - theta(j-1))
```
```  lambda = k - j + one - k*COS(th)
```
```  n = 1
```
```  rlast = lambda
```

```  DO
```
```    CALL RANDOM_NUMBER(r)
```
```    IF (r > rlast) EXIT
```
```    n = n + 1
```
```    rlast = r
```
```  END DO
```

```  IF (n .NE. 2*(n/2)) EXIT         ! is n even?
```
```  CALL RANDOM_NUMBER(r)
```
```END DO
```

```fn_val = SIGN(th, (r - rlast)/(one - rlast) - half)
```
```RETURN
```
```END FUNCTION random_von_Mises
```

```SUBROUTINE integral(a, b, result, dk)
```

```!     Gaussian integration of exp(k.cosx) from a to b.
```

```REAL (dp), INTENT(IN) :: dk
```
```REAL, INTENT(IN)      :: a, b
```
```REAL, INTENT(OUT)     :: result
```

```!     Local variables
```

```REAL (dp)  :: xmid, range, x1, x2,                                    &
```
```  x(3) = (/0.238619186083197_dp, 0.661209386466265_dp, 0.932469514203152_dp/), &
```
```  w(3) = (/0.467913934572691_dp, 0.360761573048139_dp, 0.171324492379170_dp/)
```

```xmid = (a + b)/2._dp
```

```DO i = 1, 3
```
```  x1 = xmid + x(i)*range
```
```  x2 = xmid - x(i)*range
```
```  result = result + w(i)*(EXP(dk*COS(x1)) + EXP(dk*COS(x2)))
```
```END DO
```

```RETURN
```
```END SUBROUTINE integral
```

```FUNCTION random_Cauchy() RESULT(fn_val)
```

```!     Generate a random deviate from the standard Cauchy distribution
```

```REAL     :: fn_val
```

```!     Local variables
```
```REAL     :: v(2)
```

```DO
```
```  CALL RANDOM_NUMBER(v)
```
```  v = two*(v - half)
```
```  IF (ABS(v(2)) < vsmall) CYCLE               ! Test for zero
```
```  IF (v(1)**2 + v(2)**2 < one) EXIT
```
```END DO
```
```fn_val = v(1) / v(2)
```

```RETURN
```
```END FUNCTION random_Cauchy
```

```SUBROUTINE random_order(order, n)
```

```!     Generate a random ordering of the integers 1 ... n.
```

```!     Local variables
```

```REAL    :: wk
```

```DO i = 1, n
```
```  order(i) = i
```
```END DO
```

```!     Starting at the end, swap the current last indicator with one
```
```!     randomly chosen from those preceeding it.
```

```DO i = n, 2, -1
```
```  CALL RANDOM_NUMBER(wk)
```
```  j = 1 + i * wk
```
```  IF (j < i) THEN
```
```    k = order(i)
```
```    order(i) = order(j)
```
```    order(j) = k
```
```  END IF
```
```END DO
```

```RETURN
```
```END SUBROUTINE random_order
```

```SUBROUTINE seed_random_number(iounit)
```

```! Local variables
```

```CALL RANDOM_SEED(SIZE=k)
```
```ALLOCATE( seed(k) )
```

```WRITE(*, '(a, i2, a)')' Enter ', k, ' integers for random no. seeds: '
```
```READ(*, *) seed
```
```WRITE(iounit, '(a, (7i10))') ' Random no. seeds: ', seed
```
```CALL RANDOM_SEED(PUT=seed)
```

```DEALLOCATE( seed )
```

```RETURN
```
```END SUBROUTINE seed_random_number
```

```END MODULE random
```

```.
```