[29] | 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 | } |
---|