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
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
|
#include <math.h>
#include "grand.h"
#include "core_math.h"
double fsign (double x, double y)
{
if (y >= 0.0)
{
return Abs(x);
}
else
{
return -Abs(x);
}
}
int C2F(ignpoi)(double *mu)
/*
**********************************************************************
This source code was taken in the project "freemat"(BSD license)
This source code was modified by Gaüzère Sabine according to the
modifications done by JJV
long ignpoi(float mu)
GENerate POIsson random deviate
Function
Generates a single random deviate from a Poisson
distribution with mean AV.
Arguments
av --> The mean of the Poisson distribution from which
a random deviate is to be generated.
genexp <-- The random deviate.
Method
Renames KPOIS from TOMS as slightly modified by BWB to use RANF
instead of SUNIF.
For details see:
Ahrens, J.H. and Dieter, U.
Computer Generation of Poisson Deviates
From Modified Normal Distributions.
ACM Trans. Math. Software, 8, 2
(June 1982),163-179
**********************************************************************
**********************************************************************
P O I S S O N DISTRIBUTION
**********************************************************************
**********************************************************************
FOR DETAILS SEE:
AHRENS, J.H. AND DIETER, U.
COMPUTER GENERATION OF POISSON DEVIATES
FROM MODIFIED NORMAL DISTRIBUTIONS.
ACM TRANS. MATH. SOFTWARE, 8,2 (JUNE 1982), 163 - 179.
(SLIGHTLY MODIFIED VERSION OF THE PROGRAM IN THE ABOVE ARTICLE)
**********************************************************************
INTEGER FUNCTION IGNPOI(IR,MU)
INPUT: IR=CURRENT STATE OF BASIC RANDOM NUMBER GENERATOR
MU=MEAN MU OF THE POISSON DISTRIBUTION
OUTPUT: IGNPOI=SAMPLE FROM THE POISSON-(MU)-DISTRIBUTION
MUPREV=PREVIOUS MU, MUOLD=MU AT LAST EXECUTION OF STEP P OR B.
TABLES: COEFFICIENTS A0-A7 FOR STEP F. FACTORIALS FACT
COEFFICIENTS A(K) - FOR PX = FK*V*V*SUM(A(K)*V**K)-DEL
SEPARATION OF CASES A AND B
*/
{
//extern float sign( float num, float sign );
static double a0 = -0.5;
static double a1 = 0.3333333;
static double a2 = -0.2500068;
static double a3 = 0.2000118;
static double a4 = -0.1661269;
static double a5 = 0.1421878;
static double a6 = -0.1384794;
static double a7 = 0.125006;
static double muold = -1.0E37;
static double muprev = -1.0E37;
static double fact[10] =
{
1.0, 1.0, 2.0, 6.0, 24.0, 120.0, 720.0, 5040.0, 40320.0, 362880.0
};
static int ignpoi, j, k, kflag, l, m, ll;
// JJV I added a variable 'll' here - it is the 'l' for CASE A
static double b1, b2, c, c0, c1, c2, c3, d, del, difmuk, e, fk, fx, fy, g, omega, p, p0, px, py, q, s,
t, u, v, x, xx, pp[35];
if (*mu == muprev)
{
goto S10;
}
if (*mu < 10.0)
{
goto S120;
}
/*
C A S E A. (RECALCULATION OF S,D,L IF MU HAS CHANGED)
*/
// JJV This is the case where I changed 'l' to 'll'
// JJV Here 'll' is set once and used in a comparison once
muprev = *mu;
s = sqrt(*mu);
d = 6.0 * *mu * *mu;
/*
THE POISSON PROBABILITIES PK EXCEED THE DISCRETE NORMAL
PROBABILITIES FK WHENEVER K >= M(MU). L=IFIX(MU-1.1484)
IS AN UPPER BOUND TO M(MU) FOR ALL MU >= 10 .
*/
ll = (int) (*mu - 1.1484);
S10:
/*
STEP N. NORMAL SAMPLE - SNORM(IR) FOR STANDARD NORMAL DEVIATE
*/
g = *mu + s * C2F(snorm)();
if (g < 0.0)
{
goto S20;
}
ignpoi = (int) (g);
/*
STEP I. IMMEDIATE ACCEPTANCE IF IGNPOI IS LARGE ENOUGH
*/
if (ignpoi >= ll)
{
return ignpoi;
}
/*
STEP S. SQUEEZE ACCEPTANCE - SUNIF(IR) FOR (0,1)-SAMPLE U
*/
fk = (double)ignpoi;
difmuk = *mu - fk;
u = C2F(ranf)();
if (d*u >= difmuk * difmuk * difmuk)
{
return ignpoi;
}
S20:
/*
STEP P. PREPARATIONS FOR STEPS Q AND H.
(RECALCULATIONS OF PARAMETERS IF NECESSARY)
.3989423=(2*PI)**(-.5) .416667E-1=1./24. .1428571=1./7.
THE QUANTITIES B1, B2, C3, C2, C1, C0 ARE FOR THE HERMITE
APPROXIMATIONS TO THE DISCRETE NORMAL PROBABILITIES FK.
C=.1069/MU GUARANTEES MAJORIZATION BY THE 'HAT'-FUNCTION.
*/
if (*mu == muold)
{
goto S30;
}
muold = *mu;
omega = 0.3989423 / s;
b1 = 4.166667E-2 / *mu;
b2 = 0.3 * b1 * b1;
c3 = 0.1428571 * b1 * b2;
c2 = b2 - 15.0 * c3;
c1 = b1 - 6.0 * b2 + 45.0 * c3;
c0 = 1.0 - b1 + 3.0 * b2 - 15.0 * c3;
c = 0.1069 / *mu;
S30:
if (g < 0.0)
{
goto S50;
}
/*
'SUBROUTINE' F IS CALLED (KFLAG=0 FOR CORRECT RETURN)
*/
kflag = 0;
goto S70;
S40:
/*
STEP Q. QUOTIENT ACCEPTANCE (RARE CASE)
*/
if (fy - u*fy <= py * exp(px - fx))
{
return ignpoi;
}
S50:
/*
STEP E. EXPONENTIAL SAMPLE - SEXPO(IR) FOR STANDARD EXPONENTIAL
DEVIATE E AND SAMPLE T FROM THE LAPLACE 'HAT'
(IF T <= -.6744 THEN PK < FK FOR ALL MU >= 10.)
*/
e = C2F(sexpo)();
u = C2F(ranf)();
u += (u - 1.0);
t = 1.8 + fsign(e, u);
if (t <= -0.6744)
{
goto S50;
}
ignpoi = (int) (*mu + s * t);
fk = (double)ignpoi;
difmuk = *mu - fk;
/*
'SUBROUTINE' F IS CALLED (KFLAG=1 FOR CORRECT RETURN)
*/
kflag = 1;
goto S70;
S60:
/*
STEP H. HAT ACCEPTANCE (E IS REPEATED ON REJECTION)
*/
if (c * Abs(u) > py * exp(px + e) - fy * exp(fx + e))
{
goto S50;
}
return ignpoi;
S70:
/*
STEP F. 'SUBROUTINE' F. CALCULATION OF PX,PY,FX,FY.
CASE IGNPOI .LT. 10 USES FACTORIALS FROM TABLE FACT
*/
if (ignpoi >= 10)
{
goto S80;
}
px = - *mu;
py = pow((double) * mu, (double)ignpoi) / *(fact + ignpoi);
goto S110;
S80:
/*
CASE IGNPOI .GE. 10 USES POLYNOMIAL APPROXIMATION
A0-A7 FOR ACCURACY WHEN ADVISABLE
.8333333E-1=1./12. .3989423=(2*PI)**(-.5)
*/
del = 8.333333E-2 / fk;
del -= (4.8 * del * del * del);
v = difmuk / fk;
if (Abs(v) <= 0.25)
{
goto S90;
}
px = fk * log(1.0 + v) - difmuk - del;
goto S100;
S90:
px = fk * v * v * (((((((a7 * v + a6) * v + a5) * v + a4) * v + a3) * v + a2) * v + a1) * v + a0) - del;
S100:
py = 0.3989423 / sqrt(fk);
S110:
x = (0.5 - difmuk) / s;
xx = x * x;
fx = -0.5 * xx;
fy = omega * (((c3 * xx + c2) * xx + c1) * xx + c0);
if (kflag <= 0)
{
goto S40;
}
goto S60;
S120:
/*
C A S E B. (START NEW TABLE AND CALCULATE P0 IF NECESSARY)
*/
// JJV changed MUPREV assignment from 0.0 to initial value
muprev = -1.0E37;
// Jpc 1999: the next lines seams to produce a bug
// and I finaly commented them out
// IF (mu.EQ.muold) GO TO 130
// JJV added argument checker here
// JJV added line label here
// 125 muold = mu
// if(mu == muold) goto S130;
// muold = mu;
m = Max(1L, (int) (*mu));
l = 0;
p = exp(- *mu);
q = p0 = p;
S130:
/*
STEP U. UNIFORM SAMPLE FOR INVERSION METHOD
*/
u = C2F(ranf)();
ignpoi = 0;
if (u <= p0)
{
return ignpoi;
}
/*
STEP T. TABLE COMPARISON UNTIL THE END PP(L) OF THE
PP-TABLE OF CUMULATIVE POISSON PROBABILITIES
(0.458=PP(9) FOR MU=10)
*/
if (l == 0)
{
goto S150;
}
j = 1;
if (u > 0.458)
{
j = Min(l, m);
}
for (k = j; k <= l; k++)
{
if (u <= *(pp + k - 1))
{
goto S180;
}
}
if (l == 35)
{
goto S130;
}
S150:
/*
STEP C. CREATION OF NEW POISSON PROBABILITIES P
AND THEIR CUMULATIVES Q=PP(K)
*/
l += 1;
for (k = l; k <= 35; k++)
{
p = p* *mu / (double)k;
q += p;
*(pp + k - 1) = q;
if (u <= q)
{
goto S170;
}
}
l = 35;
goto S130;
S170:
l = k;
S180:
ignpoi = k;
return ignpoi;
}
|