GammaDistribution.java

  1. /*
  2.  * Licensed to the Apache Software Foundation (ASF) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * The ASF licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *      http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */
  17. package org.apache.commons.math3.distribution;

  18. import org.apache.commons.math3.exception.NotStrictlyPositiveException;
  19. import org.apache.commons.math3.exception.util.LocalizedFormats;
  20. import org.apache.commons.math3.random.RandomGenerator;
  21. import org.apache.commons.math3.random.Well19937c;
  22. import org.apache.commons.math3.special.Gamma;
  23. import org.apache.commons.math3.util.FastMath;

  24. /**
  25.  * Implementation of the Gamma distribution.
  26.  *
  27.  * @see <a href="http://en.wikipedia.org/wiki/Gamma_distribution">Gamma distribution (Wikipedia)</a>
  28.  * @see <a href="http://mathworld.wolfram.com/GammaDistribution.html">Gamma distribution (MathWorld)</a>
  29.  */
  30. public class GammaDistribution extends AbstractRealDistribution {
  31.     /**
  32.      * Default inverse cumulative probability accuracy.
  33.      * @since 2.1
  34.      */
  35.     public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9;
  36.     /** Serializable version identifier. */
  37.     private static final long serialVersionUID = 20120524L;
  38.     /** The shape parameter. */
  39.     private final double shape;
  40.     /** The scale parameter. */
  41.     private final double scale;
  42.     /**
  43.      * The constant value of {@code shape + g + 0.5}, where {@code g} is the
  44.      * Lanczos constant {@link Gamma#LANCZOS_G}.
  45.      */
  46.     private final double shiftedShape;
  47.     /**
  48.      * The constant value of
  49.      * {@code shape / scale * sqrt(e / (2 * pi * (shape + g + 0.5))) / L(shape)},
  50.      * where {@code L(shape)} is the Lanczos approximation returned by
  51.      * {@link Gamma#lanczos(double)}. This prefactor is used in
  52.      * {@link #density(double)}, when no overflow occurs with the natural
  53.      * calculation.
  54.      */
  55.     private final double densityPrefactor1;
  56.     /**
  57.      * The constant value of
  58.      * {@code log(shape / scale * sqrt(e / (2 * pi * (shape + g + 0.5))) / L(shape))},
  59.      * where {@code L(shape)} is the Lanczos approximation returned by
  60.      * {@link Gamma#lanczos(double)}. This prefactor is used in
  61.      * {@link #logDensity(double)}, when no overflow occurs with the natural
  62.      * calculation.
  63.      */
  64.     private final double logDensityPrefactor1;
  65.     /**
  66.      * The constant value of
  67.      * {@code shape * sqrt(e / (2 * pi * (shape + g + 0.5))) / L(shape)},
  68.      * where {@code L(shape)} is the Lanczos approximation returned by
  69.      * {@link Gamma#lanczos(double)}. This prefactor is used in
  70.      * {@link #density(double)}, when overflow occurs with the natural
  71.      * calculation.
  72.      */
  73.     private final double densityPrefactor2;
  74.     /**
  75.      * The constant value of
  76.      * {@code log(shape * sqrt(e / (2 * pi * (shape + g + 0.5))) / L(shape))},
  77.      * where {@code L(shape)} is the Lanczos approximation returned by
  78.      * {@link Gamma#lanczos(double)}. This prefactor is used in
  79.      * {@link #logDensity(double)}, when overflow occurs with the natural
  80.      * calculation.
  81.      */
  82.     private final double logDensityPrefactor2;
  83.     /**
  84.      * Lower bound on {@code y = x / scale} for the selection of the computation
  85.      * method in {@link #density(double)}. For {@code y <= minY}, the natural
  86.      * calculation overflows.
  87.      */
  88.     private final double minY;
  89.     /**
  90.      * Upper bound on {@code log(y)} ({@code y = x / scale}) for the selection
  91.      * of the computation method in {@link #density(double)}. For
  92.      * {@code log(y) >= maxLogY}, the natural calculation overflows.
  93.      */
  94.     private final double maxLogY;
  95.     /** Inverse cumulative probability accuracy. */
  96.     private final double solverAbsoluteAccuracy;

  97.     /**
  98.      * Creates a new gamma distribution with specified values of the shape and
  99.      * scale parameters.
  100.      * <p>
  101.      * <b>Note:</b> this constructor will implicitly create an instance of
  102.      * {@link Well19937c} as random generator to be used for sampling only (see
  103.      * {@link #sample()} and {@link #sample(int)}). In case no sampling is
  104.      * needed for the created distribution, it is advised to pass {@code null}
  105.      * as random generator via the appropriate constructors to avoid the
  106.      * additional initialisation overhead.
  107.      *
  108.      * @param shape the shape parameter
  109.      * @param scale the scale parameter
  110.      * @throws NotStrictlyPositiveException if {@code shape <= 0} or
  111.      * {@code scale <= 0}.
  112.      */
  113.     public GammaDistribution(double shape, double scale) throws NotStrictlyPositiveException {
  114.         this(shape, scale, DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
  115.     }

  116.     /**
  117.      * Creates a new gamma distribution with specified values of the shape and
  118.      * scale parameters.
  119.      * <p>
  120.      * <b>Note:</b> this constructor will implicitly create an instance of
  121.      * {@link Well19937c} as random generator to be used for sampling only (see
  122.      * {@link #sample()} and {@link #sample(int)}). In case no sampling is
  123.      * needed for the created distribution, it is advised to pass {@code null}
  124.      * as random generator via the appropriate constructors to avoid the
  125.      * additional initialisation overhead.
  126.      *
  127.      * @param shape the shape parameter
  128.      * @param scale the scale parameter
  129.      * @param inverseCumAccuracy the maximum absolute error in inverse
  130.      * cumulative probability estimates (defaults to
  131.      * {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}).
  132.      * @throws NotStrictlyPositiveException if {@code shape <= 0} or
  133.      * {@code scale <= 0}.
  134.      * @since 2.1
  135.      */
  136.     public GammaDistribution(double shape, double scale, double inverseCumAccuracy)
  137.         throws NotStrictlyPositiveException {
  138.         this(new Well19937c(), shape, scale, inverseCumAccuracy);
  139.     }

  140.     /**
  141.      * Creates a Gamma distribution.
  142.      *
  143.      * @param rng Random number generator.
  144.      * @param shape the shape parameter
  145.      * @param scale the scale parameter
  146.      * @throws NotStrictlyPositiveException if {@code shape <= 0} or
  147.      * {@code scale <= 0}.
  148.      * @since 3.3
  149.      */
  150.     public GammaDistribution(RandomGenerator rng, double shape, double scale)
  151.         throws NotStrictlyPositiveException {
  152.         this(rng, shape, scale, DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
  153.     }

  154.     /**
  155.      * Creates a Gamma distribution.
  156.      *
  157.      * @param rng Random number generator.
  158.      * @param shape the shape parameter
  159.      * @param scale the scale parameter
  160.      * @param inverseCumAccuracy the maximum absolute error in inverse
  161.      * cumulative probability estimates (defaults to
  162.      * {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}).
  163.      * @throws NotStrictlyPositiveException if {@code shape <= 0} or
  164.      * {@code scale <= 0}.
  165.      * @since 3.1
  166.      */
  167.     public GammaDistribution(RandomGenerator rng,
  168.                              double shape,
  169.                              double scale,
  170.                              double inverseCumAccuracy)
  171.         throws NotStrictlyPositiveException {
  172.         super(rng);

  173.         if (shape <= 0) {
  174.             throw new NotStrictlyPositiveException(LocalizedFormats.SHAPE, shape);
  175.         }
  176.         if (scale <= 0) {
  177.             throw new NotStrictlyPositiveException(LocalizedFormats.SCALE, scale);
  178.         }

  179.         this.shape = shape;
  180.         this.scale = scale;
  181.         this.solverAbsoluteAccuracy = inverseCumAccuracy;
  182.         this.shiftedShape = shape + Gamma.LANCZOS_G + 0.5;
  183.         final double aux = FastMath.E / (2.0 * FastMath.PI * shiftedShape);
  184.         this.densityPrefactor2 = shape * FastMath.sqrt(aux) / Gamma.lanczos(shape);
  185.         this.logDensityPrefactor2 = FastMath.log(shape) + 0.5 * FastMath.log(aux) -
  186.                                     FastMath.log(Gamma.lanczos(shape));
  187.         this.densityPrefactor1 = this.densityPrefactor2 / scale *
  188.                 FastMath.pow(shiftedShape, -shape) *
  189.                 FastMath.exp(shape + Gamma.LANCZOS_G);
  190.         this.logDensityPrefactor1 = this.logDensityPrefactor2 - FastMath.log(scale) -
  191.                 FastMath.log(shiftedShape) * shape +
  192.                 shape + Gamma.LANCZOS_G;
  193.         this.minY = shape + Gamma.LANCZOS_G - FastMath.log(Double.MAX_VALUE);
  194.         this.maxLogY = FastMath.log(Double.MAX_VALUE) / (shape - 1.0);
  195.     }

  196.     /**
  197.      * Returns the shape parameter of {@code this} distribution.
  198.      *
  199.      * @return the shape parameter
  200.      * @deprecated as of version 3.1, {@link #getShape()} should be preferred.
  201.      * This method will be removed in version 4.0.
  202.      */
  203.     @Deprecated
  204.     public double getAlpha() {
  205.         return shape;
  206.     }

  207.     /**
  208.      * Returns the shape parameter of {@code this} distribution.
  209.      *
  210.      * @return the shape parameter
  211.      * @since 3.1
  212.      */
  213.     public double getShape() {
  214.         return shape;
  215.     }

  216.     /**
  217.      * Returns the scale parameter of {@code this} distribution.
  218.      *
  219.      * @return the scale parameter
  220.      * @deprecated as of version 3.1, {@link #getScale()} should be preferred.
  221.      * This method will be removed in version 4.0.
  222.      */
  223.     @Deprecated
  224.     public double getBeta() {
  225.         return scale;
  226.     }

  227.     /**
  228.      * Returns the scale parameter of {@code this} distribution.
  229.      *
  230.      * @return the scale parameter
  231.      * @since 3.1
  232.      */
  233.     public double getScale() {
  234.         return scale;
  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.     /** {@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.      * {@inheritDoc}
  322.      *
  323.      * The implementation of this method is based on:
  324.      * <ul>
  325.      *  <li>
  326.      *   <a href="http://mathworld.wolfram.com/Chi-SquaredDistribution.html">
  327.      *    Chi-Squared Distribution</a>, equation (9).
  328.      *  </li>
  329.      *  <li>Casella, G., & Berger, R. (1990). <i>Statistical Inference</i>.
  330.      *    Belmont, CA: Duxbury Press.
  331.      *  </li>
  332.      * </ul>
  333.      */
  334.     public double cumulativeProbability(double x) {
  335.         double ret;

  336.         if (x <= 0) {
  337.             ret = 0;
  338.         } else {
  339.             ret = Gamma.regularizedGammaP(shape, x / scale);
  340.         }

  341.         return ret;
  342.     }

  343.     /** {@inheritDoc} */
  344.     @Override
  345.     protected double getSolverAbsoluteAccuracy() {
  346.         return solverAbsoluteAccuracy;
  347.     }

  348.     /**
  349.      * {@inheritDoc}
  350.      *
  351.      * For shape parameter {@code alpha} and scale parameter {@code beta}, the
  352.      * mean is {@code alpha * beta}.
  353.      */
  354.     public double getNumericalMean() {
  355.         return shape * scale;
  356.     }

  357.     /**
  358.      * {@inheritDoc}
  359.      *
  360.      * For shape parameter {@code alpha} and scale parameter {@code beta}, the
  361.      * variance is {@code alpha * beta^2}.
  362.      *
  363.      * @return {@inheritDoc}
  364.      */
  365.     public double getNumericalVariance() {
  366.         return shape * scale * scale;
  367.     }

  368.     /**
  369.      * {@inheritDoc}
  370.      *
  371.      * The lower bound of the support is always 0 no matter the parameters.
  372.      *
  373.      * @return lower bound of the support (always 0)
  374.      */
  375.     public double getSupportLowerBound() {
  376.         return 0;
  377.     }

  378.     /**
  379.      * {@inheritDoc}
  380.      *
  381.      * The upper bound of the support is always positive infinity
  382.      * no matter the parameters.
  383.      *
  384.      * @return upper bound of the support (always Double.POSITIVE_INFINITY)
  385.      */
  386.     public double getSupportUpperBound() {
  387.         return Double.POSITIVE_INFINITY;
  388.     }

  389.     /** {@inheritDoc} */
  390.     public boolean isSupportLowerBoundInclusive() {
  391.         return true;
  392.     }

  393.     /** {@inheritDoc} */
  394.     public boolean isSupportUpperBoundInclusive() {
  395.         return false;
  396.     }

  397.     /**
  398.      * {@inheritDoc}
  399.      *
  400.      * The support of this distribution is connected.
  401.      *
  402.      * @return {@code true}
  403.      */
  404.     public boolean isSupportConnected() {
  405.         return true;
  406.     }

  407.     /**
  408.      * <p>This implementation uses the following algorithms: </p>
  409.      *
  410.      * <p>For 0 < shape < 1: <br/>
  411.      * Ahrens, J. H. and Dieter, U., <i>Computer methods for
  412.      * sampling from gamma, beta, Poisson and binomial distributions.</i>
  413.      * Computing, 12, 223-246, 1974.</p>
  414.      *
  415.      * <p>For shape >= 1: <br/>
  416.      * Marsaglia and Tsang, <i>A Simple Method for Generating
  417.      * Gamma Variables.</i> ACM Transactions on Mathematical Software,
  418.      * Volume 26 Issue 3, September, 2000.</p>
  419.      *
  420.      * @return random value sampled from the Gamma(shape, scale) distribution
  421.      */
  422.     @Override
  423.     public double sample()  {
  424.         if (shape < 1) {
  425.             // [1]: p. 228, Algorithm GS

  426.             while (true) {
  427.                 // Step 1:
  428.                 final double u = random.nextDouble();
  429.                 final double bGS = 1 + shape / FastMath.E;
  430.                 final double p = bGS * u;

  431.                 if (p <= 1) {
  432.                     // Step 2:

  433.                     final double x = FastMath.pow(p, 1 / shape);
  434.                     final double u2 = random.nextDouble();

  435.                     if (u2 > FastMath.exp(-x)) {
  436.                         // Reject
  437.                         continue;
  438.                     } else {
  439.                         return scale * x;
  440.                     }
  441.                 } else {
  442.                     // Step 3:

  443.                     final double x = -1 * FastMath.log((bGS - p) / shape);
  444.                     final double u2 = random.nextDouble();

  445.                     if (u2 > FastMath.pow(x, shape - 1)) {
  446.                         // Reject
  447.                         continue;
  448.                     } else {
  449.                         return scale * x;
  450.                     }
  451.                 }
  452.             }
  453.         }

  454.         // Now shape >= 1

  455.         final double d = shape - 0.333333333333333333;
  456.         final double c = 1 / (3 * FastMath.sqrt(d));

  457.         while (true) {
  458.             final double x = random.nextGaussian();
  459.             final double v = (1 + c * x) * (1 + c * x) * (1 + c * x);

  460.             if (v <= 0) {
  461.                 continue;
  462.             }

  463.             final double x2 = x * x;
  464.             final double u = random.nextDouble();

  465.             // Squeeze
  466.             if (u < 1 - 0.0331 * x2 * x2) {
  467.                 return scale * d * v;
  468.             }

  469.             if (FastMath.log(u) < 0.5 * x2 + d * (1 - v + FastMath.log(v))) {
  470.                 return scale * d * v;
  471.             }
  472.         }
  473.     }
  474. }