1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
|
INTEGER FUNCTION ignpoi(mu)
C**********************************************************************
C
C INTEGER 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 DOUBLE PRECISION 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 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 ..
DOUBLE PRECISION mu
C ..
C .. Local Scalars ..
DOUBLE PRECISION 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 j,k,kflag,l,ll,m
C ..
C .. Local Arrays ..
DOUBLE PRECISION fact(10),pp(35)
C ..
C .. External Functions ..
DOUBLE PRECISION ranf,sexpo,snorm
EXTERNAL ranf,sexpo,snorm
C ..
C .. Intrinsic Functions ..
INTRINSIC abs,log,exp,float,idint,max0,min0,dsign,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, 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./
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=IDINT(MU-1.1484)
C IS AN UPPER BOUND TO M(MU) FOR ALL MU >= 10 .
C
ll = idint(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 = idint(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 = dble(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.8d0 + dsign(e,u)
IF (t.LE. (-.6744)) GO TO 50
ignpoi = idint(mu+s*t)
fk = dble(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*log(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) 40,40,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
C Jpc 1999: the next lines seams to produce a bug
C and I finaly commented them out
C IF (mu.EQ.muold) GO TO 130
C JJV added argument checker here
C JJV added line label here
C 125 muold = mu
m = max0(1,idint(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
|
