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.statistics.distribution;
018
019import org.apache.commons.numbers.gamma.LogGamma;
020import org.apache.commons.numbers.gamma.RegularizedGamma;
021import org.apache.commons.rng.UniformRandomProvider;
022import org.apache.commons.rng.sampling.distribution.AhrensDieterMarsagliaTsangGammaSampler;
023
024/**
025 * Implementation of the gamma distribution.
026 *
027 * <p>The probability density function of \( X \) is:
028 *
029 * <p>\[ f(x;k,\theta) = \frac{x^{k-1}e^{-x/\theta}}{\theta^k\Gamma(k)} \]
030 *
031 * <p>for \( k &gt; 0 \) the shape, \( \theta &gt; 0 \) the scale, \( \Gamma(k) \) is the gamma function
032 * and \( x \in (0, \infty) \).
033 *
034 * @see <a href="https://en.wikipedia.org/wiki/Gamma_distribution">Gamma distribution (Wikipedia)</a>
035 * @see <a href="https://mathworld.wolfram.com/GammaDistribution.html">Gamma distribution (MathWorld)</a>
036 */
037public final class GammaDistribution extends AbstractContinuousDistribution {
038    /** Support lower bound. */
039    private static final double SUPPORT_LO = 0;
040    /** Support upper bound. */
041    private static final double SUPPORT_HI = Double.POSITIVE_INFINITY;
042
043    /** The shape parameter. */
044    private final double shape;
045    /** The scale parameter. */
046    private final double scale;
047    /** Precomputed term for the log density: {@code -log(gamma(shape)) - log(scale)}. */
048    private final double minusLogGammaShapeMinusLogScale;
049    /** Cached value for inverse probability function. */
050    private final double mean;
051    /** Cached value for inverse probability function. */
052    private final double variance;
053
054    /**
055     * @param shape Shape parameter.
056     * @param scale Scale parameter.
057     */
058    private GammaDistribution(double shape,
059                              double scale) {
060        this.shape = shape;
061        this.scale = scale;
062        this.minusLogGammaShapeMinusLogScale = -LogGamma.value(shape) - Math.log(scale);
063        mean = shape * scale;
064        variance = shape * scale * scale;
065    }
066
067    /**
068     * Creates a gamma distribution.
069     *
070     * @param shape Shape parameter.
071     * @param scale Scale parameter.
072     * @return the distribution
073     * @throws IllegalArgumentException if {@code shape <= 0} or {@code scale <= 0}.
074     */
075    public static GammaDistribution of(double shape,
076                                       double scale) {
077        if (shape <= 0) {
078            throw new DistributionException(DistributionException.NOT_STRICTLY_POSITIVE, shape);
079        }
080        if (scale <= 0) {
081            throw new DistributionException(DistributionException.NOT_STRICTLY_POSITIVE, scale);
082        }
083        return new GammaDistribution(shape, scale);
084    }
085
086    /**
087     * Gets the shape parameter of this distribution.
088     *
089     * @return the shape parameter.
090     */
091    public double getShape() {
092        return shape;
093    }
094
095    /**
096     * Gets the scale parameter of this distribution.
097     *
098     * @return the scale parameter.
099     */
100    public double getScale() {
101        return scale;
102    }
103
104    /** {@inheritDoc}
105     *
106     * <p>Returns the limit when {@code x = 0}:
107     * <ul>
108     * <li>{@code shape < 1}: Infinity
109     * <li>{@code shape == 1}: 1 / scale
110     * <li>{@code shape > 1}: 0
111     * </ul>
112     */
113    @Override
114    public double density(double x) {
115        if (x <= SUPPORT_LO ||
116            x >= SUPPORT_HI) {
117            // Special case x=0
118            if (x == SUPPORT_LO && shape <= 1) {
119                return shape == 1 ?
120                    1 / scale :
121                    Double.POSITIVE_INFINITY;
122            }
123            return 0;
124        }
125
126        return RegularizedGamma.P.derivative(shape, x / scale) / scale;
127    }
128
129    /** {@inheritDoc}
130     *
131     * <p>Returns the limit when {@code x = 0}:
132     * <ul>
133     * <li>{@code shape < 1}: Infinity
134     * <li>{@code shape == 1}: -log(scale)
135     * <li>{@code shape > 1}: -Infinity
136     * </ul>
137     */
138    @Override
139    public double logDensity(double x) {
140        if (x <= SUPPORT_LO ||
141            x >= SUPPORT_HI) {
142            // Special case x=0
143            if (x == SUPPORT_LO && shape <= 1) {
144                return shape == 1 ?
145                    -Math.log(scale) :
146                    Double.POSITIVE_INFINITY;
147            }
148            return Double.NEGATIVE_INFINITY;
149        }
150
151        final double y = x / scale;
152
153        // More accurate to log the density when it is finite.
154        // See NUMBERS-174: 'Log of the Gamma P Derivative'
155        final double p = RegularizedGamma.P.derivative(shape, y) / scale;
156        if (p <= Double.MAX_VALUE && p >= Double.MIN_NORMAL) {
157            return Math.log(p);
158        }
159        // Use the log computation
160        return minusLogGammaShapeMinusLogScale - y + Math.log(y) * (shape - 1);
161    }
162
163    /** {@inheritDoc} */
164    @Override
165    public double cumulativeProbability(double x) {
166        if (x <= SUPPORT_LO) {
167            return 0;
168        } else if (x >= SUPPORT_HI) {
169            return 1;
170        }
171        return RegularizedGamma.P.value(shape, x / scale);
172    }
173
174    /** {@inheritDoc} */
175    @Override
176    public double survivalProbability(double x) {
177        if (x <= SUPPORT_LO) {
178            return 1;
179        } else if (x >= SUPPORT_HI) {
180            return 0;
181        }
182        return RegularizedGamma.Q.value(shape, x / scale);
183    }
184
185    /**
186     * {@inheritDoc}
187     *
188     * <p>For shape parameter \( k \) and scale parameter \( \theta \), the
189     * mean is \( k \theta \).
190     */
191    @Override
192    public double getMean() {
193        return mean;
194    }
195
196    /**
197     * {@inheritDoc}
198     *
199     * <p>For shape parameter \( k \) and scale parameter \( \theta \), the
200     * variance is \( k \theta^2 \).
201     */
202    @Override
203    public double getVariance() {
204        return variance;
205    }
206
207    /**
208     * {@inheritDoc}
209     *
210     * <p>The lower bound of the support is always 0.
211     *
212     * @return 0.
213     */
214    @Override
215    public double getSupportLowerBound() {
216        return SUPPORT_LO;
217    }
218
219    /**
220     * {@inheritDoc}
221     *
222     * <p>The upper bound of the support is always positive infinity.
223     *
224     * @return {@linkplain Double#POSITIVE_INFINITY positive infinity}.
225     */
226    @Override
227    public double getSupportUpperBound() {
228        return SUPPORT_HI;
229    }
230
231    /** {@inheritDoc} */
232    @Override
233    public ContinuousDistribution.Sampler createSampler(final UniformRandomProvider rng) {
234        // Gamma distribution sampler.
235        return AhrensDieterMarsagliaTsangGammaSampler.of(rng, shape, scale)::sample;
236    }
237}