summaryrefslogtreecommitdiff
path: root/modules/randlib/src/c/ignbin.c
blob: f86c47576458b9bbe996e589326ed2eb6a96c8c0 (plain)
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;
}