blitz Version 1.0.2
Loading...
Searching...
No Matches
gamma.h
Go to the documentation of this file.
1// -*- C++ -*-
2// $Id$
3
4/*
5 * Gamma distribution
6 *
7 * Source: Ahrens, J.H. and Dieter, U., Generating Gamma variates
8 * by a modified rejection technique. Comm. ACM, 25,1 (Jan. 1982)
9 * pp. 47-54.
10 *
11 * This code has been adapted from RANDLIB.C 1.3, by
12 * Barry W. Brown, James Lovato, Kathy Russell, and John Venier.
13 * Code was originally by Ahrens and Dieter (see above).
14 *
15 * Adapter's notes:
16 * NEEDS_WORK: more precision for literals.
17 * NEEDS_WORK: ideally the normal_ member should be driven from
18 * the same IRNG as the Gamma object, in the event that independentState
19 * is used. Not clear how this could be accomplished.
20 */
21
22#ifndef BZ_RANDOM_GAMMA
23#define BZ_RANDOM_GAMMA
24
25#ifndef BZ_RANDOM_UNIFORM
26 #include <random/uniform.h>
27#endif
28
29#ifndef BZ_RANDOM_NORMAL
30 #include <random/normal.h>
31#endif
32
33#ifndef BZ_RANDOM_EXPONENTIAL
34 #include <random/exponential.h>
35#endif
36
37#ifndef BZ_NUMINQUIRE_H
38 #include <blitz/numinquire.h>
39#endif
40
41namespace ranlib {
42
43template<typename T = double, typename IRNG = defaultIRNG,
44 typename stateTag = defaultState>
45class Gamma : public UniformOpen<T,IRNG,stateTag>
46{
47public:
48 typedef T T_numtype;
49
51 {
52 setMean(1.0);
53 }
54
55 explicit Gamma(unsigned int i) :
56 UniformOpen<T,IRNG,stateTag>(i)
57 {
58 setMean(1.0);
59 };
60
61 Gamma(T mean)
62 {
63 setMean(mean);
64 }
65
66 Gamma(T mean, unsigned int i) :
67 UniformOpen<T,IRNG,stateTag>(i)
68 {
69 setMean(mean);
70 };
71
72 T random();
73
74 void setMean(T mean)
75 {
76 BZPRECONDITION(mean >= 1.0);
77 a = mean;
78 }
79
80protected:
81 T ranf()
82 {
84 }
85
87 {
88 return normal_.random();
89 }
90
92 {
93 return exponential_.random();
94 }
95
96 T fsign(T num, T sign)
97 {
98 /* Transfers sign of argument sign to argument num */
99
100 if ((sign>0.0L && num<0.0L)||(sign<0.0L && num>0.0L))
101 return -num;
102 else
103 return num;
104 }
105
108
109 T a;
110};
111
112template<typename T, typename IRNG, typename stateTag>
114{
115 /*
116 INPUT: A =PARAMETER (MEAN) OF THE STANDARD GAMMA DISTRIBUTION
117 OUTPUT: SGAMMA = SAMPLE FROM THE GAMMA-(A)-DISTRIBUTION
118 COEFFICIENTS Q(K) - FOR Q0 = SUM(Q(K)*A**(-K))
119 COEFFICIENTS A(K) - FOR Q = Q0+(T*T/2)*SUM(A(K)*V**K)
120 COEFFICIENTS E(K) - FOR EXP(Q)-1 = SUM(E(K)*Q**K)
121 PREVIOUS A PRE-SET TO ZERO - AA IS A', AAA IS A"
122 SQRT32 IS THE SQUAREROOT OF 32 = 5.656854249492380
123 */
124
125static T q1 = 4.166669E-2;
126static T q2 = 2.083148E-2;
127static T q3 = 8.01191E-3;
128static T q4 = 1.44121E-3;
129static T q5 = -7.388E-5;
130static T q6 = 2.4511E-4;
131static T q7 = 2.424E-4;
132static T a1 = 0.3333333;
133static T a2 = -0.250003;
134static T a3 = 0.2000062;
135static T a4 = -0.1662921;
136static T a5 = 0.1423657;
137static T a6 = -0.1367177;
138static T a7 = 0.1233795;
139static T e1 = 1.0;
140static T e2 = 0.4999897;
141static T e3 = 0.166829;
142static T e4 = 4.07753E-2;
143static T e5 = 1.0293E-2;
144static T aa = 0.0;
145static T aaa = 0.0;
146static T sqrt32 = 5.656854249492380195206754896838792314280;
147
148/* JJV added b0 to fix rare and subtle bug */
149static T sgamma,s2,s,d,t,x,u,r,q0,b,b0,si,c,v,q,e,w,p;
150
151 if(a == aa) goto S10;
152 if(a < 1.0) goto S120;
153/*
154 STEP 1: RECALCULATIONS OF S2,S,D IF A HAS CHANGED
155*/
156 aa = a;
157 s2 = a-0.5;
158 s = sqrt(s2);
159 d = sqrt32-12.0*s;
160S10:
161/*
162 STEP 2: T=STANDARD NORMAL DEVIATE,
163 X=(S,1/2)-NORMAL DEVIATE.
164 IMMEDIATE ACCEPTANCE (I)
165*/
166 t = snorm();
167 x = s+0.5*t;
168 sgamma = x*x;
169 if(t >= 0.0) return sgamma;
170/*
171 STEP 3: U= 0,1 -UNIFORM SAMPLE. SQUEEZE ACCEPTANCE (S)
172*/
173 u = ranf();
174 if(d*u <= t*t*t) return sgamma;
175/*
176 STEP 4: RECALCULATIONS OF Q0,B,SI,C IF NECESSARY
177*/
178 if(a == aaa) goto S40;
179 aaa = a;
180 r = 1.0/ a;
181 q0 = ((((((q7*r+q6)*r+q5)*r+q4)*r+q3)*r+q2)*r+q1)*r;
182/*
183 APPROXIMATION DEPENDING ON SIZE OF PARAMETER A
184 THE CONSTANTS IN THE EXPRESSIONS FOR B, SI AND
185 C WERE ESTABLISHED BY NUMERICAL EXPERIMENTS
186*/
187 if(a <= 3.686) goto S30;
188 if(a <= 13.022) goto S20;
189/*
190 CASE 3: A .GT. 13.022
191*/
192 b = 1.77;
193 si = 0.75;
194 c = 0.1515/s;
195 goto S40;
196S20:
197/*
198 CASE 2: 3.686 .LT. A .LE. 13.022
199*/
200 b = 1.654+7.6E-3*s2;
201 si = 1.68/s+0.275;
202 c = 6.2E-2/s+2.4E-2;
203 goto S40;
204S30:
205/*
206 CASE 1: A .LE. 3.686
207*/
208 b = 0.463+s+0.178*s2;
209 si = 1.235;
210 c = 0.195/s-7.9E-2+1.6E-1*s;
211S40:
212/*
213 STEP 5: NO QUOTIENT TEST IF X NOT POSITIVE
214*/
215 if(x <= 0.0) goto S70;
216/*
217 STEP 6: CALCULATION OF V AND QUOTIENT Q
218*/
219 v = t/(s+s);
220 if(fabs(v) <= 0.25) goto S50;
221 q = q0-s*t+0.25*t*t+(s2+s2)*log(1.0+v);
222 goto S60;
223S50:
224 q = q0+0.5*t*t*((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v;
225S60:
226/*
227 STEP 7: QUOTIENT ACCEPTANCE (Q)
228*/
229 if(log(1.0-u) <= q) return sgamma;
230S70:
231/*
232 STEP 8: E=STANDARD EXPONENTIAL DEVIATE
233 U= 0,1 -UNIFORM DEVIATE
234 T=(B,SI)-DOUBLE EXPONENTIAL (LAPLACE) SAMPLE
235*/
236 e = sexpo();
237 u = ranf();
238 u += (u-1.0);
239 t = b+fsign(si*e,u);
240/*
241 STEP 9: REJECTION IF T .LT. TAU(1) = -.71874483771719
242*/
243 if(t < -0.7187449) goto S70;
244/*
245 STEP 10: CALCULATION OF V AND QUOTIENT Q
246*/
247 v = t/(s+s);
248 if(fabs(v) <= 0.25) goto S80;
249 q = q0-s*t+0.25*t*t+(s2+s2)*log(1.0+v);
250 goto S90;
251S80:
252 q = q0+0.5*t*t*((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v;
253S90:
254/*
255 STEP 11: HAT ACCEPTANCE (H) (IF Q NOT POSITIVE GO TO STEP 8)
256*/
257 if(q <= 0.0) goto S70;
258 if(q <= 0.5) goto S100;
259/*
260 * JJV modified the code through line 115 to handle large Q case
261 */
262 if(q < 15.0) goto S95;
263/*
264 * JJV Here Q is large enough that Q = log(exp(Q) - 1.0) (for real Q)
265 * JJV so reformulate test at 110 in terms of one EXP, if not too big
266 * JJV 87.49823 is close to the largest real which can be
267 * JJV exponentiated (87.49823 = log(1.0E38))
268 */
269 if((q+e-0.5*t*t) > 87.49823) goto S115;
270 if(c*fabs(u) > exp(q+e-0.5*t*t)) goto S70;
271 goto S115;
272S95:
273 w = exp(q)-1.0;
274 goto S110;
275S100:
276 w = ((((e5*q+e4)*q+e3)*q+e2)*q+e1)*q;
277S110:
278/*
279 IF T IS REJECTED, SAMPLE AGAIN AT STEP 8
280*/
281 if(c*fabs(u) > w*exp(e-0.5*t*t)) goto S70;
282S115:
283 x = s+0.5*t;
284 sgamma = x*x;
285 return sgamma;
286S120:
287/*
288 ALTERNATE METHOD FOR PARAMETERS A BELOW 1 (.3678794=EXP(-1.))
289
290 JJV changed B to B0 (which was added to declarations for this)
291 JJV in 120 to END to fix rare and subtle bug.
292 JJV Line: 'aa = 0.0' was removed (unnecessary, wasteful).
293 JJV Reasons: the state of AA only serves to tell the A >= 1.0
294 JJV case if certain A-dependent constants need to be recalculated.
295 JJV The A < 1.0 case (here) no longer changes any of these, and
296 JJV the recalculation of B (which used to change with an
297 JJV A < 1.0 call) is governed by the state of AAA anyway.
298 aa = 0.0;
299*/
300 b0 = 1.0+0.3678794*a;
301S130:
302 p = b0*ranf();
303 if(p >= 1.0) goto S140;
304 sgamma = exp(log(p)/ a);
305 if(sexpo() < sgamma) goto S130;
306 return sgamma;
307S140:
308 sgamma = -log((b0-p)/ a);
309 if(sexpo() < (1.0-a)*log(sgamma)) goto S130;
310 return sgamma;
311
312}
313
314}
315
316#endif // BZ_RANDOM_GAMMA
Definition exponential.h:34
Definition gamma.h:46
ExponentialUnit< T, IRNG, sharedState > exponential_
Definition gamma.h:107
T fsign(T num, T sign)
Definition gamma.h:96
Gamma(T mean, unsigned int i)
Definition gamma.h:66
T sexpo()
Definition gamma.h:91
Gamma(T mean)
Definition gamma.h:61
T a
Definition gamma.h:109
NormalUnit< T, IRNG, sharedState > normal_
Definition gamma.h:106
Gamma()
Definition gamma.h:50
T ranf()
Definition gamma.h:81
T T_numtype
Definition gamma.h:48
T snorm()
Definition gamma.h:86
void setMean(T mean)
Definition gamma.h:74
Gamma(unsigned int i)
Definition gamma.h:55
T random()
Definition gamma.h:113
Definition normal.h:36
Definition uniform.h:377
T random()
Definition uniform.h:385
Definition beta.h:50
sharedState defaultState
Definition default.h:55
MersenneTwister defaultIRNG
Definition default.h:120