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