C ALGORITHM 678, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 15, NO. 4, PP. 394-397. SUBROUTINE BTPEC(N,PP,ISEED,JX) C C BINOMIAL RANDOM VARIATE GENERATOR C MEAN .LT. 30 -- INVERSE CDF C MEAN .GE. 30 -- ALGORITHM BTPE: ACCEPTANCE-REJECTION VIA C FOUR REGION COMPOSITION. THE FOUR REGIONS ARE A TRIANGLE C (SYMMETRIC IN THE CENTER), A PAIR OF PARALLELOGRAMS (ABOVE C THE TRIANGLE), AND EXPONENTIAL LEFT AND RIGHT TAILS. C C BTPE REFERS TO BINOMIAL-TRIANGLE-PARALLELOGRAM-EXPONENTIAL. C BTPEC REFERS TO BTPE AND "COMBINED." THUS BTPE IS THE C RESEARCH CONTRIBUTION AND BTPEC IS THE IMPLEMENTATION OF A C COMPLETE USABLE ALGORITHM. C REFERENCE: VORATAS KACHITVICHYANUKUL AND BRUCE SCHMEISER, C "BINOMIAL RANDOM VARIATE GENERATION," C COMMUNICATIONS OF THE ACM, VOLUME 31, NUMBER 2, 216-222, 1988 C C WRITTEN: C BTPE : SEPTEMBER 1980. C BTPEC : MAY 1985. C COMMENTS LAST MODIFIED JUNE 1988 C C REQUIRED SUBPROGRAM: RAND() -- A UNIFORM (0,1) RANDOM NUMBER C GENERATOR C ARGUMENTS C C N : NUMBER OF BERNOULLI TRIALS (INPUT) C PP : PROBABILITY OF SUCCESS IN EACH TRIAL (INPUT) C ISEED: RANDOM NUMBER SEED (INPUT AND OUTPUT) C JX: RANDOMLY GENERATED OBSERVATION (OUTPUT) C C VARIABLES C PSAVE: VALUE OF PP FROM THE LAST CALL TO BTPEC C NSAVE: VALUE OF N FROM THE LAST CALL TO BTPEC C XNP: VALUE OF THE MEAN FROM THE LAST CALL TO BTPEC C C P: PROBABILITY USED IN THE GENERATION PHASE OF BTPEC C FFM: TEMPORARY VARIABLE EQUAL TO XNP + P C M: INTEGER REPRESENTATION OF THE CURRENT MODE C FM: FLOATING POINT REPRESENTATION OF THE CURRENT MODE C XNPQ: TEMPORARY VARIABLE USED IN SETUP AND SQUEEZING STEPS C P1: AREA OF THE TRIANGLE C C: HEIGHT OF THE PARALLELOGRAMS C XM: CENTER OF THE TRIANGLE C XL: LEFT END OF THE TRIANGLE C XR: RIGHT END OF THE TRIANGLE C AL: TEMPORARY VARIABLE C XLL: RATE FOR THE LEFT EXPONENTIAL TAIL C XLR: RATE FOR THE RIGHT EXPONENTIAL TAIL C P2: AREA OF THE TWO PARALLELOGRAMS PLUS P1 C P3: AREA OF THE LEFT EXPONENTIAL TAIL PLUS P2 C P4: AREA OF THE RIGHT EXPONENTIAL TAIL PLUS P3 C U: A U(0,P4) RANDOM VARIATE USED FIRST TO SELECT ONE OF THE C FOUR REGIONS AND THEN CONDITIONALLY TO GENERATE A VALUE C FROM THE REGION C V: A U(0,1) RANDOM NUMBER USED TO GENERATE THE RANDOM VALUE C (REGION 1) OR TRANSFORMED INTO THE VARIATE TO ACCEPT OR C REJECT THE CANDIDATE VALUE C IX: INTEGER CANDIDATE VALUE C X: PRELIMINARY CONTINUOUS CANDIDATE VALUE IN REGION 2 LOGIC C AND A FLOATING POINT IX IN THE ACCEPT/REJECT LOGIC C K: ABSOLUTE VALUE OF (IX-M) C F: THE HEIGHT OF THE SCALED DENSITY FUNCTION USED IN THE C ACCEPT/REJECT DECISION WHEN BOTH M AND IX ARE SMALL C ALSO USED IN THE INVERSE TRANSFORMATION C R: THE RATIO P/Q C G: CONSTANT USED IN CALCULATION OF PROBABILITY C MP: MODE PLUS ONE, THE LOWER INDEX FOR EXPLICIT CALCULATION C OF F WHEN IX IS GREATER THAN M C IX1: CANDIDATE VALUE PLUS ONE, THE LOWER INDEX FOR EXPLICIT C CALCULATION OF F WHEN IX IS LESS THAN M C I: INDEX FOR EXPLICIT CALCULATION OF F FOR BTPE C AMAXP: MAXIMUM ERROR OF THE LOGARITHM OF THE NORMAL BOUND C YNORM: LOGARITHM OF THE NORMAL BOUND C ALV: NATURAL LOGARITHM OF THE ACCEPT/REJECT VARIATE V C X1,F1,Z,W,Z2,X2,F2, AND W2 ARE TEMPORARY VARIABLES C USED IN THE FINAL ACCEPT/REJECT TEST C QN: PROBABILITY OF NO SUCCESS IN N TRIALS C C REMARK C IX AND JX COULD LOGICALLY BE THE SAME VARIABLE, WHICH WOULD C SAVE A MEMORY POSITION AND A LINE OF CODE. HOWEVER, SOME C COMPILERS (E.G.,CDC MNF) OPTIMIZE BETTER WHEN THE ARGUMENTS C ARE NOT INVOLVED. C C ISEED NEEDS TO BE DOUBLE PRECISION IF THE IMSL ROUTINE C GGUBFS IS USED TO GENERATE THE UNIFORM RANDOM NUMBER, OTHERWISE C TYPE OF ISEED IS DICTATED BY THE UNIFORM GENERATOR C DATA PSAVE,NSAVE/-1.,-1/ C C*****DETERMINE APPROPRIATE ALGORITHM AND WHETHER SETUP IS NECESSARY C IF(PP.NE.PSAVE) GO TO 1 IF(N.NE.NSAVE) GO TO 99 IF(XNP-30.) 14,2,2 C C*****SETUP, PERFORM ONLY WHEN PARAMETERS CHANGE C 1 PSAVE = PP P = MIN (PSAVE,1.-PSAVE) Q = 1. - P 99 XNP = N * P NSAVE = N IF (XNP .LT. 30.) GO TO 13 FFM = XNP + P M = 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. + .5 * AL) AL = (XR - FFM) / (XR*Q) XLR = AL * (1. + .5 * AL) P2 = P1 * (1. + C + C) P3 = P2 + C/XLL P4 = P3 + C/XLR C C*****GENERATE VARIATE C 2 U = RAND(ISEED) * P4 V = RAND(ISEED) C C TRIANGULAR REGION C IF (U .GT. P1) GO TO 3 IX = XM - P1 * V + U GO TO 16 C C PARALLELOGRAM REGION C 3 IF (U .GT. P2) GO TO 4 X = XL + (U - P1) / C V = V * C + 1. - ABS(XM - X) / P1 IF (V .GT. 1. .OR. V .LE. 0.) GO TO 2 IX = X GO TO 6 C C LEFT TAIL C 4 IF (U .GT. P3) GO TO 5 IX = XL + LOG(V) / XLL IF (IX .LT. 0) GO TO 2 V = V * (U - P2) * XLL GO TO 6 C C RIGHT TAIL C 5 IX = XR - LOG(V) / XLR IF (IX .GT. N) GO TO 2 V = V * (U - P3) * XLR C C*****DETERMINE APPROPRIATE WAY TO PERFORM ACCEPT/REJECT TEST C 6 K = ABS(IX - M) IF (K .GT. 20 .AND. K .LT. XNPQ / 2 - 1) GO TO 12 C C EXPLICIT EVALUATION C F = 1.0 R = P / Q G = (N + 1) * R IF (M - IX) 7,11,9 7 MP = M + 1 DO 8 I = MP,IX 8 F = F * (G / I - R) GO TO 11 9 IX1 = IX + 1 DO 10 I = IX1,M 10 F = F / (G / I - R) 11 IF (V - F) 16,16,2 C C SQUEEZING USING UPPER AND LOWER BOUNDS ON ALOG(F(X)) C 12 AMAXP = (K/XNPQ) * ((K * (K/3.+.625)+.1666666666666)/XNPQ + .5) YNORM = -K * K / (2. * XNPQ) ALV = LOG(V) IF (ALV .LT. YNORM - AMAXP) GO TO 16 IF (ALV .GT. YNORM + AMAXP) GO TO 2 C C STIRLING'S FORMULA TO MACHINE ACCURACY FOR C THE FINAL ACCEPTANCE/REJECTION TEST C X1 = IX + 1 F1 = FM + 1. Z = N + 1 - FM W = N - IX + 1. Z2 = Z * Z X2 = X1 * X1 F2 = F1 * F1 W2 = W * W IF (ALV - (XM*ALOG(F1/X1) + (N-M+.5) * LOG(Z/W) & + (IX-M) * LOG (W * P / (X1*Q)) & +(13860.-(462.-(132.-(99.-140./F2)/F2)/F2)/F2)/F1/166320. & +(13860.-(462.-(132.-(99.-140./Z2)/Z2)/Z2)/Z2)/Z/166320. & +(13860.-(462.-(132.-(99.-140./X2)/X2)/X2)/X2)/X1/166320. & +(13860.-(462.-(132.-(99.-140./W2)/W2)/W2)/W2)/W/166320.))16,16,2 C C INVERSE CDF LOGIC FOR MEAN LESS THAN 30 C 13 QN = Q ** N R = P / Q G = R * (N + 1) 14 IX = 0 F = QN U = RAND(ISEED) 15 IF (U .LT. F) GO TO 16 IF (IX .GT. 110) GO TO 14 U = U - F IX = IX + 1 F = F * (G/IX - R) GO TO 15 16 IF (PSAVE .GT. 0.5 ) IX = N - IX JX = IX RETURN END * * C C EXAMPLE DRIVER PROGRAM FOR ALGORITHM BTPEC C C NS : NUMBER OF RANDOM VARIATES TO BE GENERATED C N, P: PARAMETERS OF THE BINOMIAL DISTRIBUTION C ISEED: RANDOM NUMBER SEED (TYPE DICTATED BY THE UNIFORM C GENERATOR USED) C JX : BINOMIAL RANDOM VARIABLE GENERATED C X: ADJUSTED VALUE USED TO AVOID OVERFLOW IN THE COMPUTATION C OF THE VARIANCE WHEN NS IS LARGE C SMEAN: SAMPLE MEAN C SVAR: SAMPLE VARIANCE C TMEAN: TRUE MEAN C TVAR: TRUE VARIANCE C SUM: SUM OF BINOMIAL RANDOM VARIATE JX C SUM2: SUM OF SQUARES OF JX C C 10 WRITE (6, 8000) READ (5,*,END=99) NS, N, P, ISEED TMEAN = N * P TVAR = TMEAN * (1. - P) SUM = 0.0 SUM2 = 0.0 DO 15 I = 1, NS CALL BTPEC(N,P,ISEED,JX) X = JX - TMEAN SUM = SUM + X SUM2 = SUM2 + X * X 15 CONTINUE SMEAN = SUM / NS + TMEAN SVAR = SUM2 / NS WRITE (6,9000) NS, N, P, TMEAN, SMEAN, TVAR, SVAR GO TO 10 99 STOP 8000 FORMAT(/5X,' ENTER NS, N, P, AND ISEED '/) 9000 FORMAT(/5X,'SAMPLE SIZE = ',I20,/5X, & 'N = ',I10,5X,' P = ',F20.7,/5X, & 'TRUE MEAN = ',F20.7,/5X, & 'SAMPLE MEAN = ',F20.7,/5X, & 'TRUE VARIANCE = ',F20.7,/5X, & 'SAMPLE VARIANCE = ',E20.7,/5X) END FUNCTION RAND(IX) C C UNIFORM RANDOM NUMBER GENERATOR C BY L. SCHRAGE C 'A MORE PORTABLE FORTRAN RANDOM NUMBER GENERATOR' C ACM TRANSACTIONS ON MATHEMATICAL SOFTWARE, C 5:132-138, 1979. C INTEGER A,P,IX,B15,B16,XHI,XALO,LEFTLO,FHI,K DATA A/16807/,B15/32768/,B16/65536/,P/2147483647/ XHI = IX / B16 XALO = (IX - XHI * B16) * A LEFTLO = XALO / B16 FHI = XHI * A + LEFTLO K = FHI / B15 IX = (((XALO - LEFTLO*B16) - P) + (FHI - K*B15) * B16) + K IF (IX .LT. 0) IX = IX + P RAND = FLOAT(IX) * 4.656612875E-10 RETURN END