001: /*
002: * Copyright 2003-2004 The Apache Software Foundation.
003: *
004: * Licensed under the Apache License, Version 2.0 (the "License");
005: * you may not use this file except in compliance with the License.
006: * You may obtain a copy of the License at
007: *
008: * http://www.apache.org/licenses/LICENSE-2.0
009: *
010: * Unless required by applicable law or agreed to in writing, software
011: * distributed under the License is distributed on an "AS IS" BASIS,
012: * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
013: * See the License for the specific language governing permissions and
014: * limitations under the License.
015: */
016: package org.apache.commons.math.distribution;
017:
018: import java.io.Serializable;
019:
020: import org.apache.commons.math.MathException;
021: import org.apache.commons.math.special.Gamma;
022:
023: /**
024: * The default implementation of {@link GammaDistribution}.
025: *
026: * @version $Revision: 355770 $ $Date: 2005-12-10 12:48:57 -0700 (Sat, 10 Dec 2005) $
027: */
028: public class GammaDistributionImpl extends
029: AbstractContinuousDistribution implements GammaDistribution,
030: Serializable {
031:
032: /** Serializable version identifier */
033: private static final long serialVersionUID = -3239549463135430361L;
034:
035: /** The shape parameter. */
036: private double alpha;
037:
038: /** The scale parameter. */
039: private double beta;
040:
041: /**
042: * Create a new gamma distribution with the given alpha and beta values.
043: * @param alpha the shape parameter.
044: * @param beta the scale parameter.
045: */
046: public GammaDistributionImpl(double alpha, double beta) {
047: super ();
048: setAlpha(alpha);
049: setBeta(beta);
050: }
051:
052: /**
053: * For this disbution, X, this method returns P(X < x).
054: *
055: * The implementation of this method is based on:
056: * <ul>
057: * <li>
058: * <a href="http://mathworld.wolfram.com/Chi-SquaredDistribution.html">
059: * Chi-Squared Distribution</a>, equation (9).</li>
060: * <li>Casella, G., & Berger, R. (1990). <i>Statistical Inference</i>.
061: * Belmont, CA: Duxbury Press.</li>
062: * </ul>
063: *
064: * @param x the value at which the CDF is evaluated.
065: * @return CDF for this distribution.
066: * @throws MathException if the cumulative probability can not be
067: * computed due to convergence or other numerical errors.
068: */
069: public double cumulativeProbability(double x) throws MathException {
070: double ret;
071:
072: if (x <= 0.0) {
073: ret = 0.0;
074: } else {
075: ret = Gamma.regularizedGammaP(getAlpha(), x / getBeta());
076: }
077:
078: return ret;
079: }
080:
081: /**
082: * For this distribution, X, this method returns the critical point x, such
083: * that P(X < x) = <code>p</code>.
084: * <p>
085: * Returns 0 for p=0 and <code>Double.POSITIVE_INFINITY</code> for p=1.
086: *
087: * @param p the desired probability
088: * @return x, such that P(X < x) = <code>p</code>
089: * @throws MathException if the inverse cumulative probability can not be
090: * computed due to convergence or other numerical errors.
091: * @throws IllegalArgumentException if <code>p</code> is not a valid
092: * probability.
093: */
094: public double inverseCumulativeProbability(final double p)
095: throws MathException {
096: if (p == 0) {
097: return 0d;
098: }
099: 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 < <i>lower bound</i>) < <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 < <i>upper bound</i>) > <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 < μ) > .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 < μ) > .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: }
|