blitz Version 0.10
|
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