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