random/gamma.h

Go to the documentation of this file.
00001 /*
00002  * Gamma distribution
00003  *
00004  * Source: Ahrens, J.H. and Dieter, U., Generating Gamma variates
00005  * by a modified rejection technique.  Comm. ACM, 25,1 (Jan. 1982) 
00006  * pp. 47-54.
00007  *
00008  * This code has been adapted from RANDLIB.C 1.3, by
00009  * Barry W. Brown, James Lovato, Kathy Russell, and John Venier.
00010  * Code was originally by Ahrens and Dieter (see above).
00011  *
00012  * Adapter's notes:
00013  * NEEDS_WORK: more precision for literals.
00014  * NEEDS_WORK: ideally the normal_ member should be driven from
00015  * the same IRNG as the Gamma object, in the event that independentState
00016  * is used.  Not clear how this could be accomplished.
00017  */
00018 
00019 #ifndef BZ_RANDOM_GAMMA
00020 #define BZ_RANDOM_GAMMA
00021 
00022 #ifndef BZ_RANDOM_UNIFORM
00023  #include <random/uniform.h>
00024 #endif
00025 
00026 #ifndef BZ_RANDOM_NORMAL
00027  #include <random/normal.h>
00028 #endif
00029 
00030 #ifndef BZ_RANDOM_EXPONENTIAL
00031  #include <random/exponential.h>
00032 #endif
00033 
00034 #ifndef BZ_NUMINQUIRE_H
00035  #include <blitz/numinquire.h>
00036 #endif
00037 
00038 BZ_NAMESPACE(ranlib)
00039 
00040 template<typename T = double, typename IRNG = defaultIRNG, 
00041     typename stateTag = defaultState>
00042 class Gamma : public UniformOpen<T,IRNG,stateTag>
00043 {
00044 public:
00045     typedef T T_numtype;
00046 
00047     Gamma()
00048     {
00049         setMean(1.0);
00050     }
00051 
00052     Gamma(T mean)
00053     {
00054         setMean(mean);
00055     }
00056 
00057     T random();
00058 
00059     void setMean(T mean)
00060     {
00061         BZPRECONDITION(mean >= 1.0);
00062         a = mean;
00063     }
00064 
00065 protected:
00066     T ranf() 
00067     { 
00068         return UniformOpen<T,IRNG,stateTag>::random(); 
00069     }
00070 
00071     T snorm()
00072     {
00073         return normal_.random();
00074     }
00075 
00076     T sexpo()
00077     {
00078         return exponential_.random();
00079     }
00080 
00081     T fsign(T num, T sign)
00082     {
00083         /* Transfers sign of argument sign to argument num */
00084 
00085         if ((sign>0.0L && num<0.0L)||(sign<0.0L && num>0.0L))
00086             return -num;
00087         else 
00088             return num;
00089     }
00090 
00091     NormalUnit<T,IRNG,sharedState> normal_;
00092     ExponentialUnit<T,IRNG,sharedState> exponential_;
00093 
00094     T a;
00095 };
00096 
00097 template<typename T, typename IRNG, typename stateTag>
00098 T Gamma<T,IRNG,stateTag>::random()
00099 {
00100     /*
00101      INPUT: A =PARAMETER (MEAN) OF THE STANDARD GAMMA DISTRIBUTION
00102      OUTPUT: SGAMMA = SAMPLE FROM THE GAMMA-(A)-DISTRIBUTION
00103      COEFFICIENTS Q(K) - FOR Q0 = SUM(Q(K)*A**(-K))
00104      COEFFICIENTS A(K) - FOR Q = Q0+(T*T/2)*SUM(A(K)*V**K)
00105      COEFFICIENTS E(K) - FOR EXP(Q)-1 = SUM(E(K)*Q**K)
00106      PREVIOUS A PRE-SET TO ZERO - AA IS A', AAA IS A"
00107      SQRT32 IS THE SQUAREROOT OF 32 = 5.656854249492380
00108      */
00109 
00110 static T q1 = 4.166669E-2;
00111 static T q2 = 2.083148E-2;
00112 static T q3 = 8.01191E-3;
00113 static T q4 = 1.44121E-3;
00114 static T q5 = -7.388E-5;
00115 static T q6 = 2.4511E-4;
00116 static T q7 = 2.424E-4;
00117 static T a1 = 0.3333333;
00118 static T a2 = -0.250003;
00119 static T a3 = 0.2000062;
00120 static T a4 = -0.1662921;
00121 static T a5 = 0.1423657;
00122 static T a6 = -0.1367177;
00123 static T a7 = 0.1233795;
00124 static T e1 = 1.0;
00125 static T e2 = 0.4999897;
00126 static T e3 = 0.166829;
00127 static T e4 = 4.07753E-2;
00128 static T e5 = 1.0293E-2;
00129 static T aa = 0.0;
00130 static T aaa = 0.0;
00131 static T sqrt32 = 5.656854249492380195206754896838792314280;
00132 
00133 /* JJV added b0 to fix rare and subtle bug */
00134 static T sgamma,s2,s,d,t,x,u,r,q0,b,b0,si,c,v,q,e,w,p;
00135 
00136     if(a == aa) goto S10;
00137     if(a < 1.0) goto S120;
00138 /*
00139      STEP  1:  RECALCULATIONS OF S2,S,D IF A HAS CHANGED
00140 */
00141     aa = a;
00142     s2 = a-0.5;
00143     s = sqrt(s2);
00144     d = sqrt32-12.0*s;
00145 S10:
00146 /*
00147      STEP  2:  T=STANDARD NORMAL DEVIATE,
00148                X=(S,1/2)-NORMAL DEVIATE.
00149                IMMEDIATE ACCEPTANCE (I)
00150 */
00151     t = snorm();
00152     x = s+0.5*t;
00153     sgamma = x*x;
00154     if(t >= 0.0) return sgamma;
00155 /*
00156      STEP  3:  U= 0,1 -UNIFORM SAMPLE. SQUEEZE ACCEPTANCE (S)
00157 */
00158     u = ranf();
00159     if(d*u <= t*t*t) return sgamma;
00160 /*
00161      STEP  4:  RECALCULATIONS OF Q0,B,SI,C IF NECESSARY
00162 */
00163     if(a == aaa) goto S40;
00164     aaa = a;
00165     r = 1.0/ a;
00166     q0 = ((((((q7*r+q6)*r+q5)*r+q4)*r+q3)*r+q2)*r+q1)*r;
00167 /*
00168                APPROXIMATION DEPENDING ON SIZE OF PARAMETER A
00169                THE CONSTANTS IN THE EXPRESSIONS FOR B, SI AND
00170                C WERE ESTABLISHED BY NUMERICAL EXPERIMENTS
00171 */
00172     if(a <= 3.686) goto S30;
00173     if(a <= 13.022) goto S20;
00174 /*
00175                CASE 3:  A .GT. 13.022
00176 */
00177     b = 1.77;
00178     si = 0.75;
00179     c = 0.1515/s;
00180     goto S40;
00181 S20:
00182 /*
00183                CASE 2:  3.686 .LT. A .LE. 13.022
00184 */
00185     b = 1.654+7.6E-3*s2;
00186     si = 1.68/s+0.275;
00187     c = 6.2E-2/s+2.4E-2;
00188     goto S40;
00189 S30:
00190 /*
00191                CASE 1:  A .LE. 3.686
00192 */
00193     b = 0.463+s+0.178*s2;
00194     si = 1.235;
00195     c = 0.195/s-7.9E-2+1.6E-1*s;
00196 S40:
00197 /*
00198      STEP  5:  NO QUOTIENT TEST IF X NOT POSITIVE
00199 */
00200     if(x <= 0.0) goto S70;
00201 /*
00202      STEP  6:  CALCULATION OF V AND QUOTIENT Q
00203 */
00204     v = t/(s+s);
00205     if(fabs(v) <= 0.25) goto S50;
00206     q = q0-s*t+0.25*t*t+(s2+s2)*log(1.0+v);
00207     goto S60;
00208 S50:
00209     q = q0+0.5*t*t*((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v;
00210 S60:
00211 /*
00212      STEP  7:  QUOTIENT ACCEPTANCE (Q)
00213 */
00214     if(log(1.0-u) <= q) return sgamma;
00215 S70:
00216 /*
00217      STEP  8:  E=STANDARD EXPONENTIAL DEVIATE
00218                U= 0,1 -UNIFORM DEVIATE
00219                T=(B,SI)-DOUBLE EXPONENTIAL (LAPLACE) SAMPLE
00220 */
00221     e = sexpo();
00222     u = ranf();
00223     u += (u-1.0);
00224     t = b+fsign(si*e,u);
00225 /*
00226      STEP  9:  REJECTION IF T .LT. TAU(1) = -.71874483771719
00227 */
00228     if(t < -0.7187449) goto S70;
00229 /*
00230      STEP 10:  CALCULATION OF V AND QUOTIENT Q
00231 */
00232     v = t/(s+s);
00233     if(fabs(v) <= 0.25) goto S80;
00234     q = q0-s*t+0.25*t*t+(s2+s2)*log(1.0+v);
00235     goto S90;
00236 S80:
00237     q = q0+0.5*t*t*((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v;
00238 S90:
00239 /*
00240      STEP 11:  HAT ACCEPTANCE (H) (IF Q NOT POSITIVE GO TO STEP 8)
00241 */
00242     if(q <= 0.0) goto S70;
00243     if(q <= 0.5) goto S100;
00244 /*
00245  * JJV modified the code through line 115 to handle large Q case
00246  */
00247     if(q < 15.0) goto S95;
00248 /*
00249  * JJV Here Q is large enough that Q = log(exp(Q) - 1.0) (for real Q)
00250  * JJV so reformulate test at 110 in terms of one EXP, if not too big
00251  * JJV 87.49823 is close to the largest real which can be
00252  * JJV exponentiated (87.49823 = log(1.0E38))
00253  */
00254     if((q+e-0.5*t*t) > 87.49823) goto S115;
00255     if(c*fabs(u) > exp(q+e-0.5*t*t)) goto S70;
00256     goto S115;
00257 S95:
00258     w = exp(q)-1.0;
00259     goto S110;
00260 S100:
00261     w = ((((e5*q+e4)*q+e3)*q+e2)*q+e1)*q;
00262 S110:
00263 /*
00264                IF T IS REJECTED, SAMPLE AGAIN AT STEP 8
00265 */
00266     if(c*fabs(u) > w*exp(e-0.5*t*t)) goto S70;
00267 S115:
00268     x = s+0.5*t;
00269     sgamma = x*x;
00270     return sgamma;
00271 S120:
00272 /*
00273      ALTERNATE METHOD FOR PARAMETERS A BELOW 1  (.3678794=EXP(-1.))
00274 
00275      JJV changed B to B0 (which was added to declarations for this)
00276      JJV in 120 to END to fix rare and subtle bug.
00277      JJV Line: 'aa = 0.0' was removed (unnecessary, wasteful).
00278      JJV Reasons: the state of AA only serves to tell the A >= 1.0
00279      JJV case if certain A-dependent constants need to be recalculated.
00280      JJV The A < 1.0 case (here) no longer changes any of these, and
00281      JJV the recalculation of B (which used to change with an
00282      JJV A < 1.0 call) is governed by the state of AAA anyway.
00283     aa = 0.0;
00284 */
00285     b0 = 1.0+0.3678794*a;
00286 S130:
00287     p = b0*ranf();
00288     if(p >= 1.0) goto S140;
00289     sgamma = exp(log(p)/ a);
00290     if(sexpo() < sgamma) goto S130;
00291     return sgamma;
00292 S140:
00293     sgamma = -log((b0-p)/ a);
00294     if(sexpo() < (1.0-a)*log(sgamma)) goto S130;
00295     return sgamma;
00296 
00297 }
00298 
00299 BZ_NAMESPACE_END
00300 
00301 #endif // BZ_RANDOM_GAMMA

Generated on Wed Oct 17 15:37:28 2007 for blitz by  doxygen 1.5.2