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    }