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