#include #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; }