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.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real, | intent(in) | :: | mu | |||
| logical, | intent(in) | :: | first |
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