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.special.Gamma;
024    
025    /**
026     * The default implementation of {@link GammaDistribution}.
027     *
028     * @version $Revision: 772119 $ $Date: 2009-05-06 05:43:28 -0400 (Wed, 06 May 2009) $
029     */
030    public class GammaDistributionImpl extends AbstractContinuousDistribution
031        implements GammaDistribution, Serializable  {
032    
033        /** Serializable version identifier */
034        private static final long serialVersionUID = -3239549463135430361L;
035    
036        /** The shape parameter. */
037        private double alpha;
038        
039        /** The scale parameter. */
040        private double beta;
041        
042        /**
043         * Create a new gamma distribution with the given alpha and beta values.
044         * @param alpha the shape parameter.
045         * @param beta the scale parameter.
046         */
047        public GammaDistributionImpl(double alpha, double beta) {
048            super();
049            setAlpha(alpha);
050            setBeta(beta);
051        }
052        
053        /**
054         * For this distribution, X, this method returns P(X < x).
055         * 
056         * The implementation of this method is based on:
057         * <ul>
058         * <li>
059         * <a href="http://mathworld.wolfram.com/Chi-SquaredDistribution.html">
060         * Chi-Squared Distribution</a>, equation (9).</li>
061         * <li>Casella, G., & Berger, R. (1990). <i>Statistical Inference</i>.
062         * Belmont, CA: Duxbury Press.</li>
063         * </ul>
064         * 
065         * @param x the value at which the CDF is evaluated.
066         * @return CDF for this distribution. 
067         * @throws MathException if the cumulative probability can not be
068         *            computed due to convergence or other numerical errors.
069         */
070        public double cumulativeProbability(double x) throws MathException{
071            double ret;
072        
073            if (x <= 0.0) {
074                ret = 0.0;
075            } else {
076                ret = Gamma.regularizedGammaP(getAlpha(), x / getBeta());
077            }
078        
079            return ret;
080        }
081        
082        /**
083         * For this distribution, X, this method returns the critical point x, such
084         * that P(X &lt; x) = <code>p</code>.
085         * <p>
086         * Returns 0 for p=0 and <code>Double.POSITIVE_INFINITY</code> for p=1.</p>
087         *
088         * @param p the desired probability
089         * @return x, such that P(X &lt; x) = <code>p</code>
090         * @throws MathException if the inverse cumulative probability can not be
091         *         computed due to convergence or other numerical errors.
092         * @throws IllegalArgumentException if <code>p</code> is not a valid
093         *         probability.
094         */
095        @Override
096        public double inverseCumulativeProbability(final double p) 
097        throws MathException {
098            if (p == 0) {
099                return 0d;
100            }
101            if (p == 1) {
102                return Double.POSITIVE_INFINITY;
103            }
104            return super.inverseCumulativeProbability(p);
105        }
106        
107        /**
108         * Modify the shape parameter, alpha.
109         * @param alpha the new shape parameter.
110         * @throws IllegalArgumentException if <code>alpha</code> is not positive.
111         */
112        public void setAlpha(double alpha) {
113            if (alpha <= 0.0) {
114                throw MathRuntimeException.createIllegalArgumentException(
115                      "alpha must be positive ({0})",
116                      alpha);
117            }
118            this.alpha = alpha;
119        }
120        
121        /**
122         * Access the shape parameter, alpha
123         * @return alpha.
124         */
125        public double getAlpha() {
126            return alpha;
127        }
128        
129        /**
130         * Modify the scale parameter, beta.
131         * @param beta the new scale parameter.
132         * @throws IllegalArgumentException if <code>beta</code> is not positive.
133         */
134        public void setBeta(double beta) {
135            if (beta <= 0.0) {
136                throw MathRuntimeException.createIllegalArgumentException(
137                      "beta must be positive ({0})",
138                      beta);
139            }
140            this.beta = beta;
141        }
142        
143        /**
144         * Access the scale parameter, beta
145         * @return beta.
146         */
147        public double getBeta() {
148            return beta;
149        }
150    
151        /**
152         * Return the probability density for a particular point.
153         *
154         * @param x The point at which the density should be computed.
155         * @return The pdf at point x.
156         */
157        public double density(Double x) {
158            if (x < 0) return 0;
159            return Math.pow(x / getBeta(), getAlpha() - 1) / getBeta() * Math.exp(-x / getBeta()) / Math.exp(Gamma.logGamma(getAlpha()));
160        }
161    
162        /**
163         * Access the domain value lower bound, based on <code>p</code>, used to
164         * bracket a CDF root.  This method is used by
165         * {@link #inverseCumulativeProbability(double)} to find critical values.
166         * 
167         * @param p the desired probability for the critical value
168         * @return domain value lower bound, i.e.
169         *         P(X &lt; <i>lower bound</i>) &lt; <code>p</code>
170         */
171        @Override
172        protected double getDomainLowerBound(double p) {
173            // TODO: try to improve on this estimate
174            return Double.MIN_VALUE;
175        }
176    
177        /**
178         * Access the domain value upper bound, based on <code>p</code>, used to
179         * bracket a CDF root.  This method is used by
180         * {@link #inverseCumulativeProbability(double)} to find critical values.
181         * 
182         * @param p the desired probability for the critical value
183         * @return domain value upper bound, i.e.
184         *         P(X &lt; <i>upper bound</i>) &gt; <code>p</code> 
185         */
186        @Override
187        protected double getDomainUpperBound(double p) {
188            // TODO: try to improve on this estimate
189            // NOTE: gamma is skewed to the left
190            // NOTE: therefore, P(X < &mu;) > .5
191    
192            double ret;
193    
194            if (p < .5) {
195                // use mean
196                ret = getAlpha() * getBeta();
197            } else {
198                // use max value
199                ret = Double.MAX_VALUE;
200            }
201            
202            return ret;
203        }
204    
205        /**
206         * Access the initial domain value, based on <code>p</code>, used to
207         * bracket a CDF root.  This method is used by
208         * {@link #inverseCumulativeProbability(double)} to find critical values.
209         * 
210         * @param p the desired probability for the critical value
211         * @return initial domain value
212         */
213        @Override
214        protected double getInitialDomain(double p) {
215            // TODO: try to improve on this estimate
216            // Gamma is skewed to the left, therefore, P(X < &mu;) > .5
217    
218            double ret;
219    
220            if (p < .5) {
221                // use 1/2 mean
222                ret = getAlpha() * getBeta() * .5;
223            } else {
224                // use mean
225                ret = getAlpha() * getBeta();
226            }
227            
228            return ret;
229        }
230    }