blitz Version 1.0.2
Loading...
Searching...
No Matches
beta.h
Go to the documentation of this file.
1// -*- C++ -*-
2// $Id$
3
4/*
5 * Generate Beta random deviate
6 *
7 * Returns a single random deviate from the beta distribution with
8 * parameters A and B. The density of the beta is
9 * x^(a-1) * (1-x)^(b-1) / B(a,b) for 0 < x < 1
10 *
11 * The mean is a/(a+b).
12 * The variance is ab/((a+b)^2(a+b+1))
13 * The rth moment is (a+r-1)^(r)/(a+b+r-1)^(r)
14 * where a^(r) means a * (a-1) * (a-2) * ... * (a-r+1)
15 *
16 * Method
17 * R. C. H. Cheng
18 * Generating Beta Variates with Nonintegral Shape Parameters
19 * Communications of the ACM, 21:317-322 (1978)
20 * (Algorithms BB and BC)
21 * http://www.acm.org/pubs/citations/journals/toms/1991-17-1/p98-l_ecuyer/
22 *
23 * This class has been adapted from RANDLIB.C 1.3, by
24 * Barry W. Brown, James Lovato, Kathy Russell, and John Venier.
25 *
26 * Adapter's note (TV): This code has gone from Pascal to Fortran to C.
27 * As a result it is a bit of a mess. Note also that the algorithms were
28 * originally designed for 32-bit float, and so some of the constants
29 * below have inadequate precision. This will not cause problems for
30 * casual use, but if you are generating millions of beta variates and
31 * rely on some convergence property, you may have want to worry
32 * about this.
33 *
34 * NEEDS_WORK: dig out the original paper and determine these constants
35 * to precision adequate for 128-bit float.
36 * NEEDS_WORK: turn this into structured code.
37 */
38
39#ifndef BZ_RANDOM_BETA
40#define BZ_RANDOM_BETA
41
42#ifndef BZ_RANDOM_UNIFORM
43 #include <random/uniform.h>
44#endif
45
46#ifndef BZ_NUMINQUIRE_H
47 #include <blitz/numinquire.h>
48#endif
49
50namespace ranlib {
51
52template<typename T = double, typename IRNG = defaultIRNG,
53 typename stateTag = defaultState>
54class Beta : public UniformOpen<T,IRNG,stateTag>
55{
56public:
57 typedef T T_numtype;
58
59 Beta(T a, T b)
60 {
61 setParameters(a, b);
62 }
63
64 Beta(T a, T b, unsigned int i) : UniformOpen<T, IRNG, stateTag>(i)
65 {
66 setParameters(a, b);
67 }
68
69 T random();
70
71 void setParameters(T a, T b)
72 {
73 aa = a;
74 bb = b;
75 infnty = 0.3 * blitz::huge(T());
76 minlog = 0.085 * blitz::tiny(T());
77 expmax = log(infnty);
78 }
79
80protected:
81 T ranf()
82 {
84 }
85
86 T aa, bb;
88};
89
90template<typename T, typename IRNG, typename stateTag>
92{
93/* JJV changed expmax (log(1.0E38)==87.49823), and added minlog */
94
95 // TV: The original code had infnty = 1.0E38, minlog = 1.0E-37.
96 // I have changed these to 0.3 * huge and 8.5 * tiny. For IEEE
97 // float this works out to 1.020847E38 and 0.9991702E-37.
98 // I changed expmax from 87.49823 to log(infnty), which works out
99 // to 87.518866 for float.
100
101 static T olda = minlog;
102 static T oldb = minlog;
103 static T genbet,a,alpha,b,beta,delta,gamma,k1,k2,r,s,t,u1,u2,v,w,y,z;
104 static long qsame;
105
106 // This code determines if the aa and bb parameters are unchanged.
107 // If so, some code can be skipped.
108
109 qsame = (olda == aa) && (oldb == bb);
110
111 if (!qsame)
112 {
113 BZPRECHECK((aa > minlog) && (bb > minlog),
114 "Beta RNG: parameters must be >= " << minlog << endl
115 << "Parameters are aa=" << aa << " and bb=" << bb << endl);
116 olda = aa;
117 oldb = bb;
118 }
119
120 if (!(min(aa,bb) > 1.0))
121 goto S100;
122/*
123 Alborithm BB
124 Initialize
125*/
126 if(qsame) goto S30;
127 a = min(aa,bb);
128 b = max(aa,bb);
129 alpha = a+b;
130 beta = sqrt((alpha-2.0)/(2.0*a*b-alpha));
131 gamma = a+1.0/beta;
132S30:
133S40:
134 u1 = ranf();
135/*
136 Step 1
137*/
138 u2 = ranf();
139 v = beta*log(u1/(1.0-u1));
140/* JJV altered this */
141 if(v > expmax) goto S55;
142/*
143 * JJV added checker to see if a*exp(v) will overflow
144 * JJV S50 _was_ w = a*exp(v); also note here a > 1.0
145 */
146 w = exp(v);
147 if(w > infnty/a) goto S55;
148 w *= a;
149 goto S60;
150S55:
151 w = infnty;
152S60:
153 z = pow(u1,2.0)*u2;
154 r = gamma*v-1.3862944;
155 s = a+r-w;
156/*
157 Step 2
158*/
159 if(s+2.609438 >= 5.0*z) goto S70;
160/*
161 Step 3
162*/
163 t = log(z);
164 if(s > t) goto S70;
165/*
166 * Step 4
167 *
168 * JJV added checker to see if log(alpha/(b+w)) will
169 * JJV overflow. If so, we count the log as -INF, and
170 * JJV consequently evaluate conditional as true, i.e.
171 * JJV the algorithm rejects the trial and starts over
172 * JJV May not need this here since alpha > 2.0
173 */
174 if(alpha/(b+w) < minlog) goto S40;
175 if(r+alpha*log(alpha/(b+w)) < t) goto S40;
176S70:
177/*
178 Step 5
179*/
180 if(!(aa == a)) goto S80;
181 genbet = w/(b+w);
182 goto S90;
183S80:
184 genbet = b/(b+w);
185S90:
186 goto S230;
187S100:
188/*
189 Algorithm BC
190 Initialize
191*/
192 if(qsame) goto S110;
193 a = max(aa,bb);
194 b = min(aa,bb);
195 alpha = a+b;
196 beta = 1.0/b;
197 delta = 1.0+a-b;
198 k1 = delta*(1.38889E-2+4.16667E-2*b)/(a*beta-0.777778);
199 k2 = 0.25+(0.5+0.25/delta)*b;
200S110:
201S120:
202 u1 = ranf();
203/*
204 Step 1
205*/
206 u2 = ranf();
207 if(u1 >= 0.5) goto S130;
208/*
209 Step 2
210*/
211 y = u1*u2;
212 z = u1*y;
213 if(0.25*u2+z-y >= k1) goto S120;
214 goto S170;
215S130:
216/*
217 Step 3
218*/
219 z = pow(u1,2.0)*u2;
220 if(!(z <= 0.25)) goto S160;
221 v = beta*log(u1/(1.0-u1));
222/*
223 * JJV instead of checking v > expmax at top, I will check
224 * JJV if a < 1, then check the appropriate values
225 */
226 if(a > 1.0) goto S135;
227/* JJV a < 1 so it can help out if exp(v) would overflow */
228 if(v > expmax) goto S132;
229 w = a*exp(v);
230 goto S200;
231S132:
232 w = v + log(a);
233 if(w > expmax) goto S140;
234 w = exp(w);
235 goto S200;
236S135:
237/* JJV in this case a > 1 */
238 if(v > expmax) goto S140;
239 w = exp(v);
240 if(w > infnty/a) goto S140;
241 w *= a;
242 goto S200;
243S140:
244 w = infnty;
245 goto S200;
246/*
247 * JJV old code
248 * if(!(v > expmax)) goto S140;
249 * w = infnty;
250 * goto S150;
251 *S140:
252 * w = a*exp(v);
253 *S150:
254 * goto S200;
255 */
256S160:
257 if(z >= k2) goto S120;
258S170:
259/*
260 Step 4
261 Step 5
262*/
263 v = beta*log(u1/(1.0-u1));
264/* JJV same kind of checking as above */
265 if(a > 1.0) goto S175;
266/* JJV a < 1 so it can help out if exp(v) would overflow */
267 if(v > expmax) goto S172;
268 w = a*exp(v);
269 goto S190;
270S172:
271 w = v + log(a);
272 if(w > expmax) goto S180;
273 w = exp(w);
274 goto S190;
275S175:
276/* JJV in this case a > 1.0 */
277 if(v > expmax) goto S180;
278 w = exp(v);
279 if(w > infnty/a) goto S180;
280 w *= a;
281 goto S190;
282S180:
283 w = infnty;
284/*
285 * JJV old code
286 * if(!(v > expmax)) goto S180;
287 * w = infnty;
288 * goto S190;
289 *S180:
290 * w = a*exp(v);
291 */
292S190:
293/*
294 * JJV here we also check to see if log overlows; if so, we treat it
295 * JJV as -INF, which means condition is true, i.e. restart
296 */
297 if(alpha/(b+w) < minlog) goto S120;
298 if(alpha*(log(alpha/(b+w))+v)-1.3862944 < log(z)) goto S120;
299S200:
300/*
301 Step 6
302*/
303 if(!(a == aa)) goto S210;
304 genbet = w/(b+w);
305 goto S220;
306S210:
307 genbet = b/(b+w);
308S230:
309S220:
310 return genbet;
311}
312
313}
314
315#endif // BZ_RANDOM_BETA
Definition memblock.h:307
Definition beta.h:55
T random()
Definition beta.h:91
T minlog
Definition beta.h:87
void setParameters(T a, T b)
Definition beta.h:71
Beta(T a, T b)
Definition beta.h:59
T aa
Definition beta.h:86
Beta(T a, T b, unsigned int i)
Definition beta.h:64
T T_numtype
Definition beta.h:57
T infnty
Definition beta.h:87
T ranf()
Definition beta.h:81
T expmax
Definition beta.h:87
T bb
Definition beta.h:86
Definition uniform.h:377
T random()
Definition uniform.h:385
Definition beta.h:50
MersenneTwister defaultIRNG
Definition default.h:120
Definition default.h:53