00001 INTEGER FUNCTION ignbin(n,pp) 00002 C********************************************************************** 00003 C 00004 C INTEGER FUNCTION IGNBIN( N, PP ) 00005 C 00006 C GENerate BINomial random deviate 00007 C 00008 C 00009 C Function 00010 C 00011 C 00012 C Generates a single random deviate from a binomial 00013 C distribution whose number of trials is N and whose 00014 C probability of an event in each trial is P. 00015 C 00016 C 00017 C Arguments 00018 C 00019 C 00020 C N --> The number of trials in the binomial distribution 00021 C from which a random deviate is to be generated. 00022 C INTEGER N 00023 C JJV (N >= 0) 00024 C 00025 C PP --> The probability of an event in each trial of the 00026 C binomial distribution from which a random deviate 00027 C is to be generated. 00028 C REAL PP 00029 C JJV (0.0 <= pp <= 1.0) 00030 C 00031 C IGNBIN <-- A random deviate yielding the number of events 00032 C from N independent trials, each of which has 00033 C a probability of event P. 00034 C INTEGER IGNBIN 00035 C 00036 C 00037 C Note 00038 C 00039 C 00040 C Uses RANF so the value of the seeds, ISEED1 and ISEED2 must be set 00041 C by a call similar to the following 00042 C DUM = RANSET( ISEED1, ISEED2 ) 00043 C 00044 C 00045 C Method 00046 C 00047 C 00048 C This is algorithm BTPE from: 00049 C 00050 C Kachitvichyanukul, V. and Schmeiser, B. W. 00051 C 00052 C Binomial Random Variate Generation. 00053 C Communications of the ACM, 31, 2 00054 C (February, 1988) 216. 00055 C 00056 C********************************************************************** 00057 C SUBROUTINE BTPEC(N,PP,ISEED,JX) 00058 C 00059 C BINOMIAL RANDOM VARIATE GENERATOR 00060 C MEAN .LT. 30 -- INVERSE CDF 00061 C MEAN .GE. 30 -- ALGORITHM BTPE: ACCEPTANCE-REJECTION VIA 00062 C FOUR REGION COMPOSITION. THE FOUR REGIONS ARE A TRIANGLE 00063 C (SYMMETRIC IN THE CENTER), A PAIR OF PARALLELOGRAMS (ABOVE 00064 C THE TRIANGLE), AND EXPONENTIAL LEFT AND RIGHT TAILS. 00065 C 00066 C BTPE REFERS TO BINOMIAL-TRIANGLE-PARALLELOGRAM-EXPONENTIAL. 00067 C BTPEC REFERS TO BTPE AND "COMBINED." THUS BTPE IS THE 00068 C RESEARCH AND BTPEC IS THE IMPLEMENTATION OF A COMPLETE 00069 C USABLE ALGORITHM. 00070 C REFERENCE: VORATAS KACHITVICHYANUKUL AND BRUCE SCHMEISER, 00071 C "BINOMIAL RANDOM VARIATE GENERATION," 00072 C COMMUNICATIONS OF THE ACM, FORTHCOMING 00073 C WRITTEN: SEPTEMBER 1980. 00074 C LAST REVISED: MAY 1985, JULY 1987 00075 C REQUIRED SUBPROGRAM: RAND() -- A UNIFORM (0,1) RANDOM NUMBER 00076 C GENERATOR 00077 C ARGUMENTS 00078 C 00079 C N : NUMBER OF BERNOULLI TRIALS (INPUT) 00080 C PP : PROBABILITY OF SUCCESS IN EACH TRIAL (INPUT) 00081 C ISEED: RANDOM NUMBER SEED (INPUT AND OUTPUT) 00082 C JX: RANDOMLY GENERATED OBSERVATION (OUTPUT) 00083 C 00084 C VARIABLES 00085 C PSAVE: VALUE OF PP FROM THE LAST CALL TO BTPEC 00086 C NSAVE: VALUE OF N FROM THE LAST CALL TO BTPEC 00087 C XNP: VALUE OF THE MEAN FROM THE LAST CALL TO BTPEC 00088 C 00089 C P: PROBABILITY USED IN THE GENERATION PHASE OF BTPEC 00090 C FFM: TEMPORARY VARIABLE EQUAL TO XNP + P 00091 C M: INTEGER VALUE OF THE CURRENT MODE 00092 C FM: FLOATING POINT VALUE OF THE CURRENT MODE 00093 C XNPQ: TEMPORARY VARIABLE USED IN SETUP AND SQUEEZING STEPS 00094 C P1: AREA OF THE TRIANGLE 00095 C C: HEIGHT OF THE PARALLELOGRAMS 00096 C XM: CENTER OF THE TRIANGLE 00097 C XL: LEFT END OF THE TRIANGLE 00098 C XR: RIGHT END OF THE TRIANGLE 00099 C AL: TEMPORARY VARIABLE 00100 C XLL: RATE FOR THE LEFT EXPONENTIAL TAIL 00101 C XLR: RATE FOR THE RIGHT EXPONENTIAL TAIL 00102 C P2: AREA OF THE PARALLELOGRAMS 00103 C P3: AREA OF THE LEFT EXPONENTIAL TAIL 00104 C P4: AREA OF THE RIGHT EXPONENTIAL TAIL 00105 C U: A U(0,P4) RANDOM VARIATE USED FIRST TO SELECT ONE OF THE 00106 C FOUR REGIONS AND THEN CONDITIONALLY TO GENERATE A VALUE 00107 C FROM THE REGION 00108 C V: A U(0,1) RANDOM NUMBER USED TO GENERATE THE RANDOM VALUE 00109 C (REGION 1) OR TRANSFORMED INTO THE VARIATE TO ACCEPT OR 00110 C REJECT THE CANDIDATE VALUE 00111 C IX: INTEGER CANDIDATE VALUE 00112 C X: PRELIMINARY CONTINUOUS CANDIDATE VALUE IN REGION 2 LOGIC 00113 C AND A FLOATING POINT IX IN THE ACCEPT/REJECT LOGIC 00114 C K: ABSOLUTE VALUE OF (IX-M) 00115 C F: THE HEIGHT OF THE SCALED DENSITY FUNCTION USED IN THE 00116 C ACCEPT/REJECT DECISION WHEN BOTH M AND IX ARE SMALL 00117 C ALSO USED IN THE INVERSE TRANSFORMATION 00118 C R: THE RATIO P/Q 00119 C G: CONSTANT USED IN CALCULATION OF PROBABILITY 00120 C MP: MODE PLUS ONE, THE LOWER INDEX FOR EXPLICIT CALCULATION 00121 C OF F WHEN IX IS GREATER THAN M 00122 C IX1: CANDIDATE VALUE PLUS ONE, THE LOWER INDEX FOR EXPLICIT 00123 C CALCULATION OF F WHEN IX IS LESS THAN M 00124 C I: INDEX FOR EXPLICIT CALCULATION OF F FOR BTPE 00125 C AMAXP: MAXIMUM ERROR OF THE LOGARITHM OF NORMAL BOUND 00126 C YNORM: LOGARITHM OF NORMAL BOUND 00127 C ALV: NATURAL LOGARITHM OF THE ACCEPT/REJECT VARIATE V 00128 C 00129 C X1,F1,Z,W,Z2,X2,F2, AND W2 ARE TEMPORARY VARIABLES TO BE 00130 C USED IN THE FINAL ACCEPT/REJECT TEST 00131 C 00132 C QN: PROBABILITY OF NO SUCCESS IN N TRIALS 00133 C 00134 C REMARK 00135 C IX AND JX COULD LOGICALLY BE THE SAME VARIABLE, WHICH WOULD 00136 C SAVE A MEMORY POSITION AND A LINE OF CODE. HOWEVER, SOME 00137 C COMPILERS (E.G.,CDC MNF) OPTIMIZE BETTER WHEN THE ARGUMENTS 00138 C ARE NOT INVOLVED. 00139 C 00140 C ISEED NEEDS TO BE DOUBLE PRECISION IF THE IMSL ROUTINE 00141 C GGUBFS IS USED TO GENERATE UNIFORM RANDOM NUMBER, OTHERWISE 00142 C TYPE OF ISEED SHOULD BE DICTATED BY THE UNIFORM GENERATOR 00143 C 00144 C********************************************************************** 00145 00146 C 00147 C 00148 C 00149 C*****DETERMINE APPROPRIATE ALGORITHM AND WHETHER SETUP IS NECESSARY 00150 C 00151 C .. 00152 C .. Scalar Arguments .. 00153 REAL pp 00154 INTEGER n 00155 C .. 00156 C .. Local Scalars .. 00157 REAL al,alv,amaxp,c,f,f1,f2,ffm,fm,g,p,p1,p2,p3,p4,psave,q,qn,r,u, 00158 + v,w,w2,x,x1,x2,xl,xll,xlr,xm,xnp,xnpq,xr,ynorm,z,z2 00159 INTEGER i,ix,ix1,k,m,mp,nsave 00160 C .. 00161 C .. External Functions .. 00162 REAL ranf 00163 EXTERNAL ranf 00164 C .. 00165 C .. Intrinsic Functions .. 00166 INTRINSIC abs,alog,amin1,iabs,int,sqrt 00167 C JJV .. 00168 C JJV .. Save statement .. 00169 SAVE p,q,m,fm,xnp,xnpq,p1,xm,xl,xr,c,xll,xlr,p2,p3,p4,qn,r,g, 00170 + psave,nsave 00171 C JJV I am including the variables in data statements 00172 C .. 00173 C .. Data statements .. 00174 C JJV made these ridiculous starting values - the hope is that 00175 C JJV no one will call this the first time with them as args 00176 DATA psave,nsave/-1.0E37,-214748365/ 00177 C .. 00178 C .. Executable Statements .. 00179 IF (pp.NE.psave) GO TO 10 00180 IF (n.NE.nsave) GO TO 20 00181 IF (xnp-30.0.LT.0.0) GO TO 150 00182 GO TO 30 00183 C 00184 C*****SETUP, PERFORM ONLY WHEN PARAMETERS CHANGE 00185 C 00186 00187 C JJV added the argument checker - involved only renaming 10 00188 C JJV and 20 to the checkers and adding checkers 00189 C JJV Only remaining problem - if called initially with the 00190 C JJV initial values of psave and nsave, it will hang 00191 10 IF (pp.LT.0.0) CALL XSTOPX ('PP < 0.0 in IGNBIN - ABORT!') 00192 IF (pp.GT.1.0) CALL XSTOPX ('PP > 1.0 in IGNBIN - ABORT!') 00193 psave = pp 00194 p = amin1(psave,1.-psave) 00195 q = 1. - p 00196 20 IF (n.LT.0) CALL XSTOPX ('N < 0 in IGNBIN - ABORT!') 00197 xnp = n*p 00198 nsave = n 00199 IF (xnp.LT.30.) GO TO 140 00200 ffm = xnp + p 00201 m = ffm 00202 fm = m 00203 xnpq = xnp*q 00204 p1 = int(2.195*sqrt(xnpq)-4.6*q) + 0.5 00205 xm = fm + 0.5 00206 xl = xm - p1 00207 xr = xm + p1 00208 c = 0.134 + 20.5/ (15.3+fm) 00209 al = (ffm-xl)/ (ffm-xl*p) 00210 xll = al* (1.+.5*al) 00211 al = (xr-ffm)/ (xr*q) 00212 xlr = al* (1.+.5*al) 00213 p2 = p1* (1.+c+c) 00214 p3 = p2 + c/xll 00215 p4 = p3 + c/xlr 00216 C WRITE(6,100) N,P,P1,P2,P3,P4,XL,XR,XM,FM 00217 C 100 FORMAT(I15,4F18.7/5F18.7) 00218 C 00219 C*****GENERATE VARIATE 00220 C 00221 30 u = ranf()*p4 00222 v = ranf() 00223 C 00224 C TRIANGULAR REGION 00225 C 00226 IF (u.GT.p1) GO TO 40 00227 ix = xm - p1*v + u 00228 GO TO 170 00229 C 00230 C PARALLELOGRAM REGION 00231 C 00232 40 IF (u.GT.p2) GO TO 50 00233 x = xl + (u-p1)/c 00234 v = v*c + 1. - abs(xm-x)/p1 00235 IF (v.GT.1. .OR. v.LE.0.) GO TO 30 00236 ix = x 00237 GO TO 70 00238 C 00239 C LEFT TAIL 00240 C 00241 50 IF (u.GT.p3) GO TO 60 00242 ix = xl + alog(v)/xll 00243 IF (ix.LT.0) GO TO 30 00244 v = v* (u-p2)*xll 00245 GO TO 70 00246 C 00247 C RIGHT TAIL 00248 C 00249 60 ix = xr - alog(v)/xlr 00250 IF (ix.GT.n) GO TO 30 00251 v = v* (u-p3)*xlr 00252 C 00253 C*****DETERMINE APPROPRIATE WAY TO PERFORM ACCEPT/REJECT TEST 00254 C 00255 70 k = iabs(ix-m) 00256 IF (k.GT.20 .AND. k.LT.xnpq/2-1) GO TO 130 00257 C 00258 C EXPLICIT EVALUATION 00259 C 00260 f = 1.0 00261 r = p/q 00262 g = (n+1)*r 00263 IF (m-ix.LT.0) GO TO 80 00264 IF (m-ix.EQ.0) GO TO 120 00265 GO TO 100 00266 80 mp = m + 1 00267 DO 90 i = mp,ix 00268 f = f* (g/i-r) 00269 90 CONTINUE 00270 GO TO 120 00271 00272 100 ix1 = ix + 1 00273 DO 110 i = ix1,m 00274 f = f/ (g/i-r) 00275 110 CONTINUE 00276 120 IF (v-f.LE.0) GO TO 170 00277 GO TO 30 00278 C 00279 C SQUEEZING USING UPPER AND LOWER BOUNDS ON ALOG(F(X)) 00280 C 00281 130 amaxp = (k/xnpq)* ((k* (k/3.+.625)+.1666666666666)/xnpq+.5) 00282 ynorm = -k*k/ (2.*xnpq) 00283 alv = alog(v) 00284 IF (alv.LT.ynorm-amaxp) GO TO 170 00285 IF (alv.GT.ynorm+amaxp) GO TO 30 00286 C 00287 C STIRLING'S FORMULA TO MACHINE ACCURACY FOR 00288 C THE FINAL ACCEPTANCE/REJECTION TEST 00289 C 00290 x1 = ix + 1 00291 f1 = fm + 1. 00292 z = n + 1 - fm 00293 w = n - ix + 1. 00294 z2 = z*z 00295 x2 = x1*x1 00296 f2 = f1*f1 00297 w2 = w*w 00298 IF (alv- (xm*alog(f1/x1)+ (n-m+.5)*alog(z/w)+ (ix- 00299 + m)*alog(w*p/ (x1*q))+ (13860.- (462.- (132.- (99.- 00300 + 140./f2)/f2)/f2)/f2)/f1/166320.+ (13860.- (462.- (132.- (99.- 00301 + 140./z2)/z2)/z2)/z2)/z/166320.+ (13860.- (462.- (132.- (99.- 00302 + 140./x2)/x2)/x2)/x2)/x1/166320.+ (13860.- (462.- (132.- (99.- 00303 + 140./w2)/w2)/w2)/w2)/w/166320.) .LE. 0.) GO TO 170 00304 GO TO 30 00305 C 00306 C INVERSE CDF LOGIC FOR MEAN LESS THAN 30 00307 C 00308 140 qn = q**n 00309 r = p/q 00310 g = r* (n+1) 00311 150 ix = 0 00312 f = qn 00313 u = ranf() 00314 160 IF (u.LT.f) GO TO 170 00315 IF (ix.GT.110) GO TO 150 00316 u = u - f 00317 ix = ix + 1 00318 f = f* (g/ix-r) 00319 GO TO 160 00320 00321 170 IF (psave.GT.0.5) ix = n - ix 00322 ignbin = ix 00323 RETURN 00324 00325 END
1.7.1
