blitz  Version 1.0.2
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 
50 namespace ranlib {
51 
52 template<typename T = double, typename IRNG = defaultIRNG,
53  typename stateTag = defaultState>
54 class Beta : public UniformOpen<T,IRNG,stateTag>
55 {
56 public:
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 
80 protected:
81  T ranf()
82  {
84  }
85 
86  T aa, bb;
88 };
89 
90 template<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;
132 S30:
133 S40:
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;
150 S55:
151  w = infnty;
152 S60:
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;
176 S70:
177 /*
178  Step 5
179 */
180  if(!(aa == a)) goto S80;
181  genbet = w/(b+w);
182  goto S90;
183 S80:
184  genbet = b/(b+w);
185 S90:
186  goto S230;
187 S100:
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;
200 S110:
201 S120:
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;
215 S130:
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;
231 S132:
232  w = v + log(a);
233  if(w > expmax) goto S140;
234  w = exp(w);
235  goto S200;
236 S135:
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;
243 S140:
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  */
256 S160:
257  if(z >= k2) goto S120;
258 S170:
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;
270 S172:
271  w = v + log(a);
272  if(w > expmax) goto S180;
273  w = exp(w);
274  goto S190;
275 S175:
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;
282 S180:
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  */
292 S190:
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;
299 S200:
300 /*
301  Step 6
302 */
303  if(!(a == aa)) goto S210;
304  genbet = w/(b+w);
305  goto S220;
306 S210:
307  genbet = b/(b+w);
308 S230:
309 S220:
310  return genbet;
311 }
312 
313 }
314 
315 #endif // BZ_RANDOM_BETA
T expmax
Definition: beta.h:87
_bz_global blitz::IndexPlaceholder< 9 > r
Definition: indexexpr.h:265
Definition: beta.h:50
Definition: beta.h:54
T random()
Definition: uniform.h:385
_bz_global blitz::IndexPlaceholder< 11 > t
Definition: indexexpr.h:267
_bz_global blitz::IndexPlaceholder< 0 > i
Definition: indexexpr.h:256
T bb
Definition: beta.h:86
T infnty
Definition: beta.h:87
Definition: default.h:53
MersenneTwister defaultIRNG
Definition: default.h:120
T aa
Definition: beta.h:86
Definition: uniform.h:377
N_length & a
Definition: tvecglobs.h:47
T T_numtype
Definition: beta.h:57
Beta(T a, T b, unsigned int i)
Definition: beta.h:64
const T2 & b
Definition: minmax.h:48
Beta(T a, T b)
Definition: beta.h:59
T random()
Definition: beta.h:91
T ranf()
Definition: beta.h:81
_bz_global blitz::IndexPlaceholder< 10 > s
Definition: indexexpr.h:266
T minlog
Definition: beta.h:87
void setParameters(T a, T b)
Definition: beta.h:71