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: 925812 $ $Date: 2010-03-21 11:49:31 -0400 (Sun, 21 Mar 2010) $
029 */
030 public class GammaDistributionImpl extends AbstractContinuousDistribution
031 implements GammaDistribution, Serializable {
032
033 /**
034 * Default inverse cumulative probability accuracy
035 * @since 2.1
036 */
037 public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9;
038
039 /** Serializable version identifier */
040 private static final long serialVersionUID = -3239549463135430361L;
041
042 /** The shape parameter. */
043 private double alpha;
044
045 /** The scale parameter. */
046 private double beta;
047
048 /** Inverse cumulative probability accuracy */
049 private final double solverAbsoluteAccuracy;
050
051 /**
052 * Create a new gamma distribution with the given alpha and beta values.
053 * @param alpha the shape parameter.
054 * @param beta the scale parameter.
055 */
056 public GammaDistributionImpl(double alpha, double beta) {
057 this(alpha, beta, DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
058 }
059
060 /**
061 * Create a new gamma distribution with the given alpha and beta values.
062 * @param alpha the shape parameter.
063 * @param beta the scale parameter.
064 * @param inverseCumAccuracy the maximum absolute error in inverse cumulative probability estimates
065 * (defaults to {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY})
066 * @since 2.1
067 */
068 public GammaDistributionImpl(double alpha, double beta, double inverseCumAccuracy) {
069 super();
070 setAlphaInternal(alpha);
071 setBetaInternal(beta);
072 solverAbsoluteAccuracy = inverseCumAccuracy;
073 }
074
075 /**
076 * For this distribution, X, this method returns P(X < x).
077 *
078 * The implementation of this method is based on:
079 * <ul>
080 * <li>
081 * <a href="http://mathworld.wolfram.com/Chi-SquaredDistribution.html">
082 * Chi-Squared Distribution</a>, equation (9).</li>
083 * <li>Casella, G., & Berger, R. (1990). <i>Statistical Inference</i>.
084 * Belmont, CA: Duxbury Press.</li>
085 * </ul>
086 *
087 * @param x the value at which the CDF is evaluated.
088 * @return CDF for this distribution.
089 * @throws MathException if the cumulative probability can not be
090 * computed due to convergence or other numerical errors.
091 */
092 public double cumulativeProbability(double x) throws MathException{
093 double ret;
094
095 if (x <= 0.0) {
096 ret = 0.0;
097 } else {
098 ret = Gamma.regularizedGammaP(alpha, x / beta);
099 }
100
101 return ret;
102 }
103
104 /**
105 * For this distribution, X, this method returns the critical point x, such
106 * that P(X < x) = <code>p</code>.
107 * <p>
108 * Returns 0 for p=0 and <code>Double.POSITIVE_INFINITY</code> for p=1.</p>
109 *
110 * @param p the desired probability
111 * @return x, such that P(X < x) = <code>p</code>
112 * @throws MathException if the inverse cumulative probability can not be
113 * computed due to convergence or other numerical errors.
114 * @throws IllegalArgumentException if <code>p</code> is not a valid
115 * probability.
116 */
117 @Override
118 public double inverseCumulativeProbability(final double p)
119 throws MathException {
120 if (p == 0) {
121 return 0d;
122 }
123 if (p == 1) {
124 return Double.POSITIVE_INFINITY;
125 }
126 return super.inverseCumulativeProbability(p);
127 }
128
129 /**
130 * Modify the shape parameter, alpha.
131 * @param alpha the new shape parameter.
132 * @throws IllegalArgumentException if <code>alpha</code> is not positive.
133 * @deprecated as of 2.1 (class will become immutable in 3.0)
134 */
135 @Deprecated
136 public void setAlpha(double alpha) {
137 setAlphaInternal(alpha);
138 }
139
140 /**
141 * Modify the shape parameter, alpha.
142 * @param newAlpha the new shape parameter.
143 * @throws IllegalArgumentException if <code>newAlpha</code> is not positive.
144 */
145 private void setAlphaInternal(double newAlpha) {
146 if (newAlpha <= 0.0) {
147 throw MathRuntimeException.createIllegalArgumentException(
148 "alpha must be positive ({0})",
149 newAlpha);
150 }
151 this.alpha = newAlpha;
152 }
153
154 /**
155 * Access the shape parameter, alpha
156 * @return alpha.
157 */
158 public double getAlpha() {
159 return alpha;
160 }
161
162 /**
163 * Modify the scale parameter, beta.
164 * @param newBeta the new scale parameter.
165 * @throws IllegalArgumentException if <code>newBeta</code> is not positive.
166 * @deprecated as of 2.1 (class will become immutable in 3.0)
167 */
168 @Deprecated
169 public void setBeta(double newBeta) {
170 setBetaInternal(newBeta);
171 }
172
173 /**
174 * Modify the scale parameter, beta.
175 * @param newBeta the new scale parameter.
176 * @throws IllegalArgumentException if <code>newBeta</code> is not positive.
177 */
178 private void setBetaInternal(double newBeta) {
179 if (newBeta <= 0.0) {
180 throw MathRuntimeException.createIllegalArgumentException(
181 "beta must be positive ({0})",
182 newBeta);
183 }
184 this.beta = newBeta;
185 }
186
187 /**
188 * Access the scale parameter, beta
189 * @return beta.
190 */
191 public double getBeta() {
192 return beta;
193 }
194
195 /**
196 * Returns the probability density for a particular point.
197 *
198 * @param x The point at which the density should be computed.
199 * @return The pdf at point x.
200 */
201 @Override
202 public double density(double x) {
203 if (x < 0) return 0;
204 return Math.pow(x / beta, alpha - 1) / beta * Math.exp(-x / beta) / Math.exp(Gamma.logGamma(alpha));
205 }
206
207 /**
208 * Return the probability density for a particular point.
209 *
210 * @param x The point at which the density should be computed.
211 * @return The pdf at point x.
212 * @deprecated
213 */
214 public double density(Double x) {
215 return density(x.doubleValue());
216 }
217
218 /**
219 * Access the domain value lower bound, based on <code>p</code>, used to
220 * bracket a CDF root. This method is used by
221 * {@link #inverseCumulativeProbability(double)} to find critical values.
222 *
223 * @param p the desired probability for the critical value
224 * @return domain value lower bound, i.e.
225 * P(X < <i>lower bound</i>) < <code>p</code>
226 */
227 @Override
228 protected double getDomainLowerBound(double p) {
229 // TODO: try to improve on this estimate
230 return Double.MIN_VALUE;
231 }
232
233 /**
234 * Access the domain value upper bound, based on <code>p</code>, used to
235 * bracket a CDF root. This method is used by
236 * {@link #inverseCumulativeProbability(double)} to find critical values.
237 *
238 * @param p the desired probability for the critical value
239 * @return domain value upper bound, i.e.
240 * P(X < <i>upper bound</i>) > <code>p</code>
241 */
242 @Override
243 protected double getDomainUpperBound(double p) {
244 // TODO: try to improve on this estimate
245 // NOTE: gamma is skewed to the left
246 // NOTE: therefore, P(X < μ) > .5
247
248 double ret;
249
250 if (p < .5) {
251 // use mean
252 ret = alpha * beta;
253 } else {
254 // use max value
255 ret = Double.MAX_VALUE;
256 }
257
258 return ret;
259 }
260
261 /**
262 * Access the initial domain value, based on <code>p</code>, used to
263 * bracket a CDF root. This method is used by
264 * {@link #inverseCumulativeProbability(double)} to find critical values.
265 *
266 * @param p the desired probability for the critical value
267 * @return initial domain value
268 */
269 @Override
270 protected double getInitialDomain(double p) {
271 // TODO: try to improve on this estimate
272 // Gamma is skewed to the left, therefore, P(X < μ) > .5
273
274 double ret;
275
276 if (p < .5) {
277 // use 1/2 mean
278 ret = alpha * beta * .5;
279 } else {
280 // use mean
281 ret = alpha * beta;
282 }
283
284 return ret;
285 }
286
287 /**
288 * Return the absolute accuracy setting of the solver used to estimate
289 * inverse cumulative probabilities.
290 *
291 * @return the solver absolute accuracy
292 * @since 2.1
293 */
294 @Override
295 protected double getSolverAbsoluteAccuracy() {
296 return solverAbsoluteAccuracy;
297 }
298 }