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