Mercurial > octave-nkf
annotate libcruft/ranlib/genbet.f @ 10103:0e71ead7359d
Use CALL XSTOPX instead of STOP in Fortran ranlib routines
New call allows Octave's error handler to intercept otherwise
fatal errors in Fortran code
| author | Rik <rdrider0-list@yahoo.com> |
|---|---|
| date | Wed, 13 Jan 2010 22:22:46 -0800 |
| parents | df7c57a6639d |
| children |
| rev | line source |
|---|---|
| 2329 | 1 REAL FUNCTION genbet(aa,bb) |
| 2 C********************************************************************** | |
| 3 C | |
| 4 C REAL FUNCTION GENBET( A, B ) | |
| 5 C GeNerate BETa random deviate | |
| 6 C | |
| 7 C | |
| 8 C Function | |
| 9 C | |
| 10 C | |
| 11 C Returns a single random deviate from the beta distribution with | |
| 12 C parameters A and B. The density of the beta is | |
| 13 C x^(a-1) * (1-x)^(b-1) / B(a,b) for 0 < x < 1 | |
| 14 C | |
| 15 C | |
| 16 C Arguments | |
| 17 C | |
| 18 C | |
| 19 C A --> First parameter of the beta distribution | |
| 20 C REAL A | |
| 3188 | 21 C JJV (A > 1.0E-37) |
| 2329 | 22 C |
| 23 C B --> Second parameter of the beta distribution | |
| 24 C REAL B | |
| 3188 | 25 C JJV (B > 1.0E-37) |
| 2329 | 26 C |
| 27 C | |
| 28 C Method | |
| 29 C | |
| 30 C | |
| 31 C R. C. H. Cheng | |
| 3188 | 32 C Generating Beta Variates with Nonintegral Shape Parameters |
| 2329 | 33 C Communications of the ACM, 21:317-322 (1978) |
| 34 C (Algorithms BB and BC) | |
| 35 C | |
| 36 C********************************************************************** | |
| 37 C .. Parameters .. | |
| 38 C Close to the largest number that can be exponentiated | |
| 39 REAL expmax | |
| 3188 | 40 C JJV changed this - 89 was too high, and LOG(1.0E38) = 87.49823 |
| 41 PARAMETER (expmax=87.49823) | |
| 2329 | 42 C Close to the largest representable single precision number |
| 43 REAL infnty | |
| 44 PARAMETER (infnty=1.0E38) | |
| 3188 | 45 C JJV added the parameter minlog |
| 46 C Close to the smallest number of which a LOG can be taken. | |
| 47 REAL minlog | |
| 48 PARAMETER (minlog=1.0E-37) | |
| 2329 | 49 C .. |
| 50 C .. Scalar Arguments .. | |
| 51 REAL aa,bb | |
| 52 C .. | |
| 53 C .. Local Scalars .. | |
| 54 REAL a,alpha,b,beta,delta,gamma,k1,k2,olda,oldb,r,s,t,u1,u2,v,w,y, | |
| 55 + z | |
| 56 LOGICAL qsame | |
| 57 C .. | |
| 58 C .. External Functions .. | |
| 59 REAL ranf | |
| 60 EXTERNAL ranf | |
| 61 C .. | |
| 62 C .. Intrinsic Functions .. | |
| 63 INTRINSIC exp,log,max,min,sqrt | |
| 64 C .. | |
| 65 C .. Save statement .. | |
| 3188 | 66 C JJV added a,b |
| 67 SAVE olda,oldb,alpha,beta,gamma,k1,k2,a,b | |
| 2329 | 68 C .. |
| 69 C .. Data statements .. | |
| 3188 | 70 C JJV changed these to ridiculous values |
| 71 DATA olda,oldb/-1.0E37,-1.0E37/ | |
| 2329 | 72 C .. |
| 73 C .. Executable Statements .. | |
| 74 qsame = (olda.EQ.aa) .AND. (oldb.EQ.bb) | |
| 75 IF (qsame) GO TO 20 | |
| 3188 | 76 C 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!' | |
| 2329 | 79 WRITE (*,*) ' AA: ',aa,' BB ',bb |
|
10103
0e71ead7359d
Use CALL XSTOPX instead of STOP in Fortran ranlib routines
Rik <rdrider0-list@yahoo.com>
parents:
3188
diff
changeset
|
80 CALL XSTOPX (' AA or BB too small in GENBET - Abort!') |
| 2329 | 81 |
| 82 10 olda = aa | |
| 83 oldb = bb | |
| 84 20 IF (.NOT. (min(aa,bb).GT.1.0)) GO TO 100 | |
| 85 | |
| 86 | |
| 87 C Alborithm BB | |
| 88 | |
| 89 C | |
| 90 C Initialize | |
| 91 C | |
| 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() | |
| 100 C | |
| 101 C Step 1 | |
| 102 C | |
| 103 u2 = ranf() | |
| 104 v = beta*log(u1/ (1.0-u1)) | |
| 3188 | 105 C JJV altered this |
| 106 IF (v.GT.expmax) GO TO 55 | |
| 107 C JJV added checker to see if a*exp(v) will overflow | |
| 108 C 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 | |
| 2329 | 112 GO TO 60 |
| 3188 | 113 55 w = infnty |
| 2329 | 114 |
| 115 60 z = u1**2*u2 | |
| 116 r = gamma*v - 1.3862944 | |
| 117 s = a + r - w | |
| 118 C | |
| 119 C Step 2 | |
| 120 C | |
| 121 IF ((s+2.609438).GE. (5.0*z)) GO TO 70 | |
| 122 C | |
| 123 C Step 3 | |
| 124 C | |
| 125 t = log(z) | |
| 126 IF (s.GT.t) GO TO 70 | |
| 127 C | |
| 128 C Step 4 | |
| 129 C | |
| 3188 | 130 C JJV added checker to see if log(alpha/(b+w)) will |
| 131 C JJV overflow. If so, we count the log as -INF, and | |
| 132 C JJV consequently evaluate conditional as true, i.e. | |
| 133 C JJV the algorithm rejects the trial and starts over | |
| 134 C JJV May not need this here since ALPHA > 2.0 | |
| 135 IF (alpha/(b+w).LT.minlog) GO TO 40 | |
| 136 | |
| 2329 | 137 IF ((r+alpha*log(alpha/ (b+w))).LT.t) GO TO 40 |
| 138 C | |
| 139 C Step 5 | |
| 140 C | |
| 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 | |
| 149 C Algorithm BC | |
| 150 | |
| 151 C | |
| 152 C Initialize | |
| 153 C | |
| 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() | |
| 164 C | |
| 165 C Step 1 | |
| 166 C | |
| 167 u2 = ranf() | |
| 168 IF (u1.GE.0.5) GO TO 130 | |
| 169 C | |
| 170 C Step 2 | |
| 171 C | |
| 172 y = u1*u2 | |
| 173 z = u1*y | |
| 174 IF ((0.25*u2+z-y).GE.k1) GO TO 120 | |
| 175 GO TO 170 | |
| 176 C | |
| 177 C Step 3 | |
| 178 C | |
| 179 130 z = u1**2*u2 | |
| 180 IF (.NOT. (z.LE.0.25)) GO TO 160 | |
| 181 v = beta*log(u1/ (1.0-u1)) | |
| 3188 | 182 |
| 183 C JJV instead of checking v > expmax at top, I will check | |
| 184 C JJV if a < 1, then check the appropriate values | |
| 2329 | 185 |
| 3188 | 186 IF (a.GT.1.0) GO TO 135 |
| 187 C 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 | |
| 196 C 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 | |
| 2329 | 204 |
| 205 160 IF (z.GE.k2) GO TO 120 | |
| 206 C | |
| 207 C Step 4 | |
| 208 C | |
| 209 C | |
| 210 C Step 5 | |
| 211 C | |
| 212 170 v = beta*log(u1/ (1.0-u1)) | |
| 3188 | 213 |
| 214 C JJV same kind of checking as above | |
| 215 IF (a.GT.1.0) GO TO 175 | |
| 216 C 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) | |
| 2329 | 223 GO TO 190 |
| 224 | |
| 3188 | 225 C 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 | |
| 234 C JJV here we also check to see if log overlows; if so, we treat it | |
| 235 C 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 | |
| 2329 | 238 C |
| 239 C Step 6 | |
| 240 C | |
| 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 |