1 /* 2 * Copyright 2003-2004 The Apache Software Foundation. 3 * 4 * Licensed under the Apache License, Version 2.0 (the "License"); 5 * you may not use this file except in compliance with the License. 6 * You may obtain a copy of the License at 7 * 8 * http://www.apache.org/licenses/LICENSE-2.0 9 * 10 * Unless required by applicable law or agreed to in writing, software 11 * distributed under the License is distributed on an "AS IS" BASIS, 12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 13 * See the License for the specific language governing permissions and 14 * limitations under the License. 15 */ 16 package org.apache.commons.math.distribution; 17 18 import java.io.Serializable; 19 20 import org.apache.commons.math.MathException; 21 import org.apache.commons.math.special.Gamma; 22 23 /** 24 * The default implementation of {@link GammaDistribution}. 25 * 26 * @version $Revision: 355770 $ $Date: 2005-12-10 12:48:57 -0700 (Sat, 10 Dec 2005) $ 27 */ 28 public class GammaDistributionImpl extends AbstractContinuousDistribution 29 implements GammaDistribution, Serializable { 30 31 /** Serializable version identifier */ 32 private static final long serialVersionUID = -3239549463135430361L; 33 34 /** The shape parameter. */ 35 private double alpha; 36 37 /** The scale parameter. */ 38 private double beta; 39 40 /** 41 * Create a new gamma distribution with the given alpha and beta values. 42 * @param alpha the shape parameter. 43 * @param beta the scale parameter. 44 */ 45 public GammaDistributionImpl(double alpha, double beta) { 46 super(); 47 setAlpha(alpha); 48 setBeta(beta); 49 } 50 51 /** 52 * For this disbution, X, this method returns P(X < x). 53 * 54 * The implementation of this method is based on: 55 * <ul> 56 * <li> 57 * <a href="http://mathworld.wolfram.com/Chi-SquaredDistribution.html"> 58 * Chi-Squared Distribution</a>, equation (9).</li> 59 * <li>Casella, G., & Berger, R. (1990). <i>Statistical Inference</i>. 60 * Belmont, CA: Duxbury Press.</li> 61 * </ul> 62 * 63 * @param x the value at which the CDF is evaluated. 64 * @return CDF for this distribution. 65 * @throws MathException if the cumulative probability can not be 66 * computed due to convergence or other numerical errors. 67 */ 68 public double cumulativeProbability(double x) throws MathException{ 69 double ret; 70 71 if (x <= 0.0) { 72 ret = 0.0; 73 } else { 74 ret = Gamma.regularizedGammaP(getAlpha(), x / getBeta()); 75 } 76 77 return ret; 78 } 79 80 /** 81 * For this distribution, X, this method returns the critical point x, such 82 * that P(X < x) = <code>p</code>. 83 * <p> 84 * Returns 0 for p=0 and <code>Double.POSITIVE_INFINITY</code> for p=1. 85 * 86 * @param p the desired probability 87 * @return x, such that P(X < x) = <code>p</code> 88 * @throws MathException if the inverse cumulative probability can not be 89 * computed due to convergence or other numerical errors. 90 * @throws IllegalArgumentException if <code>p</code> is not a valid 91 * probability. 92 */ 93 public double inverseCumulativeProbability(final double p) 94 throws MathException { 95 if (p == 0) { 96 return 0d; 97 } 98 if (p == 1) { 99 return Double.POSITIVE_INFINITY; 100 } 101 return super.inverseCumulativeProbability(p); 102 } 103 104 /** 105 * Modify the shape parameter, alpha. 106 * @param alpha the new shape parameter. 107 * @throws IllegalArgumentException if <code>alpha</code> is not positive. 108 */ 109 public void setAlpha(double alpha) { 110 if (alpha <= 0.0) { 111 throw new IllegalArgumentException("alpha must be positive"); 112 } 113 this.alpha = alpha; 114 } 115 116 /** 117 * Access the shape parameter, alpha 118 * @return alpha. 119 */ 120 public double getAlpha() { 121 return alpha; 122 } 123 124 /** 125 * Modify the scale parameter, beta. 126 * @param beta the new scale parameter. 127 * @throws IllegalArgumentException if <code>beta</code> is not positive. 128 */ 129 public void setBeta(double beta) { 130 if (beta <= 0.0) { 131 throw new IllegalArgumentException("beta must be positive"); 132 } 133 this.beta = beta; 134 } 135 136 /** 137 * Access the scale parameter, beta 138 * @return beta. 139 */ 140 public double getBeta() { 141 return beta; 142 } 143 144 /** 145 * Access the domain value lower bound, based on <code>p</code>, used to 146 * bracket a CDF root. This method is used by 147 * {@link #inverseCumulativeProbability(double)} to find critical values. 148 * 149 * @param p the desired probability for the critical value 150 * @return domain value lower bound, i.e. 151 * P(X < <i>lower bound</i>) < <code>p</code> 152 */ 153 protected double getDomainLowerBound(double p) { 154 // TODO: try to improve on this estimate 155 return Double.MIN_VALUE; 156 } 157 158 /** 159 * Access the domain value upper bound, based on <code>p</code>, used to 160 * bracket a CDF root. This method is used by 161 * {@link #inverseCumulativeProbability(double)} to find critical values. 162 * 163 * @param p the desired probability for the critical value 164 * @return domain value upper bound, i.e. 165 * P(X < <i>upper bound</i>) > <code>p</code> 166 */ 167 protected double getDomainUpperBound(double p) { 168 // TODO: try to improve on this estimate 169 // NOTE: gamma is skewed to the left 170 // NOTE: therefore, P(X < μ) > .5 171 172 double ret; 173 174 if (p < .5) { 175 // use mean 176 ret = getAlpha() * getBeta(); 177 } else { 178 // use max value 179 ret = Double.MAX_VALUE; 180 } 181 182 return ret; 183 } 184 185 /** 186 * Access the initial domain value, based on <code>p</code>, used to 187 * bracket a CDF root. This method is used by 188 * {@link #inverseCumulativeProbability(double)} to find critical values. 189 * 190 * @param p the desired probability for the critical value 191 * @return initial domain value 192 */ 193 protected double getInitialDomain(double p) { 194 // TODO: try to improve on this estimate 195 // Gamma is skewed to the left, therefore, P(X < μ) > .5 196 197 double ret; 198 199 if (p < .5) { 200 // use 1/2 mean 201 ret = getAlpha() * getBeta() * .5; 202 } else { 203 // use mean 204 ret = getAlpha() * getBeta(); 205 } 206 207 return ret; 208 } 209 }