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