GNU Octave  4.2.1 A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
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 ..
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
OCTAVE_EXPORT octave_value_list or N dimensional array whose elements are all equal to the base of natural logarithms The constant ex \$e satisfies the equation log(e)
int exp(void) const
Definition: DET.h:66
real function genbet(aa, bb)
Definition: genbet.f:1
may be zero for pure relative error test tem the relative tolerance must be greater than or equal to