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 = ZigguratSampler.NormalizedGaussian.of(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 = ZigguratSampler.NormalizedGaussian.of(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}