GNU Octave  6.2.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*4 FUNCTION ignpoi(mu)
2 C**********************************************************************
3 C
4 C INTEGER*4 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
32 C instead of SUNIF.
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*4 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*4 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, pp, 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  DATA pp/35*0.0/
114 C ..
115 C .. Executable Statements ..
116 
117  IF (mu.EQ.muprev) GO TO 10
118  IF (mu.LT.10.0) GO TO 120
119 C
120 C C A S E A. (RECALCULATION OF S,D,LL IF MU HAS CHANGED)
121 C
122 C JJV This is the case where I changed 'l' to 'll'
123 C JJV Here 'll' is set once and used in a comparison once
124 
125  muprev = mu
126  s = sqrt(mu)
127  d = 6.0*mu*mu
128 C
129 C THE POISSON PROBABILITIES PK EXCEED THE DISCRETE NORMAL
130 C PROBABILITIES FK WHENEVER K >= M(MU). LL=IFIX(MU-1.1484)
131 C IS AN UPPER BOUND TO M(MU) FOR ALL MU >= 10 .
132 C
133  ll = ifix(mu-1.1484)
134 C
135 C STEP N. NORMAL SAMPLE - SNORM(IR) FOR STANDARD NORMAL DEVIATE
136 C
137  10 g = mu + s*snorm()
138  IF (g.LT.0.0) GO TO 20
139  ignpoi = ifix(g)
140 C
141 C STEP I. IMMEDIATE ACCEPTANCE IF IGNPOI IS LARGE ENOUGH
142 C
143  IF (ignpoi.GE.ll) RETURN
144 C
145 C STEP S. SQUEEZE ACCEPTANCE - SUNIF(IR) FOR (0,1)-SAMPLE U
146 C
147  fk = float(ignpoi)
148  difmuk = mu - fk
149  u = ranf()
150  IF (d*u.GE.difmuk*difmuk*difmuk) RETURN
151 C
152 C STEP P. PREPARATIONS FOR STEPS Q AND H.
153 C (RECALCULATIONS OF PARAMETERS IF NECESSARY)
154 C .3989423=(2*PI)**(-.5) .416667E-1=1./24. .1428571=1./7.
155 C THE QUANTITIES B1, B2, C3, C2, C1, C0 ARE FOR THE HERMITE
156 C APPROXIMATIONS TO THE DISCRETE NORMAL PROBABILITIES FK.
157 C C=.1069/MU GUARANTEES MAJORIZATION BY THE 'HAT'-FUNCTION.
158 C
159  20 IF (mu.EQ.muold) GO TO 30
160  muold = mu
161  omega = .3989423/s
162  b1 = .4166667e-1/mu
163  b2 = .3*b1*b1
164  c3 = .1428571*b1*b2
165  c2 = b2 - 15.*c3
166  c1 = b1 - 6.*b2 + 45.*c3
167  c0 = 1. - b1 + 3.*b2 - 15.*c3
168  c = .1069/mu
169  30 IF (g.LT.0.0) GO TO 50
170 C
171 C 'SUBROUTINE' F IS CALLED (KFLAG=0 FOR CORRECT RETURN)
172 C
173  kflag = 0
174  GO TO 70
175 C
176 C STEP Q. QUOTIENT ACCEPTANCE (RARE CASE)
177 C
178  40 IF (fy-u*fy.LE.py*exp(px-fx)) RETURN
179 C
180 C STEP E. EXPONENTIAL SAMPLE - SEXPO(IR) FOR STANDARD EXPONENTIAL
181 C DEVIATE E AND SAMPLE T FROM THE LAPLACE 'HAT'
182 C (IF T <= -.6744 THEN PK < FK FOR ALL MU >= 10.)
183 C
184  50 e = sexpo()
185  u = ranf()
186  u = u + u - 1.0
187  t = 1.8 + sign(e,u)
188  IF (t.LE. (-.6744)) GO TO 50
189  ignpoi = ifix(mu+s*t)
190  fk = float(ignpoi)
191  difmuk = mu - fk
192 C
193 C 'SUBROUTINE' F IS CALLED (KFLAG=1 FOR CORRECT RETURN)
194 C
195  kflag = 1
196  GO TO 70
197 C
198 C STEP H. HAT ACCEPTANCE (E IS REPEATED ON REJECTION)
199 C
200  60 IF (c*abs(u).GT.py*exp(px+e)-fy*exp(fx+e)) GO TO 50
201  RETURN
202 C
203 C STEP F. 'SUBROUTINE' F. CALCULATION OF PX,PY,FX,FY.
204 C CASE IGNPOI .LT. 10 USES FACTORIALS FROM TABLE FACT
205 C
206  70 IF (ignpoi.GE.10) GO TO 80
207  px = -mu
208  py = mu**ignpoi/fact(ignpoi+1)
209  GO TO 110
210 C
211 C CASE IGNPOI .GE. 10 USES POLYNOMIAL APPROXIMATION
212 C A0-A7 FOR ACCURACY WHEN ADVISABLE
213 C .8333333E-1=1./12. .3989423=(2*PI)**(-.5)
214 C
215  80 del = .8333333e-1/fk
216  del = del - 4.8*del*del*del
217  v = difmuk/fk
218  IF (abs(v).LE.0.25) GO TO 90
219  px = fk*alog(1.0+v) - difmuk - del
220  GO TO 100
221 
222  90 px = fk*v*v* (((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v+a0) -
223  + del
224  100 py = .3989423/sqrt(fk)
225  110 x = (0.5-difmuk)/s
226  xx = x*x
227  fx = -0.5*xx
228  fy = omega* (((c3*xx+c2)*xx+c1)*xx+c0)
229  IF (kflag.LE.0) GO TO 40
230  GO TO 60
231 C
232 C C A S E B. (START NEW TABLE AND CALCULATE P0 IF NECESSARY)
233 C
234 C JJV changed MUPREV assignment from 0.0 to initial value
235  120 muprev = -1.0e37
236  IF (mu.EQ.muold) GO TO 130
237 C JJV added argument checker here
238  IF (mu.GE.0.0) GO TO 125
239  WRITE (*,*) 'MU < 0 in IGNPOI - ABORT'
240  WRITE (*,*) 'Value of MU: ',mu
241  CALL xstopx ('MU < 0 in IGNPOI - ABORT')
242 C JJV added line label here
243  125 muold = mu
244  m = max0(1,ifix(mu))
245  l = 0
246  p = exp(-mu)
247  q = p
248  p0 = p
249 C
250 C STEP U. UNIFORM SAMPLE FOR INVERSION METHOD
251 C
252  130 u = ranf()
253  ignpoi = 0
254  IF (u.LE.p0) RETURN
255 C
256 C STEP T. TABLE COMPARISON UNTIL THE END PP(L) OF THE
257 C PP-TABLE OF CUMULATIVE POISSON PROBABILITIES
258 C (0.458=PP(9) FOR MU=10)
259 C
260  IF (l.EQ.0) GO TO 150
261  j = 1
262  IF (u.GT.0.458) j = min0(l,m)
263  DO 140 k = j,l
264  IF (u.LE.pp(k)) GO TO 180
265  140 CONTINUE
266  IF (l.EQ.35) GO TO 130
267 C
268 C STEP C. CREATION OF NEW POISSON PROBABILITIES P
269 C AND THEIR CUMULATIVES Q=PP(K)
270 C
271  150 l = l + 1
272  DO 160 k = l,35
273  p = p*mu/float(k)
274  q = q + p
275  pp(k) = q
276  IF (u.LE.q) GO TO 170
277  160 CONTINUE
278  l = 35
279  GO TO 130
280 
281  170 l = k
282  180 ignpoi = k
283  RETURN
284 
285  END
integer *4 function ignpoi(mu)
Definition: ignpoi.f:2
static T abs(T x)
Definition: pr-output.cc:1678
real function ranf()
Definition: ranf.f:2
real function sexpo()
Definition: sexpo.f:2
real function snorm()
Definition: snorm.f:2