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 < theta < 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 theta >= 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   * @since 1.0
43   */
44  public class AhrensDieterMarsagliaTsangGammaSampler
45      extends SamplerBase
46      implements ContinuousSampler {
47      /** 1/3 */
48      private static final double ONE_THIRD = 1d / 3;
49      /** The shape parameter. */
50      private final double theta;
51      /** The alpha parameter. */
52      private final double alpha;
53      /** Inverse of "theta". */
54      private final double oneOverTheta;
55      /** Optimization (see code). */
56      private final double bGSOptim;
57      /** Optimization (see code). */
58      private final double dOptim;
59      /** Optimization (see code). */
60      private final double cOptim;
61      /** Gaussian sampling. */
62      private final NormalizedGaussianSampler gaussian;
63      /** Underlying source of randomness. */
64      private final UniformRandomProvider rng;
65  
66      /**
67       * @param rng Generator of uniformly distributed random numbers.
68       * @param alpha Alpha parameter of the distribution.
69       * @param theta Theta parameter of the distribution.
70       */
71      public AhrensDieterMarsagliaTsangGammaSampler(UniformRandomProvider rng,
72                                                    double alpha,
73                                                    double theta) {
74          super(null);
75          this.rng = rng;
76          this.alpha = alpha;
77          this.theta = theta;
78          gaussian = new ZigguratNormalizedGaussianSampler(rng);
79          oneOverTheta = 1 / theta;
80          bGSOptim = 1 + theta / Math.E;
81          dOptim = theta - ONE_THIRD;
82          cOptim = ONE_THIRD / Math.sqrt(dOptim);
83      }
84  
85      /** {@inheritDoc} */
86      @Override
87      public double sample() {
88          if (theta < 1) {
89              // [1]: p. 228, Algorithm GS.
90  
91              while (true) {
92                  // Step 1:
93                  final double u = rng.nextDouble();
94                  final double p = bGSOptim * u;
95  
96                  if (p <= 1) {
97                      // Step 2:
98  
99                      final double x = Math.pow(p, oneOverTheta);
100                     final double u2 = rng.nextDouble();
101 
102                     if (u2 > Math.exp(-x)) {
103                         // Reject.
104                         continue;
105                     } else {
106                         return alpha * x;
107                     }
108                 } else {
109                     // Step 3:
110 
111                     final double x = -Math.log((bGSOptim - p) * oneOverTheta);
112                     final double u2 = rng.nextDouble();
113 
114                     if (u2 > Math.pow(x, theta - 1)) {
115                         // Reject.
116                         continue;
117                     } else {
118                         return alpha * x;
119                     }
120                 }
121             }
122         } else {
123             while (true) {
124                 final double x = gaussian.sample();
125                 final double oPcTx = 1 + cOptim * x;
126                 final double v = oPcTx * oPcTx * oPcTx;
127 
128                 if (v <= 0) {
129                     continue;
130                 }
131 
132                 final double x2 = x * x;
133                 final double u = rng.nextDouble();
134 
135                 // Squeeze.
136                 if (u < 1 - 0.0331 * x2 * x2) {
137                     return alpha * dOptim * v;
138                 }
139 
140                 if (Math.log(u) < 0.5 * x2 + dOptim * (1 - v + Math.log(v))) {
141                     return alpha * dOptim * v;
142                 }
143             }
144         }
145     }
146 
147     /** {@inheritDoc} */
148     @Override
149     public String toString() {
150         return "Ahrens-Dieter-Marsaglia-Tsang Gamma deviate [" + rng.toString() + "]";
151     }
152 }