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}