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