1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.apache.commons.rng.sampling.distribution;
18
19 import org.apache.commons.rng.UniformRandomProvider;
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53 public class AhrensDieterMarsagliaTsangGammaSampler
54 extends SamplerBase
55 implements SharedStateContinuousSampler {
56
57 private final SharedStateContinuousSampler delegate;
58
59
60
61
62 private abstract static class BaseGammaSampler
63 implements SharedStateContinuousSampler {
64
65
66 protected final UniformRandomProvider rng;
67
68 protected final double alpha;
69
70 protected final double theta;
71
72
73
74
75
76
77
78 BaseGammaSampler(UniformRandomProvider rng,
79 double alpha,
80 double theta) {
81 if (alpha <= 0) {
82 throw new IllegalArgumentException("alpha is not strictly positive: " + alpha);
83 }
84 if (theta <= 0) {
85 throw new IllegalArgumentException("theta is not strictly positive: " + theta);
86 }
87 this.rng = rng;
88 this.alpha = alpha;
89 this.theta = theta;
90 }
91
92
93
94
95
96 BaseGammaSampler(UniformRandomProvider rng,
97 BaseGammaSampler source) {
98 this.rng = rng;
99 this.alpha = source.alpha;
100 this.theta = source.theta;
101 }
102
103
104 @Override
105 public String toString() {
106 return "Ahrens-Dieter-Marsaglia-Tsang Gamma deviate [" + rng.toString() + "]";
107 }
108 }
109
110
111
112
113
114
115
116
117
118
119 private static class AhrensDieterGammaSampler
120 extends BaseGammaSampler {
121
122
123 private final double oneOverAlpha;
124
125 private final double bGSOptim;
126
127
128
129
130
131
132
133 AhrensDieterGammaSampler(UniformRandomProvider rng,
134 double alpha,
135 double theta) {
136 super(rng, alpha, theta);
137 oneOverAlpha = 1 / alpha;
138 bGSOptim = 1 + alpha / Math.E;
139 }
140
141
142
143
144
145 AhrensDieterGammaSampler(UniformRandomProvider rng,
146 AhrensDieterGammaSampler source) {
147 super(rng, source);
148 oneOverAlpha = source.oneOverAlpha;
149 bGSOptim = source.bGSOptim;
150 }
151
152 @Override
153 public double sample() {
154
155
156 while (true) {
157
158 final double u = rng.nextDouble();
159 final double p = bGSOptim * u;
160
161 if (p <= 1) {
162
163 final double x = Math.pow(p, oneOverAlpha);
164 final double u2 = rng.nextDouble();
165
166 if (u2 > Math.exp(-x)) {
167
168 continue;
169 }
170 return theta * x;
171 }
172
173
174 final double x = -Math.log((bGSOptim - p) * oneOverAlpha);
175 final double u2 = rng.nextDouble();
176
177 if (u2 <= Math.pow(x, alpha - 1)) {
178 return theta * x;
179 }
180
181 }
182 }
183
184 @Override
185 public SharedStateContinuousSampler withUniformRandomProvider(UniformRandomProvider rng) {
186 return new AhrensDieterGammaSampler(rng, this);
187 }
188 }
189
190
191
192
193
194
195
196
197
198
199 private static class MarsagliaTsangGammaSampler
200 extends BaseGammaSampler {
201
202
203 private static final double ONE_THIRD = 1d / 3;
204
205
206 private final double dOptim;
207
208 private final double cOptim;
209
210 private final NormalizedGaussianSampler gaussian;
211
212
213
214
215
216
217
218 MarsagliaTsangGammaSampler(UniformRandomProvider rng,
219 double alpha,
220 double theta) {
221 super(rng, alpha, theta);
222 gaussian = ZigguratSampler.NormalizedGaussian.of(rng);
223 dOptim = alpha - ONE_THIRD;
224 cOptim = ONE_THIRD / Math.sqrt(dOptim);
225 }
226
227
228
229
230
231 MarsagliaTsangGammaSampler(UniformRandomProvider rng,
232 MarsagliaTsangGammaSampler source) {
233 super(rng, source);
234 gaussian = ZigguratSampler.NormalizedGaussian.of(rng);
235 dOptim = source.dOptim;
236 cOptim = source.cOptim;
237 }
238
239 @Override
240 public double sample() {
241 while (true) {
242 final double x = gaussian.sample();
243 final double oPcTx = 1 + cOptim * x;
244 final double v = oPcTx * oPcTx * oPcTx;
245
246 if (v <= 0) {
247 continue;
248 }
249
250 final double x2 = x * x;
251 final double u = rng.nextDouble();
252
253
254 if (u < 1 - 0.0331 * x2 * x2) {
255 return theta * dOptim * v;
256 }
257
258 if (Math.log(u) < 0.5 * x2 + dOptim * (1 - v + Math.log(v))) {
259 return theta * dOptim * v;
260 }
261 }
262 }
263
264 @Override
265 public SharedStateContinuousSampler withUniformRandomProvider(UniformRandomProvider rng) {
266 return new MarsagliaTsangGammaSampler(rng, this);
267 }
268 }
269
270
271
272
273
274
275
276
277
278
279 public AhrensDieterMarsagliaTsangGammaSampler(UniformRandomProvider rng,
280 double alpha,
281 double theta) {
282 super(null);
283 delegate = of(rng, alpha, theta);
284 }
285
286
287 @Override
288 public double sample() {
289 return delegate.sample();
290 }
291
292
293 @Override
294 public String toString() {
295 return delegate.toString();
296 }
297
298
299
300
301
302
303 @Override
304 public SharedStateContinuousSampler withUniformRandomProvider(UniformRandomProvider rng) {
305
306 return delegate.withUniformRandomProvider(rng);
307 }
308
309
310
311
312
313
314
315
316
317
318
319 public static SharedStateContinuousSampler of(UniformRandomProvider rng,
320 double alpha,
321 double theta) {
322
323 return alpha < 1 ?
324 new AhrensDieterGammaSampler(rng, alpha, theta) :
325 new MarsagliaTsangGammaSampler(rng, alpha, theta);
326 }
327 }