blitz Version 0.10
|
00001 // -*- C++ -*- 00002 // $Id: beta.h,v 1.4 2011/03/25 22:41:17 julianc Exp $ 00003 00004 /* 00005 * Generate Beta random deviate 00006 * 00007 * Returns a single random deviate from the beta distribution with 00008 * parameters A and B. The density of the beta is 00009 * x^(a-1) * (1-x)^(b-1) / B(a,b) for 0 < x < 1 00010 * 00011 * The mean is a/(a+b). 00012 * The variance is ab/((a+b)^2(a+b+1)) 00013 * The rth moment is (a+r-1)^(r)/(a+b+r-1)^(r) 00014 * where a^(r) means a * (a-1) * (a-2) * ... * (a-r+1) 00015 * 00016 * Method 00017 * R. C. H. Cheng 00018 * Generating Beta Variates with Nonintegral Shape Parameters 00019 * Communications of the ACM, 21:317-322 (1978) 00020 * (Algorithms BB and BC) 00021 * http://www.acm.org/pubs/citations/journals/toms/1991-17-1/p98-l_ecuyer/ 00022 * 00023 * This class has been adapted from RANDLIB.C 1.3, by 00024 * Barry W. Brown, James Lovato, Kathy Russell, and John Venier. 00025 * 00026 * Adapter's note (TV): This code has gone from Pascal to Fortran to C. 00027 * As a result it is a bit of a mess. Note also that the algorithms were 00028 * originally designed for 32-bit float, and so some of the constants 00029 * below have inadequate precision. This will not cause problems for 00030 * casual use, but if you are generating millions of beta variates and 00031 * rely on some convergence property, you may have want to worry 00032 * about this. 00033 * 00034 * NEEDS_WORK: dig out the original paper and determine these constants 00035 * to precision adequate for 128-bit float. 00036 * NEEDS_WORK: turn this into structured code. 00037 */ 00038 00039 #ifndef BZ_RANDOM_BETA 00040 #define BZ_RANDOM_BETA 00041 00042 #ifndef BZ_RANDOM_UNIFORM 00043 #include <random/uniform.h> 00044 #endif 00045 00046 #ifndef BZ_NUMINQUIRE_H 00047 #include <blitz/numinquire.h> 00048 #endif 00049 00050 BZ_NAMESPACE(ranlib) 00051 00052 template<typename T = double, typename IRNG = defaultIRNG, 00053 typename stateTag = defaultState> 00054 class Beta : public UniformOpen<T,IRNG,stateTag> 00055 { 00056 public: 00057 typedef T T_numtype; 00058 00059 Beta(T a, T b) 00060 { 00061 setParameters(a, b); 00062 } 00063 00064 Beta(T a, T b, unsigned int i) : UniformOpen<T, IRNG, stateTag>(i) 00065 { 00066 setParameters(a, b); 00067 } 00068 00069 T random(); 00070 00071 void setParameters(T a, T b) 00072 { 00073 aa = a; 00074 bb = b; 00075 infnty = 0.3 * blitz::huge(T()); 00076 minlog = 0.085 * blitz::tiny(T()); 00077 expmax = log(infnty); 00078 } 00079 00080 protected: 00081 T ranf() 00082 { 00083 return UniformOpen<T,IRNG,stateTag>::random(); 00084 } 00085 00086 T aa, bb; 00087 T infnty, minlog, expmax; 00088 }; 00089 00090 template<typename T, typename IRNG, typename stateTag> 00091 T Beta<T,IRNG,stateTag>::random() 00092 { 00093 /* JJV changed expmax (log(1.0E38)==87.49823), and added minlog */ 00094 00095 // TV: The original code had infnty = 1.0E38, minlog = 1.0E-37. 00096 // I have changed these to 0.3 * huge and 8.5 * tiny. For IEEE 00097 // float this works out to 1.020847E38 and 0.9991702E-37. 00098 // I changed expmax from 87.49823 to log(infnty), which works out 00099 // to 87.518866 for float. 00100 00101 static T olda = minlog; 00102 static T oldb = minlog; 00103 static T genbet,a,alpha,b,beta,delta,gamma,k1,k2,r,s,t,u1,u2,v,w,y,z; 00104 static long qsame; 00105 00106 // This code determines if the aa and bb parameters are unchanged. 00107 // If so, some code can be skipped. 00108 00109 qsame = (olda == aa) && (oldb == bb); 00110 00111 if (!qsame) 00112 { 00113 BZPRECHECK((aa > minlog) && (bb > minlog), 00114 "Beta RNG: parameters must be >= " << minlog << endl 00115 << "Parameters are aa=" << aa << " and bb=" << bb << endl); 00116 olda = aa; 00117 oldb = bb; 00118 } 00119 00120 if (!(min(aa,bb) > 1.0)) 00121 goto S100; 00122 /* 00123 Alborithm BB 00124 Initialize 00125 */ 00126 if(qsame) goto S30; 00127 a = min(aa,bb); 00128 b = max(aa,bb); 00129 alpha = a+b; 00130 beta = sqrt((alpha-2.0)/(2.0*a*b-alpha)); 00131 gamma = a+1.0/beta; 00132 S30: 00133 S40: 00134 u1 = ranf(); 00135 /* 00136 Step 1 00137 */ 00138 u2 = ranf(); 00139 v = beta*log(u1/(1.0-u1)); 00140 /* JJV altered this */ 00141 if(v > expmax) goto S55; 00142 /* 00143 * JJV added checker to see if a*exp(v) will overflow 00144 * JJV S50 _was_ w = a*exp(v); also note here a > 1.0 00145 */ 00146 w = exp(v); 00147 if(w > infnty/a) goto S55; 00148 w *= a; 00149 goto S60; 00150 S55: 00151 w = infnty; 00152 S60: 00153 z = pow(u1,2.0)*u2; 00154 r = gamma*v-1.3862944; 00155 s = a+r-w; 00156 /* 00157 Step 2 00158 */ 00159 if(s+2.609438 >= 5.0*z) goto S70; 00160 /* 00161 Step 3 00162 */ 00163 t = log(z); 00164 if(s > t) goto S70; 00165 /* 00166 * Step 4 00167 * 00168 * JJV added checker to see if log(alpha/(b+w)) will 00169 * JJV overflow. If so, we count the log as -INF, and 00170 * JJV consequently evaluate conditional as true, i.e. 00171 * JJV the algorithm rejects the trial and starts over 00172 * JJV May not need this here since alpha > 2.0 00173 */ 00174 if(alpha/(b+w) < minlog) goto S40; 00175 if(r+alpha*log(alpha/(b+w)) < t) goto S40; 00176 S70: 00177 /* 00178 Step 5 00179 */ 00180 if(!(aa == a)) goto S80; 00181 genbet = w/(b+w); 00182 goto S90; 00183 S80: 00184 genbet = b/(b+w); 00185 S90: 00186 goto S230; 00187 S100: 00188 /* 00189 Algorithm BC 00190 Initialize 00191 */ 00192 if(qsame) goto S110; 00193 a = max(aa,bb); 00194 b = min(aa,bb); 00195 alpha = a+b; 00196 beta = 1.0/b; 00197 delta = 1.0+a-b; 00198 k1 = delta*(1.38889E-2+4.16667E-2*b)/(a*beta-0.777778); 00199 k2 = 0.25+(0.5+0.25/delta)*b; 00200 S110: 00201 S120: 00202 u1 = ranf(); 00203 /* 00204 Step 1 00205 */ 00206 u2 = ranf(); 00207 if(u1 >= 0.5) goto S130; 00208 /* 00209 Step 2 00210 */ 00211 y = u1*u2; 00212 z = u1*y; 00213 if(0.25*u2+z-y >= k1) goto S120; 00214 goto S170; 00215 S130: 00216 /* 00217 Step 3 00218 */ 00219 z = pow(u1,2.0)*u2; 00220 if(!(z <= 0.25)) goto S160; 00221 v = beta*log(u1/(1.0-u1)); 00222 /* 00223 * JJV instead of checking v > expmax at top, I will check 00224 * JJV if a < 1, then check the appropriate values 00225 */ 00226 if(a > 1.0) goto S135; 00227 /* JJV a < 1 so it can help out if exp(v) would overflow */ 00228 if(v > expmax) goto S132; 00229 w = a*exp(v); 00230 goto S200; 00231 S132: 00232 w = v + log(a); 00233 if(w > expmax) goto S140; 00234 w = exp(w); 00235 goto S200; 00236 S135: 00237 /* JJV in this case a > 1 */ 00238 if(v > expmax) goto S140; 00239 w = exp(v); 00240 if(w > infnty/a) goto S140; 00241 w *= a; 00242 goto S200; 00243 S140: 00244 w = infnty; 00245 goto S200; 00246 /* 00247 * JJV old code 00248 * if(!(v > expmax)) goto S140; 00249 * w = infnty; 00250 * goto S150; 00251 *S140: 00252 * w = a*exp(v); 00253 *S150: 00254 * goto S200; 00255 */ 00256 S160: 00257 if(z >= k2) goto S120; 00258 S170: 00259 /* 00260 Step 4 00261 Step 5 00262 */ 00263 v = beta*log(u1/(1.0-u1)); 00264 /* JJV same kind of checking as above */ 00265 if(a > 1.0) goto S175; 00266 /* JJV a < 1 so it can help out if exp(v) would overflow */ 00267 if(v > expmax) goto S172; 00268 w = a*exp(v); 00269 goto S190; 00270 S172: 00271 w = v + log(a); 00272 if(w > expmax) goto S180; 00273 w = exp(w); 00274 goto S190; 00275 S175: 00276 /* JJV in this case a > 1.0 */ 00277 if(v > expmax) goto S180; 00278 w = exp(v); 00279 if(w > infnty/a) goto S180; 00280 w *= a; 00281 goto S190; 00282 S180: 00283 w = infnty; 00284 /* 00285 * JJV old code 00286 * if(!(v > expmax)) goto S180; 00287 * w = infnty; 00288 * goto S190; 00289 *S180: 00290 * w = a*exp(v); 00291 */ 00292 S190: 00293 /* 00294 * JJV here we also check to see if log overlows; if so, we treat it 00295 * JJV as -INF, which means condition is true, i.e. restart 00296 */ 00297 if(alpha/(b+w) < minlog) goto S120; 00298 if(alpha*(log(alpha/(b+w))+v)-1.3862944 < log(z)) goto S120; 00299 S200: 00300 /* 00301 Step 6 00302 */ 00303 if(!(a == aa)) goto S210; 00304 genbet = w/(b+w); 00305 goto S220; 00306 S210: 00307 genbet = b/(b+w); 00308 S230: 00309 S220: 00310 return genbet; 00311 } 00312 00313 BZ_NAMESPACE_END 00314 00315 #endif // BZ_RANDOM_BETA