1 | /* |
---|
2 | * Mathlib : A C Library of Special Functions |
---|
3 | * Copyright (C) 1998 Ross Ihaka |
---|
4 | * Copyright (C) 2000-2005 The R Development Core Team |
---|
5 | * |
---|
6 | * This program is free software; you can redistribute it and/or modify |
---|
7 | * it under the terms of the GNU General Public License as published by |
---|
8 | * the Free Software Foundation; either version 2 of the License, or |
---|
9 | * (at your option) any later version. |
---|
10 | * |
---|
11 | * This program is distributed in the hope that it will be useful, |
---|
12 | * but WITHOUT ANY WARRANTY; without even the implied warranty of |
---|
13 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
---|
14 | * GNU General Public License for more details. |
---|
15 | * |
---|
16 | * You should have received a copy of the GNU General Public License |
---|
17 | * along with this program; if not, a copy is available at |
---|
18 | * http://www.r-project.org/Licenses/ |
---|
19 | * |
---|
20 | * SYNOPSIS |
---|
21 | * |
---|
22 | * #include <Rmath.h> |
---|
23 | * double rgamma(double a, double scale); |
---|
24 | * |
---|
25 | * DESCRIPTION |
---|
26 | * |
---|
27 | * Random variates from the gamma distribution. |
---|
28 | * |
---|
29 | * REFERENCES |
---|
30 | * |
---|
31 | * [1] Shape parameter a >= 1. Algorithm GD in: |
---|
32 | * |
---|
33 | * Ahrens, J.H. and Dieter, U. (1982). |
---|
34 | * Generating gamma variates by a modified |
---|
35 | * rejection technique. |
---|
36 | * Comm. ACM, 25, 47-54. |
---|
37 | * |
---|
38 | * |
---|
39 | * [2] Shape parameter 0 < a < 1. Algorithm GS in: |
---|
40 | * |
---|
41 | * Ahrens, J.H. and Dieter, U. (1974). |
---|
42 | * Computer methods for sampling from gamma, beta, |
---|
43 | * poisson and binomial distributions. |
---|
44 | * Computing, 12, 223-246. |
---|
45 | * |
---|
46 | * Input: a = parameter (mean) of the standard gamma distribution. |
---|
47 | * Output: a variate from the gamma(a)-distribution |
---|
48 | */ |
---|
49 | |
---|
50 | #include "nmath.h" |
---|
51 | |
---|
52 | #define repeat for(;;) |
---|
53 | |
---|
54 | double rgamma(double a, double scale) |
---|
55 | { |
---|
56 | /* Constants : */ |
---|
57 | const static double sqrt32 = 5.656854; |
---|
58 | const static double exp_m1 = 0.36787944117144232159;/* exp(-1) = 1/e */ |
---|
59 | |
---|
60 | /* Coefficients q[k] - for q0 = sum(q[k]*a^(-k)) |
---|
61 | * Coefficients a[k] - for q = q0+(t*t/2)*sum(a[k]*v^k) |
---|
62 | * Coefficients e[k] - for exp(q)-1 = sum(e[k]*q^k) |
---|
63 | */ |
---|
64 | const static double q1 = 0.04166669; |
---|
65 | const static double q2 = 0.02083148; |
---|
66 | const static double q3 = 0.00801191; |
---|
67 | const static double q4 = 0.00144121; |
---|
68 | const static double q5 = -7.388e-5; |
---|
69 | const static double q6 = 2.4511e-4; |
---|
70 | const static double q7 = 2.424e-4; |
---|
71 | |
---|
72 | const static double a1 = 0.3333333; |
---|
73 | const static double a2 = -0.250003; |
---|
74 | const static double a3 = 0.2000062; |
---|
75 | const static double a4 = -0.1662921; |
---|
76 | const static double a5 = 0.1423657; |
---|
77 | const static double a6 = -0.1367177; |
---|
78 | const static double a7 = 0.1233795; |
---|
79 | |
---|
80 | /* State variables [FIXME for threading!] :*/ |
---|
81 | static double aa = 0.; |
---|
82 | static double aaa = 0.; |
---|
83 | static double s, s2, d; /* no. 1 (step 1) */ |
---|
84 | static double q0, b, si, c;/* no. 2 (step 4) */ |
---|
85 | |
---|
86 | double e, p, q, r, t, u, v, w, x, ret_val; |
---|
87 | |
---|
88 | if (!R_FINITE(a) || !R_FINITE(scale) || a < 0.0 || scale <= 0.0) |
---|
89 | ML_ERR_return_NAN; |
---|
90 | |
---|
91 | if (a < 1.) { /* GS algorithm for parameters a < 1 */ |
---|
92 | if(a == 0) |
---|
93 | return 0.; |
---|
94 | e = 1.0 + exp_m1 * a; |
---|
95 | repeat { |
---|
96 | p = e * unif_rand(); |
---|
97 | if (p >= 1.0) { |
---|
98 | x = -log((e - p) / a); |
---|
99 | if (exp_rand() >= (1.0 - a) * log(x)) |
---|
100 | break; |
---|
101 | } else { |
---|
102 | x = exp(log(p) / a); |
---|
103 | if (exp_rand() >= x) |
---|
104 | break; |
---|
105 | } |
---|
106 | } |
---|
107 | return scale * x; |
---|
108 | } |
---|
109 | |
---|
110 | /* --- a >= 1 : GD algorithm --- */ |
---|
111 | |
---|
112 | /* Step 1: Recalculations of s2, s, d if a has changed */ |
---|
113 | if (a != aa) { |
---|
114 | aa = a; |
---|
115 | s2 = a - 0.5; |
---|
116 | s = sqrt(s2); |
---|
117 | d = sqrt32 - s * 12.0; |
---|
118 | } |
---|
119 | /* Step 2: t = standard normal deviate, |
---|
120 | x = (s,1/2) -normal deviate. */ |
---|
121 | |
---|
122 | /* immediate acceptance (i) */ |
---|
123 | t = norm_rand(); |
---|
124 | x = s + 0.5 * t; |
---|
125 | ret_val = x * x; |
---|
126 | if (t >= 0.0) |
---|
127 | return scale * ret_val; |
---|
128 | |
---|
129 | /* Step 3: u = 0,1 - uniform sample. squeeze acceptance (s) */ |
---|
130 | u = unif_rand(); |
---|
131 | if (d * u <= t * t * t) |
---|
132 | return scale * ret_val; |
---|
133 | |
---|
134 | /* Step 4: recalculations of q0, b, si, c if necessary */ |
---|
135 | |
---|
136 | if (a != aaa) { |
---|
137 | aaa = a; |
---|
138 | r = 1.0 / a; |
---|
139 | q0 = ((((((q7 * r + q6) * r + q5) * r + q4) * r + q3) * r |
---|
140 | + q2) * r + q1) * r; |
---|
141 | |
---|
142 | /* Approximation depending on size of parameter a */ |
---|
143 | /* The constants in the expressions for b, si and c */ |
---|
144 | /* were established by numerical experiments */ |
---|
145 | |
---|
146 | if (a <= 3.686) { |
---|
147 | b = 0.463 + s + 0.178 * s2; |
---|
148 | si = 1.235; |
---|
149 | c = 0.195 / s - 0.079 + 0.16 * s; |
---|
150 | } else if (a <= 13.022) { |
---|
151 | b = 1.654 + 0.0076 * s2; |
---|
152 | si = 1.68 / s + 0.275; |
---|
153 | c = 0.062 / s + 0.024; |
---|
154 | } else { |
---|
155 | b = 1.77; |
---|
156 | si = 0.75; |
---|
157 | c = 0.1515 / s; |
---|
158 | } |
---|
159 | } |
---|
160 | /* Step 5: no quotient test if x not positive */ |
---|
161 | |
---|
162 | if (x > 0.0) { |
---|
163 | /* Step 6: calculation of v and quotient q */ |
---|
164 | v = t / (s + s); |
---|
165 | if (fabs(v) <= 0.25) |
---|
166 | q = q0 + 0.5 * t * t * ((((((a7 * v + a6) * v + a5) * v + a4) * v |
---|
167 | + a3) * v + a2) * v + a1) * v; |
---|
168 | else |
---|
169 | q = q0 - s * t + 0.25 * t * t + (s2 + s2) * log(1.0 + v); |
---|
170 | |
---|
171 | |
---|
172 | /* Step 7: quotient acceptance (q) */ |
---|
173 | if (log(1.0 - u) <= q) |
---|
174 | return scale * ret_val; |
---|
175 | } |
---|
176 | |
---|
177 | repeat { |
---|
178 | /* Step 8: e = standard exponential deviate |
---|
179 | * u = 0,1 -uniform deviate |
---|
180 | * t = (b,si)-double exponential (laplace) sample */ |
---|
181 | e = exp_rand(); |
---|
182 | u = unif_rand(); |
---|
183 | u = u + u - 1.0; |
---|
184 | if (u < 0.0) |
---|
185 | t = b - si * e; |
---|
186 | else |
---|
187 | t = b + si * e; |
---|
188 | /* Step 9: rejection if t < tau(1) = -0.71874483771719 */ |
---|
189 | if (t >= -0.71874483771719) { |
---|
190 | /* Step 10: calculation of v and quotient q */ |
---|
191 | v = t / (s + s); |
---|
192 | if (fabs(v) <= 0.25) |
---|
193 | q = q0 + 0.5 * t * t * |
---|
194 | ((((((a7 * v + a6) * v + a5) * v + a4) * v + a3) * v |
---|
195 | + a2) * v + a1) * v; |
---|
196 | else |
---|
197 | q = q0 - s * t + 0.25 * t * t + (s2 + s2) * log(1.0 + v); |
---|
198 | /* Step 11: hat acceptance (h) */ |
---|
199 | /* (if q not positive go to step 8) */ |
---|
200 | if (q > 0.0) { |
---|
201 | w = expm1(q); |
---|
202 | /* ^^^^^ original code had approximation with rel.err < 2e-7 */ |
---|
203 | /* if t is rejected sample again at step 8 */ |
---|
204 | if (c * fabs(u) <= w * exp(e - 0.5 * t * t)) |
---|
205 | break; |
---|
206 | } |
---|
207 | } |
---|
208 | } /* repeat .. until `t' is accepted */ |
---|
209 | x = s + 0.5 * t; |
---|
210 | return scale * x * x; |
---|
211 | } |
---|