GNU Octave  4.0.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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
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 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, 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
int exp(void) const
Definition: DET.h:64
T abs(T x)
Definition: pr-output.cc:3062
octave_value sqrt(void) const
Definition: ov.h:1200