| /* |
| * 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) { |
| if (alpha <= 0) { |
| throw new IllegalArgumentException("alpha is not strictly positive: " + alpha); |
| } |
| if (theta <= 0) { |
| throw new IllegalArgumentException("theta is not strictly positive: " + theta); |
| } |
| 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 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 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 = new ZigguratNormalizedGaussianSampler(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 = new ZigguratNormalizedGaussianSampler(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); |
| } |
| } |