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 < alpha < 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 alpha >= 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 * <p>Sampling uses:</p>
043 *
044 * <ul>
045 *   <li>{@link UniformRandomProvider#nextDouble()} (both algorithms)
046 *   <li>{@link UniformRandomProvider#nextLong()} (only for {@code alpha >= 1})
047 * </ul>
048 *
049 * @see <a href="http://mathworld.wolfram.com/GammaDistribution.html">MathWorld Gamma distribution</a>
050 * @see <a href="https://en.wikipedia.org/wiki/Gamma_distribution">Wikipedia Gamma distribution</a>
051 * @since 1.0
052 */
053public class AhrensDieterMarsagliaTsangGammaSampler
054    extends SamplerBase
055    implements SharedStateContinuousSampler {
056    /** The appropriate gamma sampler for the parameters. */
057    private final SharedStateContinuousSampler delegate;
058
059    /**
060     * Base class for a sampler from the Gamma distribution.
061     */
062    private abstract static class BaseGammaSampler
063        implements SharedStateContinuousSampler {
064
065        /** Underlying source of randomness. */
066        protected final UniformRandomProvider rng;
067        /** The alpha parameter. This is a shape parameter. */
068        protected final double alpha;
069        /** The theta parameter. This is a scale parameter. */
070        protected final double theta;
071
072        /**
073         * @param rng Generator of uniformly distributed random numbers.
074         * @param alpha Alpha parameter of the distribution.
075         * @param theta Theta parameter of the distribution.
076         * @throws IllegalArgumentException if {@code alpha <= 0} or {@code theta <= 0}
077         */
078        BaseGammaSampler(UniformRandomProvider rng,
079                         double alpha,
080                         double theta) {
081            // Validation before java.lang.Object constructor exits prevents partially initialized object
082            this(InternalUtils.requireStrictlyPositive(alpha, "alpha"),
083                 InternalUtils.requireStrictlyPositive(theta, "theta"),
084                 rng);
085        }
086
087        /**
088         * @param alpha Alpha parameter of the distribution.
089         * @param theta Theta parameter of the distribution.
090         * @param rng Generator of uniformly distributed random numbers.
091         */
092        private BaseGammaSampler(double alpha,
093                                 double theta,
094                                 UniformRandomProvider rng) {
095            this.rng = rng;
096            this.alpha = alpha;
097            this.theta = theta;
098        }
099
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}