1 DOUBLE PRECISION FUNCTION genbet(aa,bb) 2C********************************************************************** 3C 4C DOUBLE PRECISION FUNCTION GENBET( A, B ) 5C GeNerate BETa random deviate 6C 7C 8C Function 9C 10C 11C Returns a single random deviate from the beta distribution with 12C parameters A and B. The density of the beta is 13C x^(a-1) * (1-x)^(b-1) / B(a,b) for 0 < x < 1 14C 15C 16C Arguments 17C 18C 19C A --> First parameter of the beta distribution 20C DOUBLE PRECISION A 21C JJV (A > 1.0E-37) 22C 23C B --> Second parameter of the beta distribution 24C DOUBLE PRECISION B 25C JJV (B > 1.0E-37) 26C 27C 28C Method 29C 30C 31C R. C. H. Cheng 32C Generating Beta Variates with Nonintegral Shape Parameters 33C Communications of the ACM, 21:317-322 (1978) 34C (Algorithms BB and BC) 35C 36C********************************************************************** 37C .. Parameters .. 38C Close to the largest number that can be exponentiated 39 DOUBLE PRECISION expmax 40C JJV changed this - 89 was too high, and LOG(1.0E38) = 87.49823 41 PARAMETER (expmax=87.49823) 42C Close to the largest representable single precision number 43 DOUBLE PRECISION infnty 44 PARAMETER (infnty=1.0D38) 45C JJV added the parameter minlog 46C Close to the smallest number of which a LOG can be taken. 47 DOUBLE PRECISION minlog 48 PARAMETER (minlog=1.0D-37) 49C .. 50C .. Scalar Arguments .. 51 DOUBLE PRECISION aa,bb 52C .. 53C .. Local Scalars .. 54 DOUBLE PRECISION a,alpha,b,beta,delta,gamma,k1,k2,olda,oldb,r,s, 55 + t,u1,u2,v,w,y,z 56 LOGICAL qsame 57C .. 58C .. External Functions .. 59 DOUBLE PRECISION ranf 60 EXTERNAL ranf 61C .. 62C .. Intrinsic Functions .. 63 INTRINSIC exp,log,max,min,sqrt 64C .. 65C .. Save statement .. 66C JJV added a,b 67 SAVE olda,oldb,alpha,beta,gamma,k1,k2,a,b 68C .. 69C .. Data statements .. 70C JJV changed these to ridiculous values 71 DATA olda,oldb/-1.0E37,-1.0E37/ 72C .. 73C .. Executable Statements .. 74 qsame = (olda.EQ.aa) .AND. (oldb.EQ.bb) 75 IF (qsame) GO TO 20 76C JJV added small minimum for small log problem in calc of W 77 IF (.NOT. (aa.LT.minlog.OR.bb.LT.minlog)) GO TO 10 78 WRITE (*,*) ' AA or BB < ',minlog,' in GENBET - Abort!' 79 WRITE (*,*) ' AA: ',aa,' BB ',bb 80 STOP ' AA or BB too small in GENBET - Abort!' 81 82 10 olda = aa 83 oldb = bb 84 20 IF (.NOT. (min(aa,bb).GT.1.0)) GO TO 100 85 86 87C Alborithm BB 88 89C 90C Initialize 91C 92 IF (qsame) GO TO 30 93 a = min(aa,bb) 94 b = max(aa,bb) 95 alpha = a + b 96 beta = sqrt((alpha-2.0)/ (2.0*a*b-alpha)) 97 gamma = a + 1.0/beta 98 30 CONTINUE 99 40 u1 = ranf() 100C 101C Step 1 102C 103 u2 = ranf() 104 v = beta*log(u1/ (1.0-u1)) 105C JJV altered this 106 IF (v.GT.expmax) GO TO 55 107C JJV added checker to see if a*exp(v) will overflow 108C JJV 50 _was_ w = a*exp(v); also note here a > 1.0 109 50 w = exp(v) 110 IF (w.GT.infnty/a) GO TO 55 111 w = a*w 112 GO TO 60 113 55 w = infnty 114 115 60 z = u1**2*u2 116 r = gamma*v - 1.3862944 117 s = a + r - w 118C 119C Step 2 120C 121 IF ((s+2.609438).GE. (5.0*z)) GO TO 70 122C 123C Step 3 124C 125 t = log(z) 126 IF (s.GT.t) GO TO 70 127C 128C Step 4 129C 130C JJV added checker to see if log(alpha/(b+w)) will 131C JJV overflow. If so, we count the log as -INF, and 132C JJV consequently evaluate conditional as true, i.e. 133C JJV the algorithm rejects the trial and starts over 134C JJV May not need this here since ALPHA > 2.0 135 IF (alpha/(b+w).LT.minlog) GO TO 40 136 137 IF ((r+alpha*log(alpha/ (b+w))).LT.t) GO TO 40 138C 139C Step 5 140C 141 70 IF (.NOT. (aa.EQ.a)) GO TO 80 142 genbet = w/ (b+w) 143 GO TO 90 144 145 80 genbet = b/ (b+w) 146 90 GO TO 230 147 148 149C Algorithm BC 150 151C 152C Initialize 153C 154 100 IF (qsame) GO TO 110 155 a = max(aa,bb) 156 b = min(aa,bb) 157 alpha = a + b 158 beta = 1.0/b 159 delta = 1.0 + a - b 160 k1 = delta* (0.0138889+0.0416667*b)/ (a*beta-0.777778) 161 k2 = 0.25 + (0.5+0.25/delta)*b 162 110 CONTINUE 163 120 u1 = ranf() 164C 165C Step 1 166C 167 u2 = ranf() 168 IF (u1.GE.0.5) GO TO 130 169C 170C Step 2 171C 172 y = u1*u2 173 z = u1*y 174 IF ((0.25*u2+z-y).GE.k1) GO TO 120 175 GO TO 170 176C 177C Step 3 178C 179 130 z = u1**2*u2 180 IF (.NOT. (z.LE.0.25)) GO TO 160 181 v = beta*log(u1/ (1.0-u1)) 182 183C JJV instead of checking v > expmax at top, I will check 184C JJV if a < 1, then check the appropriate values 185 186 IF (a.GT.1.0) GO TO 135 187C JJV A < 1 so it can help out if EXP(V) would overflow 188 IF (v.GT.expmax) GO TO 132 189 w = a*exp(v) 190 GO TO 200 191 132 w = v + log(a) 192 IF (w.GT.expmax) GO TO 140 193 w = exp(w) 194 GO TO 200 195 196C JJV in this case A > 1 197 135 IF (v.GT.expmax) GO TO 140 198 w = exp(v) 199 IF (w.GT.infnty/a) GO TO 140 200 w = a*w 201 GO TO 200 202 140 w = infnty 203 GO TO 200 204 205 160 IF (z.GE.k2) GO TO 120 206C 207C Step 4 208C 209C 210C Step 5 211C 212 170 v = beta*log(u1/ (1.0-u1)) 213 214C JJV same kind of checking as above 215 IF (a.GT.1.0) GO TO 175 216C JJV A < 1 so it can help out if EXP(V) would overflow 217 IF (v.GT.expmax) GO TO 172 218 w = a*exp(v) 219 GO TO 190 220 172 w = v + log(a) 221 IF (w.GT.expmax) GO TO 180 222 w = exp(w) 223 GO TO 190 224 225C JJV in this case A > 1 226 175 IF (v.GT.expmax) GO TO 180 227 w = exp(v) 228 IF (w.GT.infnty/a) GO TO 180 229 w = a*w 230 GO TO 190 231 232 180 w = infnty 233 234C JJV here we also check to see if log overlows; if so, we treat it 235C JJV as -INF, which means condition is true, i.e. restart 236 190 IF (alpha/(b+w).LT.minlog) GO TO 120 237 IF ((alpha* (log(alpha/ (b+w))+v)-1.3862944).LT.log(z)) GO TO 120 238C 239C Step 6 240C 241 200 IF (.NOT. (a.EQ.aa)) GO TO 210 242 genbet = w/ (b+w) 243 GO TO 220 244 245 210 genbet = b/ (b+w) 246 220 CONTINUE 247 230 RETURN 248 249 END 250