1 /*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17 package org.apache.commons.rng.sampling.distribution;
18
19 import org.apache.commons.rng.UniformRandomProvider;
20
21 /**
22 * Sampling from the <a href="http://mathworld.wolfram.com/GammaDistribution.html">gamma distribution</a>.
23 * <ul>
24 * <li>
25 * For {@code 0 < alpha < 1}:
26 * <blockquote>
27 * Ahrens, J. H. and Dieter, U.,
28 * <i>Computer methods for sampling from gamma, beta, Poisson and binomial distributions,</i>
29 * Computing, 12, 223-246, 1974.
30 * </blockquote>
31 * </li>
32 * <li>
33 * For {@code alpha >= 1}:
34 * <blockquote>
35 * Marsaglia and Tsang, <i>A Simple Method for Generating
36 * Gamma Variables.</i> ACM Transactions on Mathematical Software,
37 * Volume 26 Issue 3, September, 2000.
38 * </blockquote>
39 * </li>
40 * </ul>
41 *
42 * <p>Sampling uses:</p>
43 *
44 * <ul>
45 * <li>{@link UniformRandomProvider#nextDouble()} (both algorithms)
46 * <li>{@link UniformRandomProvider#nextLong()} (only for {@code alpha >= 1})
47 * </ul>
48 *
49 * @see <a href="http://mathworld.wolfram.com/GammaDistribution.html">MathWorld Gamma distribution</a>
50 * @see <a href="https://en.wikipedia.org/wiki/Gamma_distribution">Wikipedia Gamma distribution</a>
51 * @since 1.0
52 */
53 public class AhrensDieterMarsagliaTsangGammaSampler
54 extends SamplerBase
55 implements SharedStateContinuousSampler {
56 /** The appropriate gamma sampler for the parameters. */
57 private final SharedStateContinuousSampler delegate;
58
59 /**
60 * Base class for a sampler from the Gamma distribution.
61 */
62 private abstract static class BaseGammaSampler
63 implements SharedStateContinuousSampler {
64
65 /** Underlying source of randomness. */
66 protected final UniformRandomProvider rng;
67 /** The alpha parameter. This is a shape parameter. */
68 protected final double alpha;
69 /** The theta parameter. This is a scale parameter. */
70 protected final double theta;
71
72 /**
73 * @param rng Generator of uniformly distributed random numbers.
74 * @param alpha Alpha parameter of the distribution.
75 * @param theta Theta parameter of the distribution.
76 * @throws IllegalArgumentException if {@code alpha <= 0} or {@code theta <= 0}
77 */
78 BaseGammaSampler(UniformRandomProvider rng,
79 double alpha,
80 double theta) {
81 // Validation before java.lang.Object constructor exits prevents partially initialized object
82 this(InternalUtils.requireStrictlyPositive(alpha, "alpha"),
83 InternalUtils.requireStrictlyPositive(theta, "theta"),
84 rng);
85 }
86
87 /**
88 * @param alpha Alpha parameter of the distribution.
89 * @param theta Theta parameter of the distribution.
90 * @param rng Generator of uniformly distributed random numbers.
91 */
92 private BaseGammaSampler(double alpha,
93 double theta,
94 UniformRandomProvider rng) {
95 this.rng = rng;
96 this.alpha = alpha;
97 this.theta = theta;
98 }
99
100 /**
101 * @param rng Generator of uniformly distributed random numbers.
102 * @param source Source to copy.
103 */
104 BaseGammaSampler(UniformRandomProvider rng,
105 BaseGammaSampler source) {
106 this.rng = rng;
107 this.alpha = source.alpha;
108 this.theta = source.theta;
109 }
110
111 /** {@inheritDoc} */
112 @Override
113 public String toString() {
114 return "Ahrens-Dieter-Marsaglia-Tsang Gamma deviate [" + rng.toString() + "]";
115 }
116 }
117
118 /**
119 * Class to sample from the Gamma distribution when {@code 0 < alpha < 1}.
120 *
121 * <blockquote>
122 * Ahrens, J. H. and Dieter, U.,
123 * <i>Computer methods for sampling from gamma, beta, Poisson and binomial distributions,</i>
124 * Computing, 12, 223-246, 1974.
125 * </blockquote>
126 */
127 private static final class AhrensDieterGammaSampler
128 extends BaseGammaSampler {
129
130 /** Inverse of "alpha". */
131 private final double oneOverAlpha;
132 /** Optimization (see code). */
133 private final double bGSOptim;
134
135 /**
136 * @param rng Generator of uniformly distributed random numbers.
137 * @param alpha Alpha parameter of the distribution.
138 * @param theta Theta parameter of the distribution.
139 * @throws IllegalArgumentException if {@code alpha <= 0} or {@code theta <= 0}
140 */
141 AhrensDieterGammaSampler(UniformRandomProvider rng,
142 double alpha,
143 double theta) {
144 super(rng, alpha, theta);
145 oneOverAlpha = 1 / alpha;
146 bGSOptim = 1 + alpha / Math.E;
147 }
148
149 /**
150 * @param rng Generator of uniformly distributed random numbers.
151 * @param source Source to copy.
152 */
153 AhrensDieterGammaSampler(UniformRandomProvider rng,
154 AhrensDieterGammaSampler source) {
155 super(rng, source);
156 oneOverAlpha = source.oneOverAlpha;
157 bGSOptim = source.bGSOptim;
158 }
159
160 @Override
161 public double sample() {
162 // [1]: p. 228, Algorithm GS.
163
164 while (true) {
165 // Step 1:
166 final double u = rng.nextDouble();
167 final double p = bGSOptim * u;
168
169 if (p <= 1) {
170 // Step 2:
171 final double x = Math.pow(p, oneOverAlpha);
172 final double u2 = rng.nextDouble();
173
174 if (u2 > Math.exp(-x)) {
175 // Reject.
176 continue;
177 }
178 return theta * x;
179 }
180
181 // Step 3:
182 final double x = -Math.log((bGSOptim - p) * oneOverAlpha);
183 final double u2 = rng.nextDouble();
184
185 if (u2 <= Math.pow(x, alpha - 1)) {
186 return theta * x;
187 }
188 // Reject and continue.
189 }
190 }
191
192 @Override
193 public SharedStateContinuousSampler withUniformRandomProvider(UniformRandomProvider rng) {
194 return new AhrensDieterGammaSampler(rng, this);
195 }
196 }
197
198 /**
199 * Class to sample from the Gamma distribution when the {@code alpha >= 1}.
200 *
201 * <blockquote>
202 * Marsaglia and Tsang, <i>A Simple Method for Generating
203 * Gamma Variables.</i> ACM Transactions on Mathematical Software,
204 * Volume 26 Issue 3, September, 2000.
205 * </blockquote>
206 */
207 private static final class MarsagliaTsangGammaSampler
208 extends BaseGammaSampler {
209
210 /** 1/3. */
211 private static final double ONE_THIRD = 1d / 3;
212
213 /** Optimization (see code). */
214 private final double dOptim;
215 /** Optimization (see code). */
216 private final double cOptim;
217 /** Gaussian sampling. */
218 private final NormalizedGaussianSampler gaussian;
219
220 /**
221 * @param rng Generator of uniformly distributed random numbers.
222 * @param alpha Alpha parameter of the distribution.
223 * @param theta Theta parameter of the distribution.
224 * @throws IllegalArgumentException if {@code alpha <= 0} or {@code theta <= 0}
225 */
226 MarsagliaTsangGammaSampler(UniformRandomProvider rng,
227 double alpha,
228 double theta) {
229 super(rng, alpha, theta);
230 gaussian = ZigguratSampler.NormalizedGaussian.of(rng);
231 dOptim = alpha - ONE_THIRD;
232 cOptim = ONE_THIRD / Math.sqrt(dOptim);
233 }
234
235 /**
236 * @param rng Generator of uniformly distributed random numbers.
237 * @param source Source to copy.
238 */
239 MarsagliaTsangGammaSampler(UniformRandomProvider rng,
240 MarsagliaTsangGammaSampler source) {
241 super(rng, source);
242 gaussian = ZigguratSampler.NormalizedGaussian.of(rng);
243 dOptim = source.dOptim;
244 cOptim = source.cOptim;
245 }
246
247 @Override
248 public double sample() {
249 while (true) {
250 final double x = gaussian.sample();
251 final double oPcTx = 1 + cOptim * x;
252 final double v = oPcTx * oPcTx * oPcTx;
253
254 if (v <= 0) {
255 continue;
256 }
257
258 final double x2 = x * x;
259 final double u = rng.nextDouble();
260
261 // Squeeze.
262 if (u < 1 - 0.0331 * x2 * x2) {
263 return theta * dOptim * v;
264 }
265
266 if (Math.log(u) < 0.5 * x2 + dOptim * (1 - v + Math.log(v))) {
267 return theta * dOptim * v;
268 }
269 }
270 }
271
272 @Override
273 public SharedStateContinuousSampler withUniformRandomProvider(UniformRandomProvider rng) {
274 return new MarsagliaTsangGammaSampler(rng, this);
275 }
276 }
277
278 /**
279 * This instance delegates sampling. Use the factory method
280 * {@link #of(UniformRandomProvider, double, double)} to create an optimal sampler.
281 *
282 * @param rng Generator of uniformly distributed random numbers.
283 * @param alpha Alpha parameter of the distribution (this is a shape parameter).
284 * @param theta Theta parameter of the distribution (this is a scale parameter).
285 * @throws IllegalArgumentException if {@code alpha <= 0} or {@code theta <= 0}
286 */
287 public AhrensDieterMarsagliaTsangGammaSampler(UniformRandomProvider rng,
288 double alpha,
289 double theta) {
290 super(null);
291 delegate = of(rng, alpha, theta);
292 }
293
294 /** {@inheritDoc} */
295 @Override
296 public double sample() {
297 return delegate.sample();
298 }
299
300 /** {@inheritDoc} */
301 @Override
302 public String toString() {
303 return delegate.toString();
304 }
305
306 /**
307 * {@inheritDoc}
308 *
309 * @since 1.3
310 */
311 @Override
312 public SharedStateContinuousSampler withUniformRandomProvider(UniformRandomProvider rng) {
313 // Direct return of the optimised sampler
314 return delegate.withUniformRandomProvider(rng);
315 }
316
317 /**
318 * Creates a new gamma distribution sampler.
319 *
320 * @param rng Generator of uniformly distributed random numbers.
321 * @param alpha Alpha parameter of the distribution (this is a shape parameter).
322 * @param theta Theta parameter of the distribution (this is a scale parameter).
323 * @return the sampler
324 * @throws IllegalArgumentException if {@code alpha <= 0} or {@code theta <= 0}
325 * @since 1.3
326 */
327 public static SharedStateContinuousSampler of(UniformRandomProvider rng,
328 double alpha,
329 double theta) {
330 // Each sampler should check the input arguments.
331 return alpha < 1 ?
332 new AhrensDieterGammaSampler(rng, alpha, theta) :
333 new MarsagliaTsangGammaSampler(rng, alpha, theta);
334 }
335 }