GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
genbet.f
Go to the documentation of this file.
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
21 C JJV (A > 1.0E-37)
22 C
23 C B --> Second parameter of the beta distribution
24 C REAL B
25 C JJV (B > 1.0E-37)
26 C
27 C
28 C Method
29 C
30 C
31 C R. C. H. Cheng
32 C Generating Beta Variates with Nonintegral Shape Parameters
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
40 C JJV changed this - 89 was too high, and LOG(1.0E38) = 87.49823
41  parameter(expmax=87.49823)
42 C Close to the largest representable single precision number
43  REAL infnty
44  parameter(infnty=1.0e38)
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)
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 ..
66 C JJV added a,b
67  SAVE olda,oldb,alpha,beta,gamma,k1,k2,a,b
68 C ..
69 C .. Data statements ..
70 C JJV changed these to ridiculous values
71  DATA olda,oldb/-1.0e37,-1.0e37/
72 C ..
73 C .. Executable Statements ..
74  qsame = (olda.EQ.aa) .AND. (oldb.EQ.bb)
75  IF (qsame) go to 20
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!'
79  WRITE (*,*) ' AA: ',aa,' BB ',bb
80  CALL xstopx(' 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 
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))
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
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
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
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 
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))
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
185 
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
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))
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)
223  go to 190
224 
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
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