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 < theta < 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 theta >= 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 *
042 * @since 1.0
043 */
044public class AhrensDieterMarsagliaTsangGammaSampler
045    extends SamplerBase
046    implements ContinuousSampler {
047    /** 1/3 */
048    private static final double ONE_THIRD = 1d / 3;
049    /** The shape parameter. */
050    private final double theta;
051    /** The alpha parameter. */
052    private final double alpha;
053    /** Inverse of "theta". */
054    private final double oneOverTheta;
055    /** Optimization (see code). */
056    private final double bGSOptim;
057    /** Optimization (see code). */
058    private final double dOptim;
059    /** Optimization (see code). */
060    private final double cOptim;
061    /** Gaussian sampling. */
062    private final NormalizedGaussianSampler gaussian;
063    /** Underlying source of randomness. */
064    private final UniformRandomProvider rng;
065
066    /**
067     * @param rng Generator of uniformly distributed random numbers.
068     * @param alpha Alpha parameter of the distribution.
069     * @param theta Theta parameter of the distribution.
070     */
071    public AhrensDieterMarsagliaTsangGammaSampler(UniformRandomProvider rng,
072                                                  double alpha,
073                                                  double theta) {
074        super(null);
075        this.rng = rng;
076        this.alpha = alpha;
077        this.theta = theta;
078        gaussian = new ZigguratNormalizedGaussianSampler(rng);
079        oneOverTheta = 1 / theta;
080        bGSOptim = 1 + theta / Math.E;
081        dOptim = theta - ONE_THIRD;
082        cOptim = ONE_THIRD / Math.sqrt(dOptim);
083    }
084
085    /** {@inheritDoc} */
086    @Override
087    public double sample() {
088        if (theta < 1) {
089            // [1]: p. 228, Algorithm GS.
090
091            while (true) {
092                // Step 1:
093                final double u = rng.nextDouble();
094                final double p = bGSOptim * u;
095
096                if (p <= 1) {
097                    // Step 2:
098
099                    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}