    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.math.distribution;
    018
    019    import java.io.Serializable;
    020
    021    import org.apache.commons.math.MathException;
    022    import org.apache.commons.math.MathRuntimeException;
    023    import org.apache.commons.math.exception.util.LocalizedFormats;
    024    import org.apache.commons.math.special.Gamma;
    025    import org.apache.commons.math.util.FastMath;
    026
    027    /**
    028     * The default implementation of {@link GammaDistribution}.
    029     *
    030     * @version $Revision: 1054524 $ $Date: 2011-01-03 05:59:18 +0100 (lun. 03 janv. 2011) $
    031     */
    032    public class GammaDistributionImpl extends AbstractContinuousDistribution
    033        implements GammaDistribution, Serializable  {
    034
    035        /**
    036         * Default inverse cumulative probability accuracy
    037         * @since 2.1
    038         */
    039        public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9;
    040
    041        /** Serializable version identifier */
    042        private static final long serialVersionUID = -3239549463135430361L;
    043
    044        /** The shape parameter. */
    045        private double alpha;
    046
    047        /** The scale parameter. */
    048        private double beta;
    049
    050        /** Inverse cumulative probability accuracy */
    051        private final double solverAbsoluteAccuracy;
    052
    053        /**
    054         * Create a new gamma distribution with the given alpha and beta values.
    055         * @param alpha the shape parameter.
    056         * @param beta the scale parameter.
    057         */
    058        public GammaDistributionImpl(double alpha, double beta) {
    059            this(alpha, beta, DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
    060        }
    061
    062        /**
    063         * Create a new gamma distribution with the given alpha and beta values.
    064         * @param alpha the shape parameter.
    065         * @param beta the scale parameter.
    066         * @param inverseCumAccuracy the maximum absolute error in inverse cumulative probability estimates
    067         * (defaults to {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY})
    068         * @since 2.1
    069         */
    070        public GammaDistributionImpl(double alpha, double beta, double inverseCumAccuracy) {
    071            super();
    072            setAlphaInternal(alpha);
    073            setBetaInternal(beta);
    074            solverAbsoluteAccuracy = inverseCumAccuracy;
    075        }
    076
    077        /**
    078         * For this distribution, X, this method returns P(X &lt; x).
    079         *
    080         * The implementation of this method is based on:
    081         * <ul>
    082         * <li>
    083         * <a href="http://mathworld.wolfram.com/Chi-SquaredDistribution.html">
    084         * Chi-Squared Distribution</a>, equation (9).</li>
    085         * <li>Casella, G., & Berger, R. (1990). <i>Statistical Inference</i>.
    086         * Belmont, CA: Duxbury Press.</li>
    087         * </ul>
    088         *
    089         * @param x the value at which the CDF is evaluated.
    090         * @return CDF for this distribution.
    091         * @throws MathException if the cumulative probability can not be
    092         *            computed due to convergence or other numerical errors.
    093         */
    094        public double cumulativeProbability(double x) throws MathException{
    095            double ret;
    096
    097            if (x <= 0.0) {
    098                ret = 0.0;
    099            } else {
    100                ret = Gamma.regularizedGammaP(alpha, x / beta);
    101            }
    102
    103            return ret;
    104        }
    105
    106        /**
    107         * For this distribution, X, this method returns the critical point x, such
    108         * that P(X &lt; x) = <code>p</code>.
    109         * <p>
    110         * Returns 0 for p=0 and <code>Double.POSITIVE_INFINITY</code> for p=1.</p>
    111         *
    112         * @param p the desired probability
    113         * @return x, such that P(X &lt; x) = <code>p</code>
    114         * @throws MathException if the inverse cumulative probability can not be
    115         *         computed due to convergence or other numerical errors.
    116         * @throws IllegalArgumentException if <code>p</code> is not a valid
    117         *         probability.
    118         */
    119        @Override
    120        public double inverseCumulativeProbability(final double p)
    121        throws MathException {
    122            if (p == 0) {
    123                return 0d;
    124            }
    125            if (p == 1) {
    126                return Double.POSITIVE_INFINITY;
    127            }
    128            return super.inverseCumulativeProbability(p);
    129        }
    130
    131        /**
    132         * Modify the shape parameter, alpha.
    133         * @param alpha the new shape parameter.
    134         * @throws IllegalArgumentException if <code>alpha</code> is not positive.
    135         * @deprecated as of 2.1 (class will become immutable in 3.0)
    136         */
    137        @Deprecated
    138        public void setAlpha(double alpha) {
    139            setAlphaInternal(alpha);
    140        }
    141
    142        /**
    143         * Modify the shape parameter, alpha.
    144         * @param newAlpha the new shape parameter.
    145         * @throws IllegalArgumentException if <code>newAlpha</code> is not positive.
    146         */
    147        private void setAlphaInternal(double newAlpha) {
    148            if (newAlpha <= 0.0) {
    149                throw MathRuntimeException.createIllegalArgumentException(
    150                      LocalizedFormats.NOT_POSITIVE_ALPHA,
    151                      newAlpha);
    152            }
    153            this.alpha = newAlpha;
    154        }
    155
    156        /**
    157         * Access the shape parameter, alpha
    158         * @return alpha.
    159         */
    160        public double getAlpha() {
    161            return alpha;
    162        }
    163
    164        /**
    165         * Modify the scale parameter, beta.
    166         * @param newBeta the new scale parameter.
    167         * @throws IllegalArgumentException if <code>newBeta</code> is not positive.
    168         * @deprecated as of 2.1 (class will become immutable in 3.0)
    169         */
    170        @Deprecated
    171        public void setBeta(double newBeta) {
    172            setBetaInternal(newBeta);
    173        }
    174
    175        /**
    176         * Modify the scale parameter, beta.
    177         * @param newBeta the new scale parameter.
    178         * @throws IllegalArgumentException if <code>newBeta</code> is not positive.
    179         */
    180        private void setBetaInternal(double newBeta) {
    181            if (newBeta <= 0.0) {
    182                throw MathRuntimeException.createIllegalArgumentException(
    183                      LocalizedFormats.NOT_POSITIVE_BETA,
    184                      newBeta);
    185            }
    186            this.beta = newBeta;
    187        }
    188
    189        /**
    190         * Access the scale parameter, beta
    191         * @return beta.
    192         */
    193        public double getBeta() {
    194            return beta;
    195        }
    196
    197        /**
    198         * Returns the probability density for a particular point.
    199         *
    200         * @param x The point at which the density should be computed.
    201         * @return The pdf at point x.
    202         */
    203        @Override
    204        public double density(double x) {
    205            if (x < 0) return 0;
    206            return FastMath.pow(x / beta, alpha - 1) / beta * FastMath.exp(-x / beta) / FastMath.exp(Gamma.logGamma(alpha));
    207        }
    208
    209        /**
    210         * Return the probability density for a particular point.
    211         *
    212         * @param x The point at which the density should be computed.
    213         * @return The pdf at point x.
    214         * @deprecated
    215         */
    216        @Deprecated
    217        public double density(Double x) {
    218            return density(x.doubleValue());
    219        }
    220
    221        /**
    222         * Access the domain value lower bound, based on <code>p</code>, used to
    223         * bracket a CDF root.  This method is used by
    224         * {@link #inverseCumulativeProbability(double)} to find critical values.
    225         *
    226         * @param p the desired probability for the critical value
    227         * @return domain value lower bound, i.e.
    228         *         P(X &lt; <i>lower bound</i>) &lt; <code>p</code>
    229         */
    230        @Override
    231        protected double getDomainLowerBound(double p) {
    232            // TODO: try to improve on this estimate
    233            return Double.MIN_VALUE;
    234        }
    235
    236        /**
    237         * Access the domain value upper bound, based on <code>p</code>, used to
    238         * bracket a CDF root.  This method is used by
    239         * {@link #inverseCumulativeProbability(double)} to find critical values.
    240         *
    241         * @param p the desired probability for the critical value
    242         * @return domain value upper bound, i.e.
    243         *         P(X &lt; <i>upper bound</i>) &gt; <code>p</code>
    244         */
    245        @Override
    246        protected double getDomainUpperBound(double p) {
    247            // TODO: try to improve on this estimate
    248            // NOTE: gamma is skewed to the left
    249            // NOTE: therefore, P(X < &mu;) > .5
    250
    251            double ret;
    252
    253            if (p < .5) {
    254                // use mean
    255                ret = alpha * beta;
    256            } else {
    257                // use max value
    258                ret = Double.MAX_VALUE;
    259            }
    260
    261            return ret;
    262        }
    263
    264        /**
    265         * Access the initial domain value, based on <code>p</code>, used to
    266         * bracket a CDF root.  This method is used by
    267         * {@link #inverseCumulativeProbability(double)} to find critical values.
    268         *
    269         * @param p the desired probability for the critical value
    270         * @return initial domain value
    271         */
    272        @Override
    273        protected double getInitialDomain(double p) {
    274            // TODO: try to improve on this estimate
    275            // Gamma is skewed to the left, therefore, P(X < &mu;) > .5
    276
    277            double ret;
    278
    279            if (p < .5) {
    280                // use 1/2 mean
    281                ret = alpha * beta * .5;
    282            } else {
    283                // use mean
    284                ret = alpha * beta;
    285            }
    286
    287            return ret;
    288        }
    289
    290        /**
    291         * Return the absolute accuracy setting of the solver used to estimate
    292         * inverse cumulative probabilities.
    293         *
    294         * @return the solver absolute accuracy
    295         * @since 2.1
    296         */
    297        @Override
    298        protected double getSolverAbsoluteAccuracy() {
    299            return solverAbsoluteAccuracy;
    300        }
    301
    302        /**
    303         * Returns the upper bound of the support for the distribution.
    304         *
    305         * The lower bound of the support is always 0, regardless of the parameters.
    306         *
    307         * @return lower bound of the support (always 0)
    308         * @since 2.2
    309         */
    310        public double getSupportLowerBound() {
    311            return 0;
    312        }
    313
    314        /**
    315         * Returns the upper bound of the support for the distribution.
    316         *
    317         * The upper bound of the support is always positive infinity,
    318         * regardless of the parameters.
    319         *
    320         * @return upper bound of the support (always Double.POSITIVE_INFINITY)
    321         * @since 2.2
    322         */
    323        public double getSupportUpperBound() {
    324            return Double.POSITIVE_INFINITY;
    325        }
    326
    327        /**
    328         * Returns the mean.
    329         *
    330         * For shape parameter <code>alpha</code> and scale
    331         * parameter <code>beta</code>, the mean is
    332         * <code>alpha * beta</code>
    333         *
    334         * @return the mean
    335         * @since 2.2
    336         */
    337        public double getNumericalMean() {
    338            return getAlpha() * getBeta();
    339        }
    340
    341        /**
    342         * Returns the variance.
    343         *
    344         * For shape parameter <code>alpha</code> and scale
    345         * parameter <code>beta</code>, the variance is
    346         * <code>alpha * beta^2</code>
    347         *
    348         * @return the variance
    349         * @since 2.2
    350         */
    351        public double getNumericalVariance() {
    352            final double b = getBeta();
    353            return getAlpha() * b * b;
    354        }
    355    }