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 }