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