    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    }