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