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