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            if (alpha <= 0) {
082                throw new IllegalArgumentException("alpha is not strictly positive: " + alpha);
083            }
084            if (theta <= 0) {
085                throw new IllegalArgumentException("theta is not strictly positive: " + theta);
086            }
087            this.rng = rng;
088            this.alpha = alpha;
089            this.theta = theta;
090        }
091
092        /**
093         * @param rng Generator of uniformly distributed random numbers.
094         * @param source Source to copy.
095         */
096        BaseGammaSampler(UniformRandomProvider rng,
097                         BaseGammaSampler source) {
098            this.rng = rng;
099            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 = new ZigguratNormalizedGaussianSampler(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 = new ZigguratNormalizedGaussianSampler(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}