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
330
331
332
333
334
335
336
337
338
339
340
|
#include <math.h>
#include "grand.h"
#include "core_math.h"
int C2F(ignbin)(int *n, double *pp)
/*
**********************************************************************
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 ignbin(long n,float pp)
GENerate BINomial random deviate
Function
Generates a single random deviate from a binomial
distribution whose number of trials is N and whose
probability of an event in each trial is P.
Arguments
n --> The number of trials in the binomial distribution
from which a random deviate is to be generated.
p --> The probability of an event in each trial of the
binomial distribution from which a random deviate
is to be generated.
ignbin <-- A random deviate yielding the number of events
from N independent trials, each of which has
a probability of event P.
Method
This is algorithm BTPE from:
Kachitvichyanukul, V. and Schmeiser, B. W.
Binomial Random Variate Generation.
Communications of the ACM, 31, 2
(February, 1988) 216.
**********************************************************************
SUBROUTINE BTPEC(N,PP,ISEED,JX)
BINOMIAL RANDOM VARIATE GENERATOR
MEAN .LT. 30 -- INVERSE CDF
MEAN .GE. 30 -- ALGORITHM BTPE: ACCEPTANCE-REJECTION VIA
FOUR REGION COMPOSITION. THE FOUR REGIONS ARE A TRIANGLE
(SYMMETRIC IN THE CENTER), A PAIR OF PARALLELOGRAMS (ABOVE
THE TRIANGLE), AND EXPONENTIAL LEFT AND RIGHT TAILS.
BTPE REFERS TO BINOMIAL-TRIANGLE-PARALLELOGRAM-EXPONENTIAL.
BTPEC REFERS TO BTPE AND "COMBINED." THUS BTPE IS THE
RESEARCH AND BTPEC IS THE IMPLEMENTATION OF A COMPLETE
USABLE ALGORITHM.
REFERENCE: VORATAS KACHITVICHYANUKUL AND BRUCE SCHMEISER,
"BINOMIAL RANDOM VARIATE GENERATION,"
COMMUNICATIONS OF THE ACM, FORTHCOMING
WRITTEN: SEPTEMBER 1980.
LAST REVISED: MAY 1985, JULY 1987
REQUIRED SUBPROGRAM: RAND() -- A UNIFORM (0,1) RANDOM NUMBER
GENERATOR
ARGUMENTS
N : NUMBER OF BERNOULLI TRIALS (INPUT)
PP : PROBABILITY OF SUCCESS IN EACH TRIAL (INPUT)
ISEED: RANDOM NUMBER SEED (INPUT AND OUTPUT)
JX: RANDOMLY GENERATED OBSERVATION (OUTPUT)
VARIABLES
PSAVE: VALUE OF PP FROM THE LAST CALL TO BTPEC
NSAVE: VALUE OF N FROM THE LAST CALL TO BTPEC
XNP: VALUE OF THE MEAN FROM THE LAST CALL TO BTPEC
P: PROBABILITY USED IN THE GENERATION PHASE OF BTPEC
FFM: TEMPORARY VARIABLE EQUAL TO XNP + P
M: INTEGER VALUE OF THE CURRENT MODE
FM: FLOATING POINT VALUE OF THE CURRENT MODE
XNPQ: TEMPORARY VARIABLE USED IN SETUP AND SQUEEZING STEPS
P1: AREA OF THE TRIANGLE
C: HEIGHT OF THE PARALLELOGRAMS
XM: CENTER OF THE TRIANGLE
XL: LEFT END OF THE TRIANGLE
XR: RIGHT END OF THE TRIANGLE
AL: TEMPORARY VARIABLE
XLL: RATE FOR THE LEFT EXPONENTIAL TAIL
XLR: RATE FOR THE RIGHT EXPONENTIAL TAIL
P2: AREA OF THE PARALLELOGRAMS
P3: AREA OF THE LEFT EXPONENTIAL TAIL
P4: AREA OF THE RIGHT EXPONENTIAL TAIL
U: A U(0,P4) RANDOM VARIATE USED FIRST TO SELECT ONE OF THE
FOUR REGIONS AND THEN CONDITIONALLY TO GENERATE A VALUE
FROM THE REGION
V: A U(0,1) RANDOM NUMBER USED TO GENERATE THE RANDOM VALUE
(REGION 1) OR TRANSFORMED INTO THE VARIATE TO ACCEPT OR
REJECT THE CANDIDATE VALUE
IX: INTEGER CANDIDATE VALUE
X: PRELIMINARY CONTINUOUS CANDIDATE VALUE IN REGION 2 LOGIC
AND A FLOATING POINT IX IN THE ACCEPT/REJECT LOGIC
K: ABSOLUTE VALUE OF (IX-M)
F: THE HEIGHT OF THE SCALED DENSITY FUNCTION USED IN THE
ACCEPT/REJECT DECISION WHEN BOTH M AND IX ARE SMALL
ALSO USED IN THE INVERSE TRANSFORMATION
R: THE RATIO P/Q
G: CONSTANT USED IN CALCULATION OF PROBABILITY
MP: MODE PLUS ONE, THE LOWER INDEX FOR EXPLICIT CALCULATION
OF F WHEN IX IS GREATER THAN M
IX1: CANDIDATE VALUE PLUS ONE, THE LOWER INDEX FOR EXPLICIT
CALCULATION OF F WHEN IX IS LESS THAN M
I: INDEX FOR EXPLICIT CALCULATION OF F FOR BTPE
AMAXP: MAXIMUM ERROR OF THE LOGARITHM OF NORMAL BOUND
YNORM: LOGARITHM OF NORMAL BOUND
ALV: NATURAL LOGARITHM OF THE ACCEPT/REJECT VARIATE V
X1,F1,Z,W,Z2,X2,F2, AND W2 ARE TEMPORARY VARIABLES TO BE
USED IN THE FINAL ACCEPT/REJECT TEST
QN: PROBABILITY OF NO SUCCESS IN N TRIALS
REMARK
IX AND JX COULD LOGICALLY BE THE SAME VARIABLE, WHICH WOULD
SAVE A MEMORY POSITION AND A LINE OF CODE. HOWEVER, SOME
COMPILERS (E.G.,CDC MNF) OPTIMIZE BETTER WHEN THE ARGUMENTS
ARE NOT INVOLVED.
ISEED NEEDS TO BE DOUBLE PRECISION IF THE IMSL ROUTINE
GGUBFS IS USED TO GENERATE UNIFORM RANDOM NUMBER, OTHERWISE
TYPE OF ISEED SHOULD BE DICTATED BY THE UNIFORM GENERATOR
**********************************************************************
*****DETERMINE APPROPRIATE ALGORITHM AND WHETHER SETUP IS NECESSARY
*/
{
static double psave = -1.0E37;
static int nsave = -214748365;
static int ignbin, i, ix, ix1, k, m, mp, T1;
static double al, alv, amaxp, c, f, f1, f2, ffm, fm, g, p, p1, p2, p3, p4, q, qn, r, u, v, w, w2, x, x1,
x2, xl, xll, xlr, xm, xnp, xnpq, xr, ynorm, z, z2;
if (*pp != psave)
{
goto S10;
}
if (*n != nsave)
{
goto S20;
}
if (xnp < 30.0)
{
goto S150;
}
goto S30;
S10:
/*
*****SETUP, PERFORM ONLY WHEN PARAMETERS CHANGE
*/
/* JJV added the argument checker - involved only renaming 10
JJV and 20 to the checkers and adding checkers
JJV Only remaining problem - if called initially with the
JJV initial values of psave and nsave, it will hang
*/
psave = *pp;
p = Min(psave, 1.0 - psave);
q = 1.0 - p;
S20:
xnp = *n * p;
nsave = *n;
if (xnp < 30.0)
{
goto S140;
}
ffm = xnp + p;
m = (int)ffm;
fm = m;
xnpq = xnp * q;
p1 = (int) (2.195 * sqrt(xnpq) - 4.6 * q) + 0.5;
xm = fm + 0.5;
xl = xm - p1;
xr = xm + p1;
c = 0.134 + 20.5 / (15.3 + fm);
al = (ffm - xl) / (ffm - xl * p);
xll = al * (1.0 + 0.5 * al);
al = (xr - ffm) / (xr * q);
xlr = al * (1.0 + 0.5 * al);
p2 = p1 * (1.0 + c + c);
p3 = p2 + c / xll;
p4 = p3 + c / xlr;
S30:
/*
*****GENERATE VARIATE
*/
u = C2F(ranf)() * p4;
v = C2F(ranf)();
/*
TRIANGULAR REGION
*/
if (u > p1)
{
goto S40;
}
ix = (int)(xm - p1 * v + u);
goto S170;
S40:
/*
PARALLELOGRAM REGION
*/
if (u > p2)
{
goto S50;
}
x = xl + (u - p1) / c;
v = v * c + 1.0 - Abs(xm - x) / p1;
if (v > 1.0 || v <= 0.0)
{
goto S30;
}
ix = (int)x;
goto S70;
S50:
/*
LEFT TAIL
*/
if (u > p3)
{
goto S60;
}
ix = (int)(xl + log(v) / xll);
if (ix < 0)
{
goto S30;
}
v *= ((u - p2) * xll);
goto S70;
S60:
/*
RIGHT TAIL
*/
ix = (int)(xr - log(v) / xlr);
if (ix > *n)
{
goto S30;
}
v *= ((u - p3) * xlr);
S70:
/*
*****DETERMINE APPROPRIATE WAY TO PERFORM ACCEPT/REJECT TEST
*/
k = Abs(ix - m);
if (k > 20 && k < xnpq / 2 - 1)
{
goto S130;
}
/*
EXPLICIT EVALUATION
*/
f = 1.0;
r = p / q;
g = (*n + 1) * r;
T1 = m - ix;
if (T1 < 0)
{
goto S80;
}
else if (T1 == 0)
{
goto S120;
}
else
{
goto S100;
}
S80:
mp = m + 1;
for (i = mp; i <= ix; i++)
{
f *= (g / i - r);
}
goto S120;
S100:
ix1 = ix + 1;
for (i = ix1; i <= m; i++)
{
f /= (g / i - r);
}
S120:
if (v <= f)
{
goto S170;
}
goto S30;
S130:
/*
SQUEEZING USING UPPER AND LOWER BOUNDS ON ALOG(F(X))
*/
amaxp = k / xnpq * ((k * (k / 3.0 + 0.625) + 0.1666666666666) / xnpq + 0.5);
ynorm = -(k * k / (2.0 * xnpq));
alv = log(v);
if (alv < ynorm - amaxp)
{
goto S170;
}
if (alv > ynorm + amaxp)
{
goto S30;
}
/*
STIRLING'S FORMULA TO MACHINE ACCURACY FOR
THE FINAL ACCEPTANCE/REJECTION TEST
*/
x1 = ix + 1.0;
f1 = fm + 1.0;
z = *n + 1.0 - fm;
w = *n - ix + 1.0;
z2 = z * z;
x2 = x1 * x1;
f2 = f1 * f1;
w2 = w * w;
if (alv <= xm * log(f1 / x1) + (*n - m + 0.5)*log(z / w) + (ix - m)*log(w * p / (x1 * q)) + (13860.0 -
(462.0 - (132.0 - (99.0 - 140.0 / f2) / f2) / f2) / f2) / f1 / 166320.0 + (13860.0 - (462.0 -
(132.0 - (99.0 - 140.0 / z2) / z2) / z2) / z2) / z / 166320.0 + (13860.0 - (462.0 - (132.0 -
(99.0 - 140.0 / x2) / x2) / x2) / x2) / x1 / 166320.0 + (13860.0 - (462.0 - (132.0 - (99.0
- 140.0 / w2) / w2) / w2) / w2) / w / 166320.0)
{
goto S170;
}
goto S30;
S140:
/*
INVERSE CDF LOGIC FOR MEAN LESS THAN 30
*/
qn = pow((double)q, (double) * n);
r = p / q;
g = r * (*n + 1);
S150:
ix = 0;
f = qn;
u = C2F(ranf)();
S160:
if (u < f)
{
goto S170;
}
if (ix > 110)
{
goto S150;
}
u -= f;
ix += 1;
f *= (g / ix - r);
goto S160;
S170:
if (psave > 0.5)
{
ix = *n - ix;
}
ignbin = ix;
return ignbin;
}
|