beta.h

Go to the documentation of this file.
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

Generated on Mon Oct 22 10:04:52 2007 for blitz by  doxygen 1.3.9.1