View Javadoc
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              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           * @param rng Generator of uniformly distributed random numbers.
94           * @param source Source to copy.
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         /** {@inheritDoc} */
104         @Override
105         public String toString() {
106             return "Ahrens-Dieter-Marsaglia-Tsang Gamma deviate [" + rng.toString() + "]";
107         }
108     }
109 
110     /**
111      * Class to sample from the Gamma distribution when {@code 0 < alpha < 1}.
112      *
113      * <blockquote>
114      *  Ahrens, J. H. and Dieter, U.,
115      *  <i>Computer methods for sampling from gamma, beta, Poisson and binomial distributions,</i>
116      *  Computing, 12, 223-246, 1974.
117      * </blockquote>
118      */
119     private static class AhrensDieterGammaSampler
120         extends BaseGammaSampler {
121 
122         /** Inverse of "alpha". */
123         private final double oneOverAlpha;
124         /** Optimization (see code). */
125         private final double bGSOptim;
126 
127         /**
128          * @param rng Generator of uniformly distributed random numbers.
129          * @param alpha Alpha parameter of the distribution.
130          * @param theta Theta parameter of the distribution.
131          * @throws IllegalArgumentException if {@code alpha <= 0} or {@code theta <= 0}
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          * @param rng Generator of uniformly distributed random numbers.
143          * @param source Source to copy.
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             // [1]: p. 228, Algorithm GS.
155 
156             while (true) {
157                 // Step 1:
158                 final double u = rng.nextDouble();
159                 final double p = bGSOptim * u;
160 
161                 if (p <= 1) {
162                     // Step 2:
163                     final double x = Math.pow(p, oneOverAlpha);
164                     final double u2 = rng.nextDouble();
165 
166                     if (u2 > Math.exp(-x)) {
167                         // Reject.
168                         continue;
169                     }
170                     return theta * x;
171                 }
172 
173                 // Step 3:
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                 // Reject and continue.
181             }
182         }
183 
184         @Override
185         public SharedStateContinuousSampler withUniformRandomProvider(UniformRandomProvider rng) {
186             return new AhrensDieterGammaSampler(rng, this);
187         }
188     }
189 
190     /**
191      * Class to sample from the Gamma distribution when the {@code alpha >= 1}.
192      *
193      * <blockquote>
194      *  Marsaglia and Tsang, <i>A Simple Method for Generating
195      *  Gamma Variables.</i> ACM Transactions on Mathematical Software,
196      *  Volume 26 Issue 3, September, 2000.
197      * </blockquote>
198      */
199     private static class MarsagliaTsangGammaSampler
200         extends BaseGammaSampler {
201 
202         /** 1/3. */
203         private static final double ONE_THIRD = 1d / 3;
204 
205         /** Optimization (see code). */
206         private final double dOptim;
207         /** Optimization (see code). */
208         private final double cOptim;
209         /** Gaussian sampling. */
210         private final NormalizedGaussianSampler gaussian;
211 
212         /**
213          * @param rng Generator of uniformly distributed random numbers.
214          * @param alpha Alpha parameter of the distribution.
215          * @param theta Theta parameter of the distribution.
216          * @throws IllegalArgumentException if {@code alpha <= 0} or {@code theta <= 0}
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          * @param rng Generator of uniformly distributed random numbers.
229          * @param source Source to copy.
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                 // Squeeze.
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      * This instance delegates sampling. Use the factory method
272      * {@link #of(UniformRandomProvider, double, double)} to create an optimal sampler.
273      *
274      * @param rng Generator of uniformly distributed random numbers.
275      * @param alpha Alpha parameter of the distribution (this is a shape parameter).
276      * @param theta Theta parameter of the distribution (this is a scale parameter).
277      * @throws IllegalArgumentException if {@code alpha <= 0} or {@code theta <= 0}
278      */
279     public AhrensDieterMarsagliaTsangGammaSampler(UniformRandomProvider rng,
280                                                   double alpha,
281                                                   double theta) {
282         super(null);
283         delegate = of(rng, alpha, theta);
284     }
285 
286     /** {@inheritDoc} */
287     @Override
288     public double sample() {
289         return delegate.sample();
290     }
291 
292     /** {@inheritDoc} */
293     @Override
294     public String toString() {
295         return delegate.toString();
296     }
297 
298     /**
299      * {@inheritDoc}
300      *
301      * @since 1.3
302      */
303     @Override
304     public SharedStateContinuousSampler withUniformRandomProvider(UniformRandomProvider rng) {
305         // Direct return of the optimised sampler
306         return delegate.withUniformRandomProvider(rng);
307     }
308 
309     /**
310      * Creates a new gamma distribution sampler.
311      *
312      * @param rng Generator of uniformly distributed random numbers.
313      * @param alpha Alpha parameter of the distribution (this is a shape parameter).
314      * @param theta Theta parameter of the distribution (this is a scale parameter).
315      * @return the sampler
316      * @throws IllegalArgumentException if {@code alpha <= 0} or {@code theta <= 0}
317      * @since 1.3
318      */
319     public static SharedStateContinuousSampler of(UniformRandomProvider rng,
320                                                   double alpha,
321                                                   double theta) {
322         // Each sampler should check the input arguments.
323         return alpha < 1 ?
324                 new AhrensDieterGammaSampler(rng, alpha, theta) :
325                 new MarsagliaTsangGammaSampler(rng, alpha, theta);
326     }
327 }