AhrensDieterMarsagliaTsangGammaSampler.java

/*
 * Licensed to the Apache Software Foundation (ASF) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * The ASF licenses this file to You under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance with
 * the License.  You may obtain a copy of the License at
 *
 *      http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
package org.apache.commons.rng.sampling.distribution;

import org.apache.commons.rng.UniformRandomProvider;

/**
 * Sampling from the <a href="http://mathworld.wolfram.com/GammaDistribution.html">gamma distribution</a>.
 * <ul>
 *  <li>
 *   For {@code 0 < alpha < 1}:
 *   <blockquote>
 *    Ahrens, J. H. and Dieter, U.,
 *    <i>Computer methods for sampling from gamma, beta, Poisson and binomial distributions,</i>
 *    Computing, 12, 223-246, 1974.
 *   </blockquote>
 *  </li>
 *  <li>
 *  For {@code alpha >= 1}:
 *   <blockquote>
 *   Marsaglia and Tsang, <i>A Simple Method for Generating
 *   Gamma Variables.</i> ACM Transactions on Mathematical Software,
 *   Volume 26 Issue 3, September, 2000.
 *   </blockquote>
 *  </li>
 * </ul>
 *
 * <p>Sampling uses:</p>
 *
 * <ul>
 *   <li>{@link UniformRandomProvider#nextDouble()} (both algorithms)
 *   <li>{@link UniformRandomProvider#nextLong()} (only for {@code alpha >= 1})
 * </ul>
 *
 * @see <a href="http://mathworld.wolfram.com/GammaDistribution.html">MathWorld Gamma distribution</a>
 * @see <a href="https://en.wikipedia.org/wiki/Gamma_distribution">Wikipedia Gamma distribution</a>
 * @since 1.0
 */
public class AhrensDieterMarsagliaTsangGammaSampler
    extends SamplerBase
    implements SharedStateContinuousSampler {
    /** The appropriate gamma sampler for the parameters. */
    private final SharedStateContinuousSampler delegate;

    /**
     * Base class for a sampler from the Gamma distribution.
     */
    private abstract static class BaseGammaSampler
        implements SharedStateContinuousSampler {

        /** Underlying source of randomness. */
        protected final UniformRandomProvider rng;
        /** The alpha parameter. This is a shape parameter. */
        protected final double alpha;
        /** The theta parameter. This is a scale parameter. */
        protected final double theta;

        /**
         * @param rng Generator of uniformly distributed random numbers.
         * @param alpha Alpha parameter of the distribution.
         * @param theta Theta parameter of the distribution.
         * @throws IllegalArgumentException if {@code alpha <= 0} or {@code theta <= 0}
         */
        BaseGammaSampler(UniformRandomProvider rng,
                         double alpha,
                         double theta) {
            // Validation before java.lang.Object constructor exits prevents partially initialized object
            this(InternalUtils.requireStrictlyPositive(alpha, "alpha"),
                 InternalUtils.requireStrictlyPositive(theta, "theta"),
                 rng);
        }

        /**
         * @param alpha Alpha parameter of the distribution.
         * @param theta Theta parameter of the distribution.
         * @param rng Generator of uniformly distributed random numbers.
         */
        private BaseGammaSampler(double alpha,
                                 double theta,
                                 UniformRandomProvider rng) {
            this.rng = rng;
            this.alpha = alpha;
            this.theta = theta;
        }

        /**
         * @param rng Generator of uniformly distributed random numbers.
         * @param source Source to copy.
         */
        BaseGammaSampler(UniformRandomProvider rng,
                         BaseGammaSampler source) {
            this.rng = rng;
            this.alpha = source.alpha;
            this.theta = source.theta;
        }

        /** {@inheritDoc} */
        @Override
        public String toString() {
            return "Ahrens-Dieter-Marsaglia-Tsang Gamma deviate [" + rng.toString() + "]";
        }
    }

    /**
     * Class to sample from the Gamma distribution when {@code 0 < alpha < 1}.
     *
     * <blockquote>
     *  Ahrens, J. H. and Dieter, U.,
     *  <i>Computer methods for sampling from gamma, beta, Poisson and binomial distributions,</i>
     *  Computing, 12, 223-246, 1974.
     * </blockquote>
     */
    private static final class AhrensDieterGammaSampler
        extends BaseGammaSampler {

        /** Inverse of "alpha". */
        private final double oneOverAlpha;
        /** Optimization (see code). */
        private final double bGSOptim;

        /**
         * @param rng Generator of uniformly distributed random numbers.
         * @param alpha Alpha parameter of the distribution.
         * @param theta Theta parameter of the distribution.
         * @throws IllegalArgumentException if {@code alpha <= 0} or {@code theta <= 0}
         */
        AhrensDieterGammaSampler(UniformRandomProvider rng,
                                 double alpha,
                                 double theta) {
            super(rng, alpha, theta);
            oneOverAlpha = 1 / alpha;
            bGSOptim = 1 + alpha / Math.E;
        }

        /**
         * @param rng Generator of uniformly distributed random numbers.
         * @param source Source to copy.
         */
        AhrensDieterGammaSampler(UniformRandomProvider rng,
                                 AhrensDieterGammaSampler source) {
            super(rng, source);
            oneOverAlpha = source.oneOverAlpha;
            bGSOptim = source.bGSOptim;
        }

        @Override
        public double sample() {
            // [1]: p. 228, Algorithm GS.

            while (true) {
                // Step 1:
                final double u = rng.nextDouble();
                final double p = bGSOptim * u;

                if (p <= 1) {
                    // Step 2:
                    final double x = Math.pow(p, oneOverAlpha);
                    final double u2 = rng.nextDouble();

                    if (u2 > Math.exp(-x)) {
                        // Reject.
                        continue;
                    }
                    return theta * x;
                }

                // Step 3:
                final double x = -Math.log((bGSOptim - p) * oneOverAlpha);
                final double u2 = rng.nextDouble();

                if (u2 <= Math.pow(x, alpha - 1)) {
                    return theta * x;
                }
                // Reject and continue.
            }
        }

        @Override
        public SharedStateContinuousSampler withUniformRandomProvider(UniformRandomProvider rng) {
            return new AhrensDieterGammaSampler(rng, this);
        }
    }

    /**
     * Class to sample from the Gamma distribution when the {@code alpha >= 1}.
     *
     * <blockquote>
     *  Marsaglia and Tsang, <i>A Simple Method for Generating
     *  Gamma Variables.</i> ACM Transactions on Mathematical Software,
     *  Volume 26 Issue 3, September, 2000.
     * </blockquote>
     */
    private static final class MarsagliaTsangGammaSampler
        extends BaseGammaSampler {

        /** 1/3. */
        private static final double ONE_THIRD = 1d / 3;

        /** Optimization (see code). */
        private final double dOptim;
        /** Optimization (see code). */
        private final double cOptim;
        /** Gaussian sampling. */
        private final NormalizedGaussianSampler gaussian;

        /**
         * @param rng Generator of uniformly distributed random numbers.
         * @param alpha Alpha parameter of the distribution.
         * @param theta Theta parameter of the distribution.
         * @throws IllegalArgumentException if {@code alpha <= 0} or {@code theta <= 0}
         */
        MarsagliaTsangGammaSampler(UniformRandomProvider rng,
                                   double alpha,
                                   double theta) {
            super(rng, alpha, theta);
            gaussian = ZigguratSampler.NormalizedGaussian.of(rng);
            dOptim = alpha - ONE_THIRD;
            cOptim = ONE_THIRD / Math.sqrt(dOptim);
        }

        /**
         * @param rng Generator of uniformly distributed random numbers.
         * @param source Source to copy.
         */
        MarsagliaTsangGammaSampler(UniformRandomProvider rng,
                                   MarsagliaTsangGammaSampler source) {
            super(rng, source);
            gaussian = ZigguratSampler.NormalizedGaussian.of(rng);
            dOptim = source.dOptim;
            cOptim = source.cOptim;
        }

        @Override
        public double sample() {
            while (true) {
                final double x = gaussian.sample();
                final double oPcTx = 1 + cOptim * x;
                final double v = oPcTx * oPcTx * oPcTx;

                if (v <= 0) {
                    continue;
                }

                final double x2 = x * x;
                final double u = rng.nextDouble();

                // Squeeze.
                if (u < 1 - 0.0331 * x2 * x2) {
                    return theta * dOptim * v;
                }

                if (Math.log(u) < 0.5 * x2 + dOptim * (1 - v + Math.log(v))) {
                    return theta * dOptim * v;
                }
            }
        }

        @Override
        public SharedStateContinuousSampler withUniformRandomProvider(UniformRandomProvider rng) {
            return new MarsagliaTsangGammaSampler(rng, this);
        }
    }

    /**
     * This instance delegates sampling. Use the factory method
     * {@link #of(UniformRandomProvider, double, double)} to create an optimal sampler.
     *
     * @param rng Generator of uniformly distributed random numbers.
     * @param alpha Alpha parameter of the distribution (this is a shape parameter).
     * @param theta Theta parameter of the distribution (this is a scale parameter).
     * @throws IllegalArgumentException if {@code alpha <= 0} or {@code theta <= 0}
     */
    public AhrensDieterMarsagliaTsangGammaSampler(UniformRandomProvider rng,
                                                  double alpha,
                                                  double theta) {
        super(null);
        delegate = of(rng, alpha, theta);
    }

    /** {@inheritDoc} */
    @Override
    public double sample() {
        return delegate.sample();
    }

    /** {@inheritDoc} */
    @Override
    public String toString() {
        return delegate.toString();
    }

    /**
     * {@inheritDoc}
     *
     * @since 1.3
     */
    @Override
    public SharedStateContinuousSampler withUniformRandomProvider(UniformRandomProvider rng) {
        // Direct return of the optimised sampler
        return delegate.withUniformRandomProvider(rng);
    }

    /**
     * Creates a new gamma distribution sampler.
     *
     * @param rng Generator of uniformly distributed random numbers.
     * @param alpha Alpha parameter of the distribution (this is a shape parameter).
     * @param theta Theta parameter of the distribution (this is a scale parameter).
     * @return the sampler
     * @throws IllegalArgumentException if {@code alpha <= 0} or {@code theta <= 0}
     * @since 1.3
     */
    public static SharedStateContinuousSampler of(UniformRandomProvider rng,
                                                  double alpha,
                                                  double theta) {
        // Each sampler should check the input arguments.
        return alpha < 1 ?
                new AhrensDieterGammaSampler(rng, alpha, theta) :
                new MarsagliaTsangGammaSampler(rng, alpha, theta);
    }
}