001/*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements.  See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License.  You may obtain a copy of the License at
008 *
009 *      http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017package org.apache.commons.rng.sampling.distribution;
018
019import org.apache.commons.rng.UniformRandomProvider;
020
021/**
022 * Sampling from the <a href="http://mathworld.wolfram.com/GammaDistribution.html">Gamma distribution</a>.
023 * <ul>
024 *  <li>
025 *   For {@code 0 < shape < 1}:
026 *   <blockquote>
027 *    Ahrens, J. H. and Dieter, U.,
028 *    <i>Computer methods for sampling from gamma, beta, Poisson and binomial distributions,</i>
029 *    Computing, 12, 223-246, 1974.
030 *   </blockquote>
031 *  </li>
032 *  <li>
033 *  For {@code shape >= 1}:
034 *   <blockquote>
035 *   Marsaglia and Tsang, <i>A Simple Method for Generating
036 *   Gamma Variables.</i> ACM Transactions on Mathematical Software,
037 *   Volume 26 Issue 3, September, 2000.
038 *   </blockquote>
039 *  </li>
040 * </ul>
041 */
042public class AhrensDieterMarsagliaTsangGammaSampler
043    extends SamplerBase
044    implements ContinuousSampler {
045    /** The shape parameter. */
046    private final double theta;
047    /** The alpha parameter. */
048    private final double alpha;
049    /** Gaussian sampling. */
050    private final BoxMullerGaussianSampler gaussian;
051
052    /**
053     * @param rng Generator of uniformly distributed random numbers.
054     * @param alpha Alpha parameter of the distribution.
055     * @param theta Theta parameter of the distribution.
056     */
057    public AhrensDieterMarsagliaTsangGammaSampler(UniformRandomProvider rng,
058                                                  double alpha,
059                                                  double theta) {
060        super(rng);
061        this.alpha = alpha;
062        this.theta = theta;
063        gaussian = new BoxMullerGaussianSampler(rng, 0, 1);
064    }
065
066    /** {@inheritDoc} */
067    @Override
068    public double sample() {
069        if (theta < 1) {
070            // [1]: p. 228, Algorithm GS.
071
072            while (true) {
073                // Step 1:
074                final double u = nextDouble();
075                final double bGS = 1 + theta / Math.E;
076                final double p = bGS * u;
077
078                if (p <= 1) {
079                    // Step 2:
080
081                    final double x = Math.pow(p, 1 / theta);
082                    final double u2 = nextDouble();
083
084                    if (u2 > Math.exp(-x)) {
085                        // Reject.
086                        continue;
087                    } else {
088                        return alpha * x;
089                    }
090                } else {
091                    // Step 3:
092
093                    final double x = -1 * Math.log((bGS - p) / theta);
094                    final double u2 = nextDouble();
095
096                    if (u2 > Math.pow(x, theta - 1)) {
097                        // Reject.
098                        continue;
099                    } else {
100                        return alpha * x;
101                    }
102                }
103            }
104        }
105
106        // Now theta >= 1.
107
108        final double d = theta - 0.333333333333333333;
109        final double c = 1 / (3 * Math.sqrt(d));
110
111        while (true) {
112            final double x = gaussian.sample();
113            final double v = (1 + c * x) * (1 + c * x) * (1 + c * x);
114
115            if (v <= 0) {
116                continue;
117            }
118
119            final double x2 = x * x;
120            final double u = nextDouble();
121
122            // Squeeze.
123            if (u < 1 - 0.0331 * x2 * x2) {
124                return alpha * d * v;
125            }
126
127            if (Math.log(u) < 0.5 * x2 + d * (1 - v + Math.log(v))) {
128                return alpha * d * v;
129            }
130        }
131    }
132
133    /** {@inheritDoc} */
134    @Override
135    public String toString() {
136        return "Ahrens-Dieter-Marsaglia-Tsang Gamma deviate [" + super.toString() + "]";
137    }
138}