random_Poisson Function

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
  INTEGER             :: ival
  !     ..
  !     .. 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
  INTEGER       :: j, k, kflag
  LOGICAL, SAVE :: full_init
  INTEGER, SAVE :: l, m
  !     ..
  !     .. 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 ..
  IF (mu > 10.0) THEN
     !     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