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 < shape < 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 shape >= 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 */ 042public class AhrensDieterMarsagliaTsangGammaSampler 043 extends SamplerBase 044 implements ContinuousSampler { 045 /** The shape parameter. */ 046 private final double theta; 047 /** The alpha parameter. */ 048 private final double alpha; 049 /** Gaussian sampling. */ 050 private final BoxMullerGaussianSampler gaussian; 051 052 /** 053 * @param rng Generator of uniformly distributed random numbers. 054 * @param alpha Alpha parameter of the distribution. 055 * @param theta Theta parameter of the distribution. 056 */ 057 public AhrensDieterMarsagliaTsangGammaSampler(UniformRandomProvider rng, 058 double alpha, 059 double theta) { 060 super(rng); 061 this.alpha = alpha; 062 this.theta = theta; 063 gaussian = new BoxMullerGaussianSampler(rng, 0, 1); 064 } 065 066 /** {@inheritDoc} */ 067 @Override 068 public double sample() { 069 if (theta < 1) { 070 // [1]: p. 228, Algorithm GS. 071 072 while (true) { 073 // Step 1: 074 final double u = nextDouble(); 075 final double bGS = 1 + theta / Math.E; 076 final double p = bGS * u; 077 078 if (p <= 1) { 079 // Step 2: 080 081 final double x = Math.pow(p, 1 / theta); 082 final double u2 = nextDouble(); 083 084 if (u2 > Math.exp(-x)) { 085 // Reject. 086 continue; 087 } else { 088 return alpha * x; 089 } 090 } else { 091 // Step 3: 092 093 final double x = -1 * Math.log((bGS - p) / theta); 094 final double u2 = nextDouble(); 095 096 if (u2 > Math.pow(x, theta - 1)) { 097 // Reject. 098 continue; 099 } else { 100 return alpha * x; 101 } 102 } 103 } 104 } 105 106 // Now theta >= 1. 107 108 final double d = theta - 0.333333333333333333; 109 final double c = 1 / (3 * Math.sqrt(d)); 110 111 while (true) { 112 final double x = gaussian.sample(); 113 final double v = (1 + c * x) * (1 + c * x) * (1 + c * x); 114 115 if (v <= 0) { 116 continue; 117 } 118 119 final double x2 = x * x; 120 final double u = nextDouble(); 121 122 // Squeeze. 123 if (u < 1 - 0.0331 * x2 * x2) { 124 return alpha * d * v; 125 } 126 127 if (Math.log(u) < 0.5 * x2 + d * (1 - v + Math.log(v))) { 128 return alpha * d * v; 129 } 130 } 131 } 132 133 /** {@inheritDoc} */ 134 @Override 135 public String toString() { 136 return "Ahrens-Dieter-Marsaglia-Tsang Gamma deviate [" + super.toString() + "]"; 137 } 138}