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 */
017 package org.apache.commons.math3.distribution;
018
019 import org.apache.commons.math3.exception.NotStrictlyPositiveException;
020 import org.apache.commons.math3.exception.util.LocalizedFormats;
021 import org.apache.commons.math3.special.Gamma;
022 import org.apache.commons.math3.util.FastMath;
023 import org.apache.commons.math3.random.RandomGenerator;
024 import org.apache.commons.math3.random.Well19937c;
025
026 /**
027 * Implementation of the Gamma distribution.
028 *
029 * @see <a href="http://en.wikipedia.org/wiki/Gamma_distribution">Gamma distribution (Wikipedia)</a>
030 * @see <a href="http://mathworld.wolfram.com/GammaDistribution.html">Gamma distribution (MathWorld)</a>
031 * @version $Id: GammaDistribution.java 1422195 2012-12-15 06:45:18Z psteitz $
032 */
033 public class GammaDistribution extends AbstractRealDistribution {
034 /**
035 * Default inverse cumulative probability accuracy.
036 * @since 2.1
037 */
038 public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9;
039 /** Serializable version identifier. */
040 private static final long serialVersionUID = 20120524L;
041 /** The shape parameter. */
042 private final double shape;
043 /** The scale parameter. */
044 private final double scale;
045 /**
046 * The constant value of {@code shape + g + 0.5}, where {@code g} is the
047 * Lanczos constant {@link Gamma#LANCZOS_G}.
048 */
049 private final double shiftedShape;
050 /**
051 * The constant value of
052 * {@code shape / scale * sqrt(e / (2 * pi * (shape + g + 0.5))) / L(shape)},
053 * where {@code L(shape)} is the Lanczos approximation returned by
054 * {@link Gamma#lanczos(double)}. This prefactor is used in
055 * {@link #density(double)}, when no overflow occurs with the natural
056 * calculation.
057 */
058 private final double densityPrefactor1;
059 /**
060 * The constant value of
061 * {@code shape * sqrt(e / (2 * pi * (shape + g + 0.5))) / L(shape)},
062 * where {@code L(shape)} is the Lanczos approximation returned by
063 * {@link Gamma#lanczos(double)}. This prefactor is used in
064 * {@link #density(double)}, when overflow occurs with the natural
065 * calculation.
066 */
067 private final double densityPrefactor2;
068 /**
069 * Lower bound on {@code y = x / scale} for the selection of the computation
070 * method in {@link #density(double)}. For {@code y <= minY}, the natural
071 * calculation overflows.
072 */
073 private final double minY;
074 /**
075 * Upper bound on {@code log(y)} ({@code y = x / scale}) for the selection
076 * of the computation method in {@link #density(double)}. For
077 * {@code log(y) >= maxLogY}, the natural calculation overflows.
078 */
079 private final double maxLogY;
080 /** Inverse cumulative probability accuracy. */
081 private final double solverAbsoluteAccuracy;
082
083 /**
084 * Creates a new gamma distribution with specified values of the shape and
085 * scale parameters.
086 *
087 * @param shape the shape parameter
088 * @param scale the scale parameter
089 * @throws NotStrictlyPositiveException if {@code shape <= 0} or
090 * {@code scale <= 0}.
091 */
092 public GammaDistribution(double shape, double scale) throws NotStrictlyPositiveException {
093 this(shape, scale, DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
094 }
095
096 /**
097 * Creates a new gamma distribution with specified values of the shape and
098 * scale parameters.
099 *
100 * @param shape the shape parameter
101 * @param scale the scale parameter
102 * @param inverseCumAccuracy the maximum absolute error in inverse
103 * cumulative probability estimates (defaults to
104 * {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}).
105 * @throws NotStrictlyPositiveException if {@code shape <= 0} or
106 * {@code scale <= 0}.
107 * @since 2.1
108 */
109 public GammaDistribution(double shape, double scale, double inverseCumAccuracy)
110 throws NotStrictlyPositiveException {
111 this(new Well19937c(), shape, scale, inverseCumAccuracy);
112 }
113
114 /**
115 * Creates a Gamma distribution.
116 *
117 * @param rng Random number generator.
118 * @param shape the shape parameter
119 * @param scale the scale parameter
120 * @param inverseCumAccuracy the maximum absolute error in inverse
121 * cumulative probability estimates (defaults to
122 * {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}).
123 * @throws NotStrictlyPositiveException if {@code shape <= 0} or
124 * {@code scale <= 0}.
125 * @since 3.1
126 */
127 public GammaDistribution(RandomGenerator rng,
128 double shape,
129 double scale,
130 double inverseCumAccuracy)
131 throws NotStrictlyPositiveException {
132 super(rng);
133
134 if (shape <= 0) {
135 throw new NotStrictlyPositiveException(LocalizedFormats.SHAPE, shape);
136 }
137 if (scale <= 0) {
138 throw new NotStrictlyPositiveException(LocalizedFormats.SCALE, scale);
139 }
140
141 this.shape = shape;
142 this.scale = scale;
143 this.solverAbsoluteAccuracy = inverseCumAccuracy;
144 this.shiftedShape = shape + Gamma.LANCZOS_G + 0.5;
145 final double aux = FastMath.E / (2.0 * FastMath.PI * shiftedShape);
146 this.densityPrefactor2 = shape * FastMath.sqrt(aux) / Gamma.lanczos(shape);
147 this.densityPrefactor1 = this.densityPrefactor2 / scale *
148 FastMath.pow(shiftedShape, -shape) *
149 FastMath.exp(shape + Gamma.LANCZOS_G);
150 this.minY = shape + Gamma.LANCZOS_G - FastMath.log(Double.MAX_VALUE);
151 this.maxLogY = FastMath.log(Double.MAX_VALUE) / (shape - 1.0);
152 }
153
154 /**
155 * Returns the shape parameter of {@code this} distribution.
156 *
157 * @return the shape parameter
158 * @deprecated as of version 3.1, {@link #getShape()} should be preferred.
159 * This method will be removed in version 4.0.
160 */
161 @Deprecated
162 public double getAlpha() {
163 return shape;
164 }
165
166 /**
167 * Returns the shape parameter of {@code this} distribution.
168 *
169 * @return the shape parameter
170 * @since 3.1
171 */
172 public double getShape() {
173 return shape;
174 }
175
176 /**
177 * Returns the scale parameter of {@code this} distribution.
178 *
179 * @return the scale parameter
180 * @deprecated as of version 3.1, {@link #getScale()} should be preferred.
181 * This method will be removed in version 4.0.
182 */
183 @Deprecated
184 public double getBeta() {
185 return scale;
186 }
187
188 /**
189 * Returns the scale parameter of {@code this} distribution.
190 *
191 * @return the scale parameter
192 * @since 3.1
193 */
194 public double getScale() {
195 return scale;
196 }
197
198 /** {@inheritDoc} */
199 public double density(double x) {
200 /* The present method must return the value of
201 *
202 * 1 x a - x
203 * ---------- (-) exp(---)
204 * x Gamma(a) b b
205 *
206 * where a is the shape parameter, and b the scale parameter.
207 * Substituting the Lanczos approximation of Gamma(a) leads to the
208 * following expression of the density
209 *
210 * a e 1 y a
211 * - sqrt(------------------) ---- (-----------) exp(a - y + g),
212 * x 2 pi (a + g + 0.5) L(a) a + g + 0.5
213 *
214 * where y = x / b. The above formula is the "natural" computation, which
215 * is implemented when no overflow is likely to occur. If overflow occurs
216 * with the natural computation, the following identity is used. It is
217 * based on the BOOST library
218 * http://www.boost.org/doc/libs/1_35_0/libs/math/doc/sf_and_dist/html/math_toolkit/special/sf_gamma/igamma.html
219 * Formula (15) needs adaptations, which are detailed below.
220 *
221 * y a
222 * (-----------) exp(a - y + g)
223 * a + g + 0.5
224 * y - a - g - 0.5 y (g + 0.5)
225 * = exp(a log1pm(---------------) - ----------- + g),
226 * a + g + 0.5 a + g + 0.5
227 *
228 * where log1pm(z) = log(1 + z) - z. Therefore, the value to be
229 * returned is
230 *
231 * a e 1
232 * - sqrt(------------------) ----
233 * x 2 pi (a + g + 0.5) L(a)
234 * y - a - g - 0.5 y (g + 0.5)
235 * * exp(a log1pm(---------------) - ----------- + g).
236 * a + g + 0.5 a + g + 0.5
237 */
238 if (x < 0) {
239 return 0;
240 }
241 final double y = x / scale;
242 if ((y <= minY) || (FastMath.log(y) >= maxLogY)) {
243 /*
244 * Overflow.
245 */
246 final double aux1 = (y - shiftedShape) / shiftedShape;
247 final double aux2 = shape * (FastMath.log1p(aux1) - aux1);
248 final double aux3 = -y * (Gamma.LANCZOS_G + 0.5) / shiftedShape +
249 Gamma.LANCZOS_G + aux2;
250 return densityPrefactor2 / x * FastMath.exp(aux3);
251 }
252 /*
253 * Natural calculation.
254 */
255 return densityPrefactor1 * FastMath.exp(-y) *
256 FastMath.pow(y, shape - 1);
257 }
258
259 /**
260 * {@inheritDoc}
261 *
262 * The implementation of this method is based on:
263 * <ul>
264 * <li>
265 * <a href="http://mathworld.wolfram.com/Chi-SquaredDistribution.html">
266 * Chi-Squared Distribution</a>, equation (9).
267 * </li>
268 * <li>Casella, G., & Berger, R. (1990). <i>Statistical Inference</i>.
269 * Belmont, CA: Duxbury Press.
270 * </li>
271 * </ul>
272 */
273 public double cumulativeProbability(double x) {
274 double ret;
275
276 if (x <= 0) {
277 ret = 0;
278 } else {
279 ret = Gamma.regularizedGammaP(shape, x / scale);
280 }
281
282 return ret;
283 }
284
285 /** {@inheritDoc} */
286 @Override
287 protected double getSolverAbsoluteAccuracy() {
288 return solverAbsoluteAccuracy;
289 }
290
291 /**
292 * {@inheritDoc}
293 *
294 * For shape parameter {@code alpha} and scale parameter {@code beta}, the
295 * mean is {@code alpha * beta}.
296 */
297 public double getNumericalMean() {
298 return shape * scale;
299 }
300
301 /**
302 * {@inheritDoc}
303 *
304 * For shape parameter {@code alpha} and scale parameter {@code beta}, the
305 * variance is {@code alpha * beta^2}.
306 *
307 * @return {@inheritDoc}
308 */
309 public double getNumericalVariance() {
310 return shape * scale * scale;
311 }
312
313 /**
314 * {@inheritDoc}
315 *
316 * The lower bound of the support is always 0 no matter the parameters.
317 *
318 * @return lower bound of the support (always 0)
319 */
320 public double getSupportLowerBound() {
321 return 0;
322 }
323
324 /**
325 * {@inheritDoc}
326 *
327 * The upper bound of the support is always positive infinity
328 * no matter the parameters.
329 *
330 * @return upper bound of the support (always Double.POSITIVE_INFINITY)
331 */
332 public double getSupportUpperBound() {
333 return Double.POSITIVE_INFINITY;
334 }
335
336 /** {@inheritDoc} */
337 public boolean isSupportLowerBoundInclusive() {
338 return true;
339 }
340
341 /** {@inheritDoc} */
342 public boolean isSupportUpperBoundInclusive() {
343 return false;
344 }
345
346 /**
347 * {@inheritDoc}
348 *
349 * The support of this distribution is connected.
350 *
351 * @return {@code true}
352 */
353 public boolean isSupportConnected() {
354 return true;
355 }
356
357 /**
358 * <p>This implementation uses the following algorithms: </p>
359 *
360 * <p>For 0 < shape < 1: <br/>
361 * Ahrens, J. H. and Dieter, U., <i>Computer methods for
362 * sampling from gamma, beta, Poisson and binomial distributions.</i>
363 * Computing, 12, 223-246, 1974.</p>
364 *
365 * <p>For shape >= 1: <br/>
366 * Marsaglia and Tsang, <i>A Simple Method for Generating
367 * Gamma Variables.</i> ACM Transactions on Mathematical Software,
368 * Volume 26 Issue 3, September, 2000.</p>
369 *
370 * @return random value sampled from the Gamma(shape, scale) distribution
371 */
372 @Override
373 public double sample() {
374 if (shape < 1) {
375 // [1]: p. 228, Algorithm GS
376
377 while (true) {
378 // Step 1:
379 final double u = random.nextDouble();
380 final double bGS = 1 + shape / FastMath.E;
381 final double p = bGS * u;
382
383 if (p <= 1) {
384 // Step 2:
385
386 final double x = FastMath.pow(p, 1 / shape);
387 final double u2 = random.nextDouble();
388
389 if (u2 > FastMath.exp(-x)) {
390 // Reject
391 continue;
392 } else {
393 return scale * x;
394 }
395 } else {
396 // Step 3:
397
398 final double x = -1 * FastMath.log((bGS - p) / shape);
399 final double u2 = random.nextDouble();
400
401 if (u2 > FastMath.pow(x, shape - 1)) {
402 // Reject
403 continue;
404 } else {
405 return scale * x;
406 }
407 }
408 }
409 }
410
411 // Now shape >= 1
412
413 final double d = shape - 0.333333333333333333;
414 final double c = 1 / (3 * FastMath.sqrt(d));
415
416 while (true) {
417 final double x = random.nextGaussian();
418 final double v = (1 + c * x) * (1 + c * x) * (1 + c * x);
419
420 if (v <= 0) {
421 continue;
422 }
423
424 final double x2 = x * x;
425 final double u = random.nextDouble();
426
427 // Squeeze
428 if (u < 1 - 0.0331 * x2 * x2) {
429 return scale * d * v;
430 }
431
432 if (FastMath.log(u) < 0.5 * x2 + d * (1 - v + FastMath.log(v))) {
433 return scale * d * v;
434 }
435 }
436 }
437 }