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.MathRuntimeException; 23 import org.apache.commons.math.special.Gamma; 24 25 /** 26 * The default implementation of {@link GammaDistribution}. 27 * 28 * @version $Revision: 772119 $ $Date: 2009-05-06 05:43:28 -0400 (Wed, 06 May 2009) $ 29 */ 30 public class GammaDistributionImpl extends AbstractContinuousDistribution 31 implements GammaDistribution, Serializable { 32 33 /** Serializable version identifier */ 34 private static final long serialVersionUID = -3239549463135430361L; 35 36 /** The shape parameter. */ 37 private double alpha; 38 39 /** The scale parameter. */ 40 private double beta; 41 42 /** 43 * Create a new gamma distribution with the given alpha and beta values. 44 * @param alpha the shape parameter. 45 * @param beta the scale parameter. 46 */ 47 public GammaDistributionImpl(double alpha, double beta) { 48 super(); 49 setAlpha(alpha); 50 setBeta(beta); 51 } 52 53 /** 54 * For this distribution, X, this method returns P(X < x). 55 * 56 * The implementation of this method is based on: 57 * <ul> 58 * <li> 59 * <a href="http://mathworld.wolfram.com/Chi-SquaredDistribution.html"> 60 * Chi-Squared Distribution</a>, equation (9).</li> 61 * <li>Casella, G., & Berger, R. (1990). <i>Statistical Inference</i>. 62 * Belmont, CA: Duxbury Press.</li> 63 * </ul> 64 * 65 * @param x the value at which the CDF is evaluated. 66 * @return CDF for this distribution. 67 * @throws MathException if the cumulative probability can not be 68 * computed due to convergence or other numerical errors. 69 */ 70 public double cumulativeProbability(double x) throws MathException{ 71 double ret; 72 73 if (x <= 0.0) { 74 ret = 0.0; 75 } else { 76 ret = Gamma.regularizedGammaP(getAlpha(), x / getBeta()); 77 } 78 79 return ret; 80 } 81 82 /** 83 * For this distribution, X, this method returns the critical point x, such 84 * that P(X < x) = <code>p</code>. 85 * <p> 86 * Returns 0 for p=0 and <code>Double.POSITIVE_INFINITY</code> for p=1.</p> 87 * 88 * @param p the desired probability 89 * @return x, such that P(X < x) = <code>p</code> 90 * @throws MathException if the inverse cumulative probability can not be 91 * computed due to convergence or other numerical errors. 92 * @throws IllegalArgumentException if <code>p</code> is not a valid 93 * probability. 94 */ 95 @Override 96 public double inverseCumulativeProbability(final double p) 97 throws MathException { 98 if (p == 0) { 99 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 }