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.LanczosApproximation;
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 <a href="http://en.wikipedia.org/wiki/Gamma_distribution">Gamma distribution</a>.
026 */
027public class GammaDistribution extends AbstractContinuousDistribution {
028    /** Support lower bound. */
029    private static final double SUPPORT_LO = 0;
030    /** Support upper bound. */
031    private static final double SUPPORT_HI = Double.POSITIVE_INFINITY;
032    /** Lanczos constant. */
033    private static final double LANCZOS_G = LanczosApproximation.g();
034    /** The shape parameter. */
035    private final double shape;
036    /** The scale parameter. */
037    private final double scale;
038    /**
039     * The constant value of {@code shape + g + 0.5}, where {@code g} is the
040     * Lanczos constant {@link LanczosApproximation#g()}.
041     */
042    private final double shiftedShape;
043    /**
044     * The constant value of
045     * {@code shape / scale * sqrt(e / (2 * pi * (shape + g + 0.5))) / L(shape)},
046     * where {@code L(shape)} is the Lanczos approximation returned by
047     * {@link LanczosApproximation#value(double)}. This prefactor is used in
048     * {@link #density(double)}, when no overflow occurs with the natural
049     * calculation.
050     */
051    private final double densityPrefactor1;
052    /**
053     * The constant value of
054     * {@code log(shape / scale * sqrt(e / (2 * pi * (shape + g + 0.5))) / L(shape))},
055     * where {@code L(shape)} is the Lanczos approximation returned by
056     * {@link LanczosApproximation#value(double)}. This prefactor is used in
057     * {@link #logDensity(double)}, when no overflow occurs with the natural
058     * calculation.
059     */
060    private final double logDensityPrefactor1;
061    /**
062     * The constant value of
063     * {@code shape * sqrt(e / (2 * pi * (shape + g + 0.5))) / L(shape)},
064     * where {@code L(shape)} is the Lanczos approximation returned by
065     * {@link LanczosApproximation#value(double)}. This prefactor is used in
066     * {@link #density(double)}, when overflow occurs with the natural
067     * calculation.
068     */
069    private final double densityPrefactor2;
070    /**
071     * The constant value of
072     * {@code log(shape * sqrt(e / (2 * pi * (shape + g + 0.5))) / L(shape))},
073     * where {@code L(shape)} is the Lanczos approximation returned by
074     * {@link LanczosApproximation#value(double)}. This prefactor is used in
075     * {@link #logDensity(double)}, when overflow occurs with the natural
076     * calculation.
077     */
078    private final double logDensityPrefactor2;
079    /**
080     * Lower bound on {@code y = x / scale} for the selection of the computation
081     * method in {@link #density(double)}. For {@code y <= minY}, the natural
082     * calculation overflows.
083     */
084    private final double minY;
085    /**
086     * Upper bound on {@code log(y)} ({@code y = x / scale}) for the selection
087     * of the computation method in {@link #density(double)}. For
088     * {@code log(y) >= maxLogY}, the natural calculation overflows.
089     */
090    private final double maxLogY;
091
092    /**
093     * Creates a distribution.
094     *
095     * @param shape the shape parameter
096     * @param scale the scale parameter
097     * @throws IllegalArgumentException if {@code shape <= 0} or
098     * {@code scale <= 0}.
099     */
100    public GammaDistribution(double shape,
101                             double scale) {
102        if (shape <= 0) {
103            throw new DistributionException(DistributionException.NEGATIVE, shape);
104        }
105        if (scale <= 0) {
106            throw new DistributionException(DistributionException.NEGATIVE, scale);
107        }
108
109        this.shape = shape;
110        this.scale = scale;
111        this.shiftedShape = shape + LANCZOS_G + 0.5;
112        final double aux = Math.E / (2.0 * Math.PI * shiftedShape);
113        this.densityPrefactor2 = shape * Math.sqrt(aux) / LanczosApproximation.value(shape);
114        this.logDensityPrefactor2 = Math.log(shape) + 0.5 * Math.log(aux) -
115            Math.log(LanczosApproximation.value(shape));
116        this.densityPrefactor1 = this.densityPrefactor2 / scale *
117            Math.pow(shiftedShape, -shape) *  // XXX FastMath vs Math
118            Math.exp(shape + LANCZOS_G);
119        this.logDensityPrefactor1 = this.logDensityPrefactor2 - Math.log(scale) -
120            Math.log(shiftedShape) * shape +
121            shape + LANCZOS_G;
122        this.minY = shape + LANCZOS_G - Math.log(Double.MAX_VALUE);
123        this.maxLogY = Math.log(Double.MAX_VALUE) / (shape - 1.0);
124    }
125
126    /**
127     * Returns the shape parameter of {@code this} distribution.
128     *
129     * @return the shape parameter
130     */
131    public double getShape() {
132        return shape;
133    }
134
135    /**
136     * Returns the scale parameter of {@code this} distribution.
137     *
138     * @return the scale parameter
139     */
140    public double getScale() {
141        return scale;
142    }
143
144    /** {@inheritDoc} */
145    @Override
146    public double density(double x) {
147       /* The present method must return the value of
148        *
149        *     1       x a     - x
150        * ---------- (-)  exp(---)
151        * x Gamma(a)  b        b
152        *
153        * where a is the shape parameter, and b the scale parameter.
154        * Substituting the Lanczos approximation of Gamma(a) leads to the
155        * following expression of the density
156        *
157        * a              e            1         y      a
158        * - sqrt(------------------) ---- (-----------)  exp(a - y + g),
159        * x      2 pi (a + g + 0.5)  L(a)  a + g + 0.5
160        *
161        * where y = x / b. The above formula is the "natural" computation, which
162        * is implemented when no overflow is likely to occur. If overflow occurs
163        * with the natural computation, the following identity is used. It is
164        * based on the BOOST library
165        * http://www.boost.org/doc/libs/1_35_0/libs/math/doc/sf_and_dist/html/math_toolkit/special/sf_gamma/igamma.html
166        * Formula (15) needs adaptations, which are detailed below.
167        *
168        *       y      a
169        * (-----------)  exp(a - y + g)
170        *  a + g + 0.5
171        *                              y - a - g - 0.5    y (g + 0.5)
172        *               = exp(a log1pm(---------------) - ----------- + g),
173        *                                a + g + 0.5      a + g + 0.5
174        *
175        *  where log1pm(z) = log(1 + z) - z. Therefore, the value to be
176        *  returned is
177        *
178        * a              e            1
179        * - sqrt(------------------) ----
180        * x      2 pi (a + g + 0.5)  L(a)
181        *                              y - a - g - 0.5    y (g + 0.5)
182        *               * exp(a log1pm(---------------) - ----------- + g).
183        *                                a + g + 0.5      a + g + 0.5
184        */
185        if (x <= SUPPORT_LO ||
186            x >= SUPPORT_HI) {
187            return 0;
188        }
189
190        final double y = x / scale;
191        if ((y <= minY) || (Math.log(y) >= maxLogY)) {
192            /*
193             * Overflow.
194             */
195            final double aux1 = (y - shiftedShape) / shiftedShape;
196            final double aux2 = shape * (Math.log1p(aux1) - aux1); // XXX FastMath vs Math
197            final double aux3 = -y * (LANCZOS_G + 0.5) / shiftedShape + LANCZOS_G + aux2;
198            return densityPrefactor2 / x * Math.exp(aux3);
199        }
200        /*
201         * Natural calculation.
202         */
203        return densityPrefactor1 * Math.exp(-y) * Math.pow(y, shape - 1);
204    }
205
206    /** {@inheritDoc} **/
207    @Override
208    public double logDensity(double x) {
209        /*
210         * see the comment in {@link #density(double)} for computation details
211         */
212        if (x <= SUPPORT_LO ||
213            x >= SUPPORT_HI) {
214            return Double.NEGATIVE_INFINITY;
215        }
216
217        final double y = x / scale;
218        if ((y <= minY) || (Math.log(y) >= maxLogY)) {
219            /*
220             * Overflow.
221             */
222            final double aux1 = (y - shiftedShape) / shiftedShape;
223            final double aux2 = shape * (Math.log1p(aux1) - aux1);
224            final double aux3 = -y * (LANCZOS_G + 0.5) / shiftedShape + LANCZOS_G + aux2;
225            return logDensityPrefactor2 - Math.log(x) + aux3;
226        }
227        /*
228         * Natural calculation.
229         */
230        return logDensityPrefactor1 - y + Math.log(y) * (shape - 1);
231    }
232
233    /**
234     * {@inheritDoc}
235     *
236     * The implementation of this method is based on:
237     * <ul>
238     *  <li>
239     *   <a href="http://mathworld.wolfram.com/Chi-SquaredDistribution.html">
240     *    Chi-Squared Distribution</a>, equation (9).
241     *  </li>
242     *  <li>Casella, G., &amp; Berger, R. (1990). <i>Statistical Inference</i>.
243     *    Belmont, CA: Duxbury Press.
244     *  </li>
245     * </ul>
246     */
247    @Override
248    public double cumulativeProbability(double x) {
249        if (x <= SUPPORT_LO) {
250            return 0;
251        } else if (x >= SUPPORT_HI) {
252            return 1;
253        }
254
255        return RegularizedGamma.P.value(shape, x / scale);
256    }
257
258    /**
259     * {@inheritDoc}
260     *
261     * For shape parameter {@code alpha} and scale parameter {@code beta}, the
262     * mean is {@code alpha * beta}.
263     */
264    @Override
265    public double getMean() {
266        return shape * scale;
267    }
268
269    /**
270     * {@inheritDoc}
271     *
272     * For shape parameter {@code alpha} and scale parameter {@code beta}, the
273     * variance is {@code alpha * beta^2}.
274     *
275     * @return {@inheritDoc}
276     */
277    @Override
278    public double getVariance() {
279        return shape * scale * scale;
280    }
281
282    /**
283     * {@inheritDoc}
284     *
285     * The lower bound of the support is always 0 no matter the parameters.
286     *
287     * @return lower bound of the support (always 0)
288     */
289    @Override
290    public double getSupportLowerBound() {
291        return SUPPORT_LO;
292    }
293
294    /**
295     * {@inheritDoc}
296     *
297     * The upper bound of the support is always positive infinity
298     * no matter the parameters.
299     *
300     * @return upper bound of the support (always Double.POSITIVE_INFINITY)
301     */
302    @Override
303    public double getSupportUpperBound() {
304        return SUPPORT_HI;
305    }
306
307    /**
308     * {@inheritDoc}
309     *
310     * The support of this distribution is connected.
311     *
312     * @return {@code true}
313     */
314    @Override
315    public boolean isSupportConnected() {
316        return true;
317    }
318
319    /**
320     * {@inheritDoc}
321     *
322     * <p>
323     * Sampling algorithms:
324     * <ul>
325     *  <li>
326     *   For {@code 0 < shape < 1}:
327     *   <blockquote>
328     *    Ahrens, J. H. and Dieter, U.,
329     *    <i>Computer methods for sampling from gamma, beta, Poisson and binomial distributions,</i>
330     *    Computing, 12, 223-246, 1974.
331     *   </blockquote>
332     *  </li>
333     *  <li>
334     *  For {@code shape >= 1}:
335     *   <blockquote>
336     *   Marsaglia and Tsang, <i>A Simple Method for Generating
337     *   Gamma Variables.</i> ACM Transactions on Mathematical Software,
338     *   Volume 26 Issue 3, September, 2000.
339     *   </blockquote>
340     *  </li>
341     * </ul>
342     */
343    @Override
344    public ContinuousDistribution.Sampler createSampler(final UniformRandomProvider rng) {
345        // Gamma distribution sampler.
346        return new AhrensDieterMarsagliaTsangGammaSampler(rng, shape, scale)::sample;
347    }
348}