d0/d25/ignpoi_8f_source.html

       INTEGER*4 FUNCTION ignpoi(mu)

 C**********************************************************************

 C

 C     INTEGER*4 FUNCTION IGNPOI( MU )

 C

 C                    GENerate POIsson random deviate

 C

 C

 C                              Function

 C

 C

 C     Generates a single random deviate from a Poisson

 C     distribution with mean MU.

 C

 C

 C                              Arguments

 C

 C

 C     MU --> The mean of the Poisson distribution from which

 C            a random deviate is to be generated.

 C                              REAL MU

 C     JJV                    (MU >= 0.0)

 C

 C     IGNPOI <-- The random deviate.

 C                              INTEGER IGNPOI (non-negative)

 C

 C

 C                              Method

 C

 C

 C     Renames KPOIS from TOMS as slightly modified by BWB to use RANF

 C     instead of SUNIF.

 C

 C     For details see:

 C

 C               Ahrens, J.H. and Dieter, U.

 C               Computer Generation of Poisson Deviates

 C               From Modified Normal Distributions.

 C               ACM Trans. Math. Software, 8, 2

 C               (June 1982),163-179

 C

 C**********************************************************************

 C**********************************************************************C

 C**********************************************************************C

 C                                                                      C

 C                                                                      C

 C     P O I S S O N  DISTRIBUTION                                      C

 C                                                                      C

 C                                                                      C

 C**********************************************************************C

 C**********************************************************************C

 C                                                                      C

 C     FOR DETAILS SEE:                                                 C

 C                                                                      C

 C               AHRENS, J.H. AND DIETER, U.                            C

 C               COMPUTER GENERATION OF POISSON DEVIATES                C

 C               FROM MODIFIED NORMAL DISTRIBUTIONS.                    C

 C               ACM TRANS. MATH. SOFTWARE, 8,2 (JUNE 1982), 163 - 179. C

 C                                                                      C

 C     (SLIGHTLY MODIFIED VERSION OF THE PROGRAM IN THE ABOVE ARTICLE)  C

 C                                                                      C

 C**********************************************************************C

 C

 C      INTEGER*4 FUNCTION IGNPOI(IR,MU)

 C

 C     INPUT:  IR=CURRENT STATE OF BASIC RANDOM NUMBER GENERATOR

 C             MU=MEAN MU OF THE POISSON DISTRIBUTION

 C     OUTPUT: IGNPOI=SAMPLE FROM THE POISSON-(MU)-DISTRIBUTION

 C

 C

 C

 C     MUPREV=PREVIOUS MU, MUOLD=MU AT LAST EXECUTION OF STEP P OR CASE B

 C     TABLES: COEFFICIENTS A0-A7 FOR STEP F. FACTORIALS FACT

 C     COEFFICIENTS A(K) - FOR PX = FK*V*V*SUM(A(K)*V**K)-DEL

 C

 C

 C

 C     SEPARATION OF CASES A AND B

 C

 C     .. Scalar Arguments ..

       REAL mu

 C     ..

 C     .. Local Scalars ..

       REAL a0,a1,a2,a3,a4,a5,a6,a7,b1,b2,c,c0,c1,c2,c3,d,del,difmuk,e,

      +     fk,fx,fy,g,muold,muprev,omega,p,p0,px,py,q,s,t,u,v,x,xx

 C     JJV I added a variable 'll' here - it is the 'l' for CASE A

       INTEGER*4 j,k,kflag,l,ll,m

 C     ..

 C     .. Local Arrays ..

       REAL fact(10),pp(35)

 C     ..

 C     .. External Functions ..

       REAL ranf,sexpo,snorm

       EXTERNAL ranf,sexpo,snorm

 C     ..

 C     .. Intrinsic Functions ..

       INTRINSIC abs,alog,exp,float,ifix,max0,min0,sign,sqrt

 C     ..

 C     JJV added this for case: mu unchanged

 C     .. Save statement ..

       SAVE s, d, l, ll, omega, c3, c2, c1, c0, c, m, p, q, p0,

      +     a0, a1, a2, a3, a4, a5, a6, a7, fact, pp, muprev, muold

 C     ..

 C     JJV end addition - I am including vars in Data statements

 C     .. Data statements ..

 C     JJV changed initial values of MUPREV and MUOLD to -1.0E37

 C     JJV if no one calls IGNPOI with MU = -1.0E37 the first time,

 C     JJV the code shouldn't break

       DATA muprev,muold/-1.0e37,-1.0e37/

       DATA a0,a1,a2,a3,a4,a5,a6,a7/-.5,.3333333,-.2500068,.2000118,

      +     -.1661269,.1421878,-.1384794,.1250060/

       DATA fact/1.,1.,2.,6.,24.,120.,720.,5040.,40320.,362880./

       DATA pp/35*0.0/

 C     ..

 C     .. Executable Statements ..


       IF (mu.EQ.muprev) GO TO 10

       IF (mu.LT.10.0) GO TO 120

 C

 C     C A S E  A. (RECALCULATION OF S,D,LL IF MU HAS CHANGED)

 C

 C     JJV This is the case where I changed 'l' to 'll'

 C     JJV Here 'll' is set once and used in a comparison once


       muprev = mu

       s = sqrt(mu)

       d = 6.0*mu*mu

 C

 C             THE POISSON PROBABILITIES PK EXCEED THE DISCRETE NORMAL

 C             PROBABILITIES FK WHENEVER K >= M(MU). LL=IFIX(MU-1.1484)

 C             IS AN UPPER BOUND TO M(MU) FOR ALL MU >= 10 .

 C

       ll = ifix(mu-1.1484)

 C

 C     STEP N. NORMAL SAMPLE - SNORM(IR) FOR STANDARD NORMAL DEVIATE

 C

    10 g = mu + s*snorm()

       IF (g.LT.0.0) GO TO 20

       ignpoi = ifix(g)

 C

 C     STEP I. IMMEDIATE ACCEPTANCE IF IGNPOI IS LARGE ENOUGH

 C

       IF (ignpoi.GE.ll) RETURN

 C

 C     STEP S. SQUEEZE ACCEPTANCE - SUNIF(IR) FOR (0,1)-SAMPLE U

 C

       fk = float(ignpoi)

       difmuk = mu - fk

       u = ranf()

       IF (d*u.GE.difmuk*difmuk*difmuk) RETURN

 C

 C     STEP P. PREPARATIONS FOR STEPS Q AND H.

 C             (RECALCULATIONS OF PARAMETERS IF NECESSARY)

 C             .3989423=(2*PI)**(-.5)  .416667E-1=1./24.  .1428571=1./7.

 C             THE QUANTITIES B1, B2, C3, C2, C1, C0 ARE FOR THE HERMITE

 C             APPROXIMATIONS TO THE DISCRETE NORMAL PROBABILITIES FK.

 C             C=.1069/MU GUARANTEES MAJORIZATION BY THE 'HAT'-FUNCTION.

 C

    20 IF (mu.EQ.muold) GO TO 30

       muold = mu

       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

    30 IF (g.LT.0.0) GO TO 50

 C

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

 C

       kflag = 0

       GO TO 70

 C

 C     STEP Q. QUOTIENT ACCEPTANCE (RARE CASE)

 C

    40 IF (fy-u*fy.LE.py*exp(px-fx)) RETURN

 C

 C     STEP E. EXPONENTIAL SAMPLE - SEXPO(IR) FOR STANDARD EXPONENTIAL

 C             DEVIATE E AND SAMPLE T FROM THE LAPLACE 'HAT'

 C             (IF T <= -.6744 THEN PK < FK FOR ALL MU >= 10.)

 C

    50 e = sexpo()

       u = ranf()

       u = u + u - 1.0

       t = 1.8 + sign(e,u)

       IF (t.LE. (-.6744)) GO TO 50

       ignpoi = ifix(mu+s*t)

       fk = float(ignpoi)

       difmuk = mu - fk

 C

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

 C

       kflag = 1

       GO TO 70

 C

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

 C

    60 IF (c*abs(u).GT.py*exp(px+e)-fy*exp(fx+e)) GO TO 50

       RETURN

 C

 C     STEP F. 'SUBROUTINE' F. CALCULATION OF PX,PY,FX,FY.

 C             CASE IGNPOI .LT. 10 USES FACTORIALS FROM TABLE FACT

 C

    70 IF (ignpoi.GE.10) GO TO 80

       px = -mu

       py = mu**ignpoi/fact(ignpoi+1)

       GO TO 110

 C

 C             CASE IGNPOI .GE. 10 USES POLYNOMIAL APPROXIMATION

 C             A0-A7 FOR ACCURACY WHEN ADVISABLE

 C             .8333333E-1=1./12.  .3989423=(2*PI)**(-.5)

 C

    80 del = .8333333e-1/fk

       del = del - 4.8*del*del*del

       v = difmuk/fk

       IF (abs(v).LE.0.25) GO TO 90

       px = fk*alog(1.0+v) - difmuk - del

       GO TO 100


    90 px = fk*v*v* (((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v+a0) -

      +     del

   100 py = .3989423/sqrt(fk)

   110 x = (0.5-difmuk)/s

       xx = x*x

       fx = -0.5*xx

       fy = omega* (((c3*xx+c2)*xx+c1)*xx+c0)

       IF (kflag.LE.0) GO TO 40

       GO TO 60

 C

 C     C A S E  B. (START NEW TABLE AND CALCULATE P0 IF NECESSARY)

 C

 C     JJV changed MUPREV assignment from 0.0 to initial value

   120 muprev = -1.0e37

       IF (mu.EQ.muold) GO TO 130

 C     JJV added argument checker here

       IF (mu.GE.0.0) GO TO 125

       WRITE (*,*) 'MU < 0 in IGNPOI - ABORT'

       WRITE (*,*) 'Value of MU: ',mu

       CALL xstopx ('MU < 0 in IGNPOI - ABORT')

 C     JJV added line label here

  125  muold = mu

       m = max0(1,ifix(mu))

       l = 0

       p = exp(-mu)

       q = p

       p0 = p

 C

 C     STEP U. UNIFORM SAMPLE FOR INVERSION METHOD

 C

   130 u = ranf()

       ignpoi = 0

       IF (u.LE.p0) RETURN

 C

 C     STEP T. TABLE COMPARISON UNTIL THE END PP(L) OF THE

 C             PP-TABLE OF CUMULATIVE POISSON PROBABILITIES

 C             (0.458=PP(9) FOR MU=10)

 C

       IF (l.EQ.0) GO TO 150

       j = 1

       IF (u.GT.0.458) j = min0(l,m)

       DO 140 k = j,l

           IF (u.LE.pp(k)) GO TO 180

   140 CONTINUE

       IF (l.EQ.35) GO TO 130

 C

 C     STEP C. CREATION OF NEW POISSON PROBABILITIES P

 C             AND THEIR CUMULATIVES Q=PP(K)

 C

   150 l = l + 1

       DO 160 k = l,35

           p = p*mu/float(k)

           q = q + p

           pp(k) = q

           IF (u.LE.q) GO TO 170

   160 CONTINUE

       l = 35

       GO TO 130


   170 l = k

   180 ignpoi = k

       RETURN


       END

ignpoi
integer *4 function ignpoi(mu)
Definition: ignpoi.f:2

abs
static T abs(T x)
Definition: pr-output.cc:1678

ranf
real function ranf()
Definition: ranf.f:2

sexpo
real function sexpo()
Definition: sexpo.f:2

snorm
real function snorm()
Definition: snorm.f:2