blitz Version 0.10
random/beta.h
Go to the documentation of this file.
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
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Defines