View Javadoc

1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  package org.apache.commons.math.distribution;
18  
19  import java.io.Serializable;
20  
21  import org.apache.commons.math.MathException;
22  import org.apache.commons.math.MathRuntimeException;
23  import org.apache.commons.math.special.Gamma;
24  
25  /**
26   * The default implementation of {@link GammaDistribution}.
27   *
28   * @version $Revision: 772119 $ $Date: 2009-05-06 05:43:28 -0400 (Wed, 06 May 2009) $
29   */
30  public class GammaDistributionImpl extends AbstractContinuousDistribution
31      implements GammaDistribution, Serializable  {
32  
33      /** Serializable version identifier */
34      private static final long serialVersionUID = -3239549463135430361L;
35  
36      /** The shape parameter. */
37      private double alpha;
38      
39      /** The scale parameter. */
40      private double beta;
41      
42      /**
43       * Create a new gamma distribution with the given alpha and beta values.
44       * @param alpha the shape parameter.
45       * @param beta the scale parameter.
46       */
47      public GammaDistributionImpl(double alpha, double beta) {
48          super();
49          setAlpha(alpha);
50          setBeta(beta);
51      }
52      
53      /**
54       * For this distribution, X, this method returns P(X < x).
55       * 
56       * The implementation of this method is based on:
57       * <ul>
58       * <li>
59       * <a href="http://mathworld.wolfram.com/Chi-SquaredDistribution.html">
60       * Chi-Squared Distribution</a>, equation (9).</li>
61       * <li>Casella, G., & Berger, R. (1990). <i>Statistical Inference</i>.
62       * Belmont, CA: Duxbury Press.</li>
63       * </ul>
64       * 
65       * @param x the value at which the CDF is evaluated.
66       * @return CDF for this distribution. 
67       * @throws MathException if the cumulative probability can not be
68       *            computed due to convergence or other numerical errors.
69       */
70      public double cumulativeProbability(double x) throws MathException{
71          double ret;
72      
73          if (x <= 0.0) {
74              ret = 0.0;
75          } else {
76              ret = Gamma.regularizedGammaP(getAlpha(), x / getBeta());
77          }
78      
79          return ret;
80      }
81      
82      /**
83       * For this distribution, X, this method returns the critical point x, such
84       * that P(X &lt; x) = <code>p</code>.
85       * <p>
86       * Returns 0 for p=0 and <code>Double.POSITIVE_INFINITY</code> for p=1.</p>
87       *
88       * @param p the desired probability
89       * @return x, such that P(X &lt; x) = <code>p</code>
90       * @throws MathException if the inverse cumulative probability can not be
91       *         computed due to convergence or other numerical errors.
92       * @throws IllegalArgumentException if <code>p</code> is not a valid
93       *         probability.
94       */
95      @Override
96      public double inverseCumulativeProbability(final double p) 
97      throws MathException {
98          if (p == 0) {
99              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 }