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 < 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 < 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 < <i>lower bound</i>) < <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 < <i>upper bound</i>) > <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 < μ) > .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 < μ) > .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 }