View Javadoc

1   /*
2    * Copyright 2003-2004 The Apache Software Foundation.
3    *
4    * Licensed under the Apache License, Version 2.0 (the "License");
5    * you may not use this file except in compliance with the License.
6    * You may obtain a copy of the License at
7    *
8    *      http://www.apache.org/licenses/LICENSE-2.0
9    *
10   * Unless required by applicable law or agreed to in writing, software
11   * distributed under the License is distributed on an "AS IS" BASIS,
12   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13   * See the License for the specific language governing permissions and
14   * limitations under the License.
15   */
16  package org.apache.commons.math.distribution;
17  
18  import java.io.Serializable;
19  
20  import org.apache.commons.math.MathException;
21  import org.apache.commons.math.special.Gamma;
22  
23  /**
24   * The default implementation of {@link GammaDistribution}.
25   *
26   * @version $Revision: 355770 $ $Date: 2005-12-10 12:48:57 -0700 (Sat, 10 Dec 2005) $
27   */
28  public class GammaDistributionImpl extends AbstractContinuousDistribution
29      implements GammaDistribution, Serializable  {
30  
31      /** Serializable version identifier */
32      private static final long serialVersionUID = -3239549463135430361L;
33  
34      /** The shape parameter. */
35      private double alpha;
36      
37      /** The scale parameter. */
38      private double beta;
39      
40      /**
41       * Create a new gamma distribution with the given alpha and beta values.
42       * @param alpha the shape parameter.
43       * @param beta the scale parameter.
44       */
45      public GammaDistributionImpl(double alpha, double beta) {
46          super();
47          setAlpha(alpha);
48          setBeta(beta);
49      }
50      
51      /**
52       * For this disbution, X, this method returns P(X < x).
53       * 
54       * The implementation of this method is based on:
55       * <ul>
56       * <li>
57       * <a href="http://mathworld.wolfram.com/Chi-SquaredDistribution.html">
58       * Chi-Squared Distribution</a>, equation (9).</li>
59       * <li>Casella, G., & Berger, R. (1990). <i>Statistical Inference</i>.
60       * Belmont, CA: Duxbury Press.</li>
61       * </ul>
62       * 
63       * @param x the value at which the CDF is evaluated.
64       * @return CDF for this distribution. 
65       * @throws MathException if the cumulative probability can not be
66       *            computed due to convergence or other numerical errors.
67       */
68      public double cumulativeProbability(double x) throws MathException{
69          double ret;
70      
71          if (x <= 0.0) {
72              ret = 0.0;
73          } else {
74              ret = Gamma.regularizedGammaP(getAlpha(), x / getBeta());
75          }
76      
77          return ret;
78      }
79      
80      /**
81       * For this distribution, X, this method returns the critical point x, such
82       * that P(X &lt; x) = <code>p</code>.
83       * <p>
84       * Returns 0 for p=0 and <code>Double.POSITIVE_INFINITY</code> for p=1.
85       *
86       * @param p the desired probability
87       * @return x, such that P(X &lt; x) = <code>p</code>
88       * @throws MathException if the inverse cumulative probability can not be
89       *         computed due to convergence or other numerical errors.
90       * @throws IllegalArgumentException if <code>p</code> is not a valid
91       *         probability.
92       */
93      public double inverseCumulativeProbability(final double p) 
94      throws MathException {
95          if (p == 0) {
96              return 0d;
97          }
98          if (p == 1) {
99              return Double.POSITIVE_INFINITY;
100         }
101         return super.inverseCumulativeProbability(p);
102     }
103     
104     /**
105      * Modify the shape parameter, alpha.
106      * @param alpha the new shape parameter.
107      * @throws IllegalArgumentException if <code>alpha</code> is not positive.
108      */
109     public void setAlpha(double alpha) {
110         if (alpha <= 0.0) {
111             throw new IllegalArgumentException("alpha must be positive");
112         }
113         this.alpha = alpha;
114     }
115     
116     /**
117      * Access the shape parameter, alpha
118      * @return alpha.
119      */
120     public double getAlpha() {
121         return alpha;
122     }
123     
124     /**
125      * Modify the scale parameter, beta.
126      * @param beta the new scale parameter.
127      * @throws IllegalArgumentException if <code>beta</code> is not positive.
128      */
129     public void setBeta(double beta) {
130         if (beta <= 0.0) {
131             throw new IllegalArgumentException("beta must be positive");
132         }
133         this.beta = beta;
134     }
135     
136     /**
137      * Access the scale parameter, beta
138      * @return beta.
139      */
140     public double getBeta() {
141         return beta;
142     }
143     
144     /**
145      * Access the domain value lower bound, based on <code>p</code>, used to
146      * bracket a CDF root.  This method is used by
147      * {@link #inverseCumulativeProbability(double)} to find critical values.
148      * 
149      * @param p the desired probability for the critical value
150      * @return domain value lower bound, i.e.
151      *         P(X &lt; <i>lower bound</i>) &lt; <code>p</code>
152      */
153     protected double getDomainLowerBound(double p) {
154         // TODO: try to improve on this estimate
155         return Double.MIN_VALUE;
156     }
157 
158     /**
159      * Access the domain value upper bound, based on <code>p</code>, used to
160      * bracket a CDF root.  This method is used by
161      * {@link #inverseCumulativeProbability(double)} to find critical values.
162      * 
163      * @param p the desired probability for the critical value
164      * @return domain value upper bound, i.e.
165      *         P(X &lt; <i>upper bound</i>) &gt; <code>p</code> 
166      */
167     protected double getDomainUpperBound(double p) {
168         // TODO: try to improve on this estimate
169         // NOTE: gamma is skewed to the left
170         // NOTE: therefore, P(X < &mu;) > .5
171 
172         double ret;
173 
174         if (p < .5) {
175             // use mean
176             ret = getAlpha() * getBeta();
177         } else {
178             // use max value
179             ret = Double.MAX_VALUE;
180         }
181         
182         return ret;
183     }
184 
185     /**
186      * Access the initial domain value, based on <code>p</code>, used to
187      * bracket a CDF root.  This method is used by
188      * {@link #inverseCumulativeProbability(double)} to find critical values.
189      * 
190      * @param p the desired probability for the critical value
191      * @return initial domain value
192      */
193     protected double getInitialDomain(double p) {
194         // TODO: try to improve on this estimate
195         // Gamma is skewed to the left, therefore, P(X < &mu;) > .5
196 
197         double ret;
198 
199         if (p < .5) {
200             // use 1/2 mean
201             ret = getAlpha() * getBeta() * .5;
202         } else {
203             // use mean
204             ret = getAlpha() * getBeta();
205         }
206         
207         return ret;
208     }
209 }