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