GNU Octave  3.8.0 A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
ignpoi.f
Go to the documentation of this file.
1  INTEGER FUNCTION ignpoi(mu)
2 C**********************************************************************
3 C
4 C INTEGER FUNCTION IGNPOI( MU )
5 C
6 C GENerate POIsson random deviate
7 C
8 C
9 C Function
10 C
11 C
12 C Generates a single random deviate from a Poisson
13 C distribution with mean MU.
14 C
15 C
16 C Arguments
17 C
18 C
19 C MU --> The mean of the Poisson distribution from which
20 C a random deviate is to be generated.
21 C REAL MU
22 C JJV (MU >= 0.0)
23 C
24 C IGNPOI <-- The random deviate.
25 C INTEGER IGNPOI (non-negative)
26 C
27 C
28 C Method
29 C
30 C
31 C Renames KPOIS from TOMS as slightly modified by BWB to use RANF
33 C
34 C For details see:
35 C
36 C Ahrens, J.H. and Dieter, U.
37 C Computer Generation of Poisson Deviates
38 C From Modified Normal Distributions.
39 C ACM Trans. Math. Software, 8, 2
40 C (June 1982),163-179
41 C
42 C**********************************************************************
43 C**********************************************************************C
44 C**********************************************************************C
45 C C
46 C C
47 C P O I S S O N DISTRIBUTION C
48 C C
49 C C
50 C**********************************************************************C
51 C**********************************************************************C
52 C C
53 C FOR DETAILS SEE: C
54 C C
55 C AHRENS, J.H. AND DIETER, U. C
56 C COMPUTER GENERATION OF POISSON DEVIATES C
57 C FROM MODIFIED NORMAL DISTRIBUTIONS. C
58 C ACM TRANS. MATH. SOFTWARE, 8,2 (JUNE 1982), 163 - 179. C
59 C C
60 C (SLIGHTLY MODIFIED VERSION OF THE PROGRAM IN THE ABOVE ARTICLE) C
61 C C
62 C**********************************************************************C
63 C
64 C INTEGER FUNCTION IGNPOI(IR,MU)
65 C
66 C INPUT: IR=CURRENT STATE OF BASIC RANDOM NUMBER GENERATOR
67 C MU=MEAN MU OF THE POISSON DISTRIBUTION
68 C OUTPUT: IGNPOI=SAMPLE FROM THE POISSON-(MU)-DISTRIBUTION
69 C
70 C
71 C
72 C MUPREV=PREVIOUS MU, MUOLD=MU AT LAST EXECUTION OF STEP P OR CASE B
73 C TABLES: COEFFICIENTS A0-A7 FOR STEP F. FACTORIALS FACT
74 C COEFFICIENTS A(K) - FOR PX = FK*V*V*SUM(A(K)*V**K)-DEL
75 C
76 C
77 C
78 C SEPARATION OF CASES A AND B
79 C
80 C .. Scalar Arguments ..
81  REAL mu
82 C ..
83 C .. Local Scalars ..
84  REAL a0,a1,a2,a3,a4,a5,a6,a7,b1,b2,c,c0,c1,c2,c3,d,del,difmuk,e,
85  + fk,fx,fy,g,muold,muprev,omega,p,p0,px,py,q,s,t,u,v,x,xx
86 C JJV I added a variable 'll' here - it is the 'l' for CASE A
87  INTEGER j,k,kflag,l,ll,m
88 C ..
89 C .. Local Arrays ..
90  REAL fact(10),pp(35)
91 C ..
92 C .. External Functions ..
93  REAL ranf,sexpo,snorm
94  EXTERNAL ranf,sexpo,snorm
95 C ..
96 C .. Intrinsic Functions ..
97  INTRINSIC abs,alog,exp,float,ifix,max0,min0,sign,sqrt
98 C ..
99 C JJV added this for case: mu unchanged
100 C .. Save statement ..
101  SAVE s, d, l, ll, omega, c3, c2, c1, c0, c, m, p, q, p0,
102  + a0, a1, a2, a3, a4, a5, a6, a7, fact, muprev, muold
103 C ..
104 C JJV end addition - I am including vars in Data statements
105 C .. Data statements ..
106 C JJV changed initial values of MUPREV and MUOLD to -1.0E37
107 C JJV if no one calls IGNPOI with MU = -1.0E37 the first time,
108 C JJV the code shouldn't break
109  DATA muprev,muold/-1.0e37,-1.0e37/
110  DATA a0,a1,a2,a3,a4,a5,a6,a7/-.5,.3333333,-.2500068,.2000118,
111  + -.1661269,.1421878,-.1384794,.1250060/
112  DATA fact/1.,1.,2.,6.,24.,120.,720.,5040.,40320.,362880./
113 C ..
114 C .. Executable Statements ..
115
116  IF (mu.EQ.muprev) go to 10
117  IF (mu.LT.10.0) go to 120
118 C
119 C C A S E A. (RECALCULATION OF S,D,LL IF MU HAS CHANGED)
120 C
121 C JJV This is the case where I changed 'l' to 'll'
122 C JJV Here 'll' is set once and used in a comparison once
123
124  muprev = mu
125  s = sqrt(mu)
126  d = 6.0*mu*mu
127 C
128 C THE POISSON PROBABILITIES PK EXCEED THE DISCRETE NORMAL
129 C PROBABILITIES FK WHENEVER K >= M(MU). LL=IFIX(MU-1.1484)
130 C IS AN UPPER BOUND TO M(MU) FOR ALL MU >= 10 .
131 C
132  ll = ifix(mu-1.1484)
133 C
134 C STEP N. NORMAL SAMPLE - SNORM(IR) FOR STANDARD NORMAL DEVIATE
135 C
136  10 g = mu + s*snorm()
137  IF (g.LT.0.0) go to 20
138  ignpoi = ifix(g)
139 C
140 C STEP I. IMMEDIATE ACCEPTANCE IF IGNPOI IS LARGE ENOUGH
141 C
142  IF (ignpoi.GE.ll) RETURN
143 C
144 C STEP S. SQUEEZE ACCEPTANCE - SUNIF(IR) FOR (0,1)-SAMPLE U
145 C
146  fk = float(ignpoi)
147  difmuk = mu - fk
148  u = ranf()
149  IF (d*u.GE.difmuk*difmuk*difmuk) RETURN
150 C
151 C STEP P. PREPARATIONS FOR STEPS Q AND H.
152 C (RECALCULATIONS OF PARAMETERS IF NECESSARY)
153 C .3989423=(2*PI)**(-.5) .416667E-1=1./24. .1428571=1./7.
154 C THE QUANTITIES B1, B2, C3, C2, C1, C0 ARE FOR THE HERMITE
155 C APPROXIMATIONS TO THE DISCRETE NORMAL PROBABILITIES FK.
156 C C=.1069/MU GUARANTEES MAJORIZATION BY THE 'HAT'-FUNCTION.
157 C
158  20 IF (mu.EQ.muold) go to 30
159  muold = mu
160  omega = .3989423/s
161  b1 = .4166667e-1/mu
162  b2 = .3*b1*b1
163  c3 = .1428571*b1*b2
164  c2 = b2 - 15.*c3
165  c1 = b1 - 6.*b2 + 45.*c3
166  c0 = 1. - b1 + 3.*b2 - 15.*c3
167  c = .1069/mu
168  30 IF (g.LT.0.0) go to 50
169 C
170 C 'SUBROUTINE' F IS CALLED (KFLAG=0 FOR CORRECT RETURN)
171 C
172  kflag = 0
173  go to 70
174 C
175 C STEP Q. QUOTIENT ACCEPTANCE (RARE CASE)
176 C
177  40 IF (fy-u*fy.LE.py*exp(px-fx)) RETURN
178 C
179 C STEP E. EXPONENTIAL SAMPLE - SEXPO(IR) FOR STANDARD EXPONENTIAL
180 C DEVIATE E AND SAMPLE T FROM THE LAPLACE 'HAT'
181 C (IF T <= -.6744 THEN PK < FK FOR ALL MU >= 10.)
182 C
183  50 e = sexpo()
184  u = ranf()
185  u = u + u - 1.0
186  t = 1.8 + sign(e,u)
187  IF (t.LE. (-.6744)) go to 50
188  ignpoi = ifix(mu+s*t)
189  fk = float(ignpoi)
190  difmuk = mu - fk
191 C
192 C 'SUBROUTINE' F IS CALLED (KFLAG=1 FOR CORRECT RETURN)
193 C
194  kflag = 1
195  go to 70
196 C
197 C STEP H. HAT ACCEPTANCE (E IS REPEATED ON REJECTION)
198 C
199  60 IF (c*abs(u).GT.py*exp(px+e)-fy*exp(fx+e)) go to 50
200  RETURN
201 C
202 C STEP F. 'SUBROUTINE' F. CALCULATION OF PX,PY,FX,FY.
203 C CASE IGNPOI .LT. 10 USES FACTORIALS FROM TABLE FACT
204 C
205  70 IF (ignpoi.GE.10) go to 80
206  px = -mu
207  py = mu**ignpoi/fact(ignpoi+1)
208  go to 110
209 C
210 C CASE IGNPOI .GE. 10 USES POLYNOMIAL APPROXIMATION
211 C A0-A7 FOR ACCURACY WHEN ADVISABLE
212 C .8333333E-1=1./12. .3989423=(2*PI)**(-.5)
213 C
214  80 del = .8333333e-1/fk
215  del = del - 4.8*del*del*del
216  v = difmuk/fk
217  IF (abs(v).LE.0.25) go to 90
218  px = fk*alog(1.0+v) - difmuk - del
219  go to 100
220
221  90 px = fk*v*v* (((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v+a0) -
222  + del
223  100 py = .3989423/sqrt(fk)
224  110 x = (0.5-difmuk)/s
225  xx = x*x
226  fx = -0.5*xx
227  fy = omega* (((c3*xx+c2)*xx+c1)*xx+c0)
228  IF (kflag.LE.0) go to 40
229  go to 60
230 C
231 C C A S E B. (START NEW TABLE AND CALCULATE P0 IF NECESSARY)
232 C
233 C JJV changed MUPREV assignment from 0.0 to initial value
234  120 muprev = -1.0e37
235  IF (mu.EQ.muold) go to 130
236 C JJV added argument checker here
237  IF (mu.GE.0.0) go to 125
238  WRITE (*,*) 'MU < 0 in IGNPOI - ABORT'
239  WRITE (*,*) 'Value of MU: ',mu
240  CALL xstopx('MU < 0 in IGNPOI - ABORT')
241 C JJV added line label here
242  125 muold = mu
243  m = max0(1,ifix(mu))
244  l = 0
245  p = exp(-mu)
246  q = p
247  p0 = p
248 C
249 C STEP U. UNIFORM SAMPLE FOR INVERSION METHOD
250 C
251  130 u = ranf()
252  ignpoi = 0
253  IF (u.LE.p0) RETURN
254 C
255 C STEP T. TABLE COMPARISON UNTIL THE END PP(L) OF THE
256 C PP-TABLE OF CUMULATIVE POISSON PROBABILITIES
257 C (0.458=PP(9) FOR MU=10)
258 C
259  IF (l.EQ.0) go to 150
260  j = 1
261  IF (u.GT.0.458) j = min0(l,m)
262  DO 140 k = j,l
263  IF (u.LE.pp(k)) go to 180
264  140 CONTINUE
265  IF (l.EQ.35) go to 130
266 C
267 C STEP C. CREATION OF NEW POISSON PROBABILITIES P
268 C AND THEIR CUMULATIVES Q=PP(K)
269 C
270  150 l = l + 1
271  DO 160 k = l,35
272  p = p*mu/float(k)
273  q = q + p
274  pp(k) = q
275  IF (u.LE.q) go to 170
276  160 CONTINUE
277  l = 35
278  go to 130
279
280  170 l = k
281  180 ignpoi = k
282  RETURN
283
284  END