1 | // |
---|
2 | // C++ Implementation: itpp_ext |
---|
3 | // |
---|
4 | // Description: |
---|
5 | // |
---|
6 | // |
---|
7 | // Author: smidl <smidl@utia.cas.cz>, (C) 2008 |
---|
8 | // |
---|
9 | // Copyright: See COPYING file that comes with this distribution |
---|
10 | // |
---|
11 | // |
---|
12 | |
---|
13 | #include <itpp/itbase.h> |
---|
14 | #include "itpp_ext.h" |
---|
15 | |
---|
16 | namespace itpp { |
---|
17 | Array<int> to_Arr(const ivec &indices) |
---|
18 | { |
---|
19 | Array<int> a(indices.size()); |
---|
20 | for (int i = 0; i < a.size(); i++) { |
---|
21 | a(i) = indices(i); |
---|
22 | } |
---|
23 | return a; |
---|
24 | } |
---|
25 | |
---|
26 | //Gamma |
---|
27 | |
---|
28 | Gamma_RNG::Gamma_RNG(double a, double b) |
---|
29 | { |
---|
30 | setup(a,b); |
---|
31 | } |
---|
32 | |
---|
33 | #define log std::log |
---|
34 | #define exp std::exp |
---|
35 | #define sqrt std::sqrt |
---|
36 | #define R_FINITE std::isfinite |
---|
37 | |
---|
38 | double Gamma_RNG::sample() |
---|
39 | { |
---|
40 | //A copy of rgamma code from the R package!! |
---|
41 | // |
---|
42 | |
---|
43 | /* Constants : */ |
---|
44 | const static double sqrt32 = 5.656854; |
---|
45 | const static double exp_m1 = 0.36787944117144232159;/* exp(-1) = 1/e */ |
---|
46 | |
---|
47 | /* Coefficients q[k] - for q0 = sum(q[k]*a^(-k)) |
---|
48 | * Coefficients a[k] - for q = q0+(t*t/2)*sum(a[k]*v^k) |
---|
49 | * Coefficients e[k] - for exp(q)-1 = sum(e[k]*q^k) |
---|
50 | */ |
---|
51 | const static double q1 = 0.04166669; |
---|
52 | const static double q2 = 0.02083148; |
---|
53 | const static double q3 = 0.00801191; |
---|
54 | const static double q4 = 0.00144121; |
---|
55 | const static double q5 = -7.388e-5; |
---|
56 | const static double q6 = 2.4511e-4; |
---|
57 | const static double q7 = 2.424e-4; |
---|
58 | |
---|
59 | const static double a1 = 0.3333333; |
---|
60 | const static double a2 = -0.250003; |
---|
61 | const static double a3 = 0.2000062; |
---|
62 | const static double a4 = -0.1662921; |
---|
63 | const static double a5 = 0.1423657; |
---|
64 | const static double a6 = -0.1367177; |
---|
65 | const static double a7 = 0.1233795; |
---|
66 | |
---|
67 | /* State variables [FIXME for threading!] :*/ |
---|
68 | static double aa = 0.; |
---|
69 | static double aaa = 0.; |
---|
70 | static double s, s2, d; /* no. 1 (step 1) */ |
---|
71 | static double q0, b, si, c;/* no. 2 (step 4) */ |
---|
72 | |
---|
73 | double e, p, q, r, t, u, v, w, x, ret_val; |
---|
74 | double a=alpha; |
---|
75 | double scale=1.0/beta; |
---|
76 | |
---|
77 | if (!R_FINITE(a) || !R_FINITE(scale) || a < 0.0 || scale <= 0.0) |
---|
78 | it_error("Gamma_RNG wrong parameters"); |
---|
79 | |
---|
80 | if (a < 1.) { /* GS algorithm for parameters a < 1 */ |
---|
81 | if(a == 0) |
---|
82 | return 0.; |
---|
83 | e = 1.0 + exp_m1 * a; |
---|
84 | for(;;) { //VS repeat |
---|
85 | p = e * unif_rand(); |
---|
86 | if (p >= 1.0) { |
---|
87 | x = -log((e - p) / a); |
---|
88 | if (exp_rand() >= (1.0 - a) * log(x)) |
---|
89 | break; |
---|
90 | } else { |
---|
91 | x = exp(log(p) / a); |
---|
92 | if (exp_rand() >= x) |
---|
93 | break; |
---|
94 | } |
---|
95 | } |
---|
96 | return scale * x; |
---|
97 | } |
---|
98 | |
---|
99 | /* --- a >= 1 : GD algorithm --- */ |
---|
100 | |
---|
101 | /* Step 1: Recalculations of s2, s, d if a has changed */ |
---|
102 | if (a != aa) { |
---|
103 | aa = a; |
---|
104 | s2 = a - 0.5; |
---|
105 | s = sqrt(s2); |
---|
106 | d = sqrt32 - s * 12.0; |
---|
107 | } |
---|
108 | /* Step 2: t = standard normal deviate, |
---|
109 | x = (s,1/2) -normal deviate. */ |
---|
110 | |
---|
111 | /* immediate acceptance (i) */ |
---|
112 | t = norm_rand(); |
---|
113 | x = s + 0.5 * t; |
---|
114 | ret_val = x * x; |
---|
115 | if (t >= 0.0) |
---|
116 | return scale * ret_val; |
---|
117 | |
---|
118 | /* Step 3: u = 0,1 - uniform sample. squeeze acceptance (s) */ |
---|
119 | u = unif_rand(); |
---|
120 | if ((d * u) <= (t * t * t)) |
---|
121 | return scale * ret_val; |
---|
122 | |
---|
123 | /* Step 4: recalculations of q0, b, si, c if necessary */ |
---|
124 | |
---|
125 | if (a != aaa) { |
---|
126 | aaa = a; |
---|
127 | r = 1.0 / a; |
---|
128 | q0 = ((((((q7 * r + q6) * r + q5) * r + q4) * r + q3) * r |
---|
129 | + q2) * r + q1) * r; |
---|
130 | |
---|
131 | /* Approximation depending on size of parameter a */ |
---|
132 | /* The constants in the expressions for b, si and c */ |
---|
133 | /* were established by numerical experiments */ |
---|
134 | |
---|
135 | if (a <= 3.686) { |
---|
136 | b = 0.463 + s + 0.178 * s2; |
---|
137 | si = 1.235; |
---|
138 | c = 0.195 / s - 0.079 + 0.16 * s; |
---|
139 | } else if (a <= 13.022) { |
---|
140 | b = 1.654 + 0.0076 * s2; |
---|
141 | si = 1.68 / s + 0.275; |
---|
142 | c = 0.062 / s + 0.024; |
---|
143 | } else { |
---|
144 | b = 1.77; |
---|
145 | si = 0.75; |
---|
146 | c = 0.1515 / s; |
---|
147 | } |
---|
148 | } |
---|
149 | /* Step 5: no quotient test if x not positive */ |
---|
150 | |
---|
151 | if (x > 0.0) { |
---|
152 | /* Step 6: calculation of v and quotient q */ |
---|
153 | v = t / (s + s); |
---|
154 | if (fabs(v) <= 0.25) |
---|
155 | q = q0 + 0.5 * t * t * ((((((a7 * v + a6) * v + a5) * v + a4) * v |
---|
156 | + a3) * v + a2) * v + a1) * v; |
---|
157 | else |
---|
158 | q = q0 - s * t + 0.25 * t * t + (s2 + s2) * log(1.0 + v); |
---|
159 | |
---|
160 | |
---|
161 | /* Step 7: quotient acceptance (q) */ |
---|
162 | if (log(1.0 - u) <= q) |
---|
163 | return scale * ret_val; |
---|
164 | } |
---|
165 | |
---|
166 | for(;;) { //VS repeat |
---|
167 | /* Step 8: e = standard exponential deviate |
---|
168 | * u = 0,1 -uniform deviate |
---|
169 | * t = (b,si)-double exponential (laplace) sample */ |
---|
170 | e = exp_rand(); |
---|
171 | u = unif_rand(); |
---|
172 | u = u + u - 1.0; |
---|
173 | if (u < 0.0) |
---|
174 | t = b - si * e; |
---|
175 | else |
---|
176 | t = b + si * e; |
---|
177 | /* Step 9: rejection if t < tau(1) = -0.71874483771719 */ |
---|
178 | if (t >= -0.71874483771719) { |
---|
179 | /* Step 10: calculation of v and quotient q */ |
---|
180 | v = t / (s + s); |
---|
181 | if (fabs(v) <= 0.25) |
---|
182 | q = q0 + 0.5 * t * t * |
---|
183 | ((((((a7 * v + a6) * v + a5) * v + a4) * v + a3) * v |
---|
184 | + a2) * v + a1) * v; |
---|
185 | else |
---|
186 | q = q0 - s * t + 0.25 * t * t + (s2 + s2) * log(1.0 + v); |
---|
187 | /* Step 11: hat acceptance (h) */ |
---|
188 | /* (if q not positive go to step 8) */ |
---|
189 | if (q > 0.0) { |
---|
190 | w = expm1(q); |
---|
191 | /* ^^^^^ original code had approximation with rel.err < 2e-7 */ |
---|
192 | /* if t is rejected sample again at step 8 */ |
---|
193 | if ((c * fabs(u)) <= (w * exp(e - 0.5 * t * t))) |
---|
194 | break; |
---|
195 | } |
---|
196 | } |
---|
197 | } /* repeat .. until `t' is accepted */ |
---|
198 | x = s + 0.5 * t; |
---|
199 | return scale * x * x; |
---|
200 | } |
---|
201 | |
---|
202 | } |
---|