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.exception.util.LocalizedFormats; 24 import org.apache.commons.math.special.Gamma; 25 import org.apache.commons.math.util.FastMath; 26 27 /** 28 * The default implementation of {@link GammaDistribution}. 29 * 30 * @version $Revision: 1054524 $ $Date: 2011-01-03 05:59:18 +0100 (lun. 03 janv. 2011) $ 31 */ 32 public class GammaDistributionImpl extends AbstractContinuousDistribution 33 implements GammaDistribution, Serializable { 34 35 /** 36 * Default inverse cumulative probability accuracy 37 * @since 2.1 38 */ 39 public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9; 40 41 /** Serializable version identifier */ 42 private static final long serialVersionUID = -3239549463135430361L; 43 44 /** The shape parameter. */ 45 private double alpha; 46 47 /** The scale parameter. */ 48 private double beta; 49 50 /** Inverse cumulative probability accuracy */ 51 private final double solverAbsoluteAccuracy; 52 53 /** 54 * Create a new gamma distribution with the given alpha and beta values. 55 * @param alpha the shape parameter. 56 * @param beta the scale parameter. 57 */ 58 public GammaDistributionImpl(double alpha, double beta) { 59 this(alpha, beta, DEFAULT_INVERSE_ABSOLUTE_ACCURACY); 60 } 61 62 /** 63 * Create a new gamma distribution with the given alpha and beta values. 64 * @param alpha the shape parameter. 65 * @param beta the scale parameter. 66 * @param inverseCumAccuracy the maximum absolute error in inverse cumulative probability estimates 67 * (defaults to {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}) 68 * @since 2.1 69 */ 70 public GammaDistributionImpl(double alpha, double beta, double inverseCumAccuracy) { 71 super(); 72 setAlphaInternal(alpha); 73 setBetaInternal(beta); 74 solverAbsoluteAccuracy = inverseCumAccuracy; 75 } 76 77 /** 78 * For this distribution, X, this method returns P(X < x). 79 * 80 * The implementation of this method is based on: 81 * <ul> 82 * <li> 83 * <a href="http://mathworld.wolfram.com/Chi-SquaredDistribution.html"> 84 * Chi-Squared Distribution</a>, equation (9).</li> 85 * <li>Casella, G., & Berger, R. (1990). <i>Statistical Inference</i>. 86 * Belmont, CA: Duxbury Press.</li> 87 * </ul> 88 * 89 * @param x the value at which the CDF is evaluated. 90 * @return CDF for this distribution. 91 * @throws MathException if the cumulative probability can not be 92 * computed due to convergence or other numerical errors. 93 */ 94 public double cumulativeProbability(double x) throws MathException{ 95 double ret; 96 97 if (x <= 0.0) { 98 ret = 0.0; 99 } else { 100 ret = Gamma.regularizedGammaP(alpha, x / beta); 101 } 102 103 return ret; 104 } 105 106 /** 107 * For this distribution, X, this method returns the critical point x, such 108 * that P(X < x) = <code>p</code>. 109 * <p> 110 * Returns 0 for p=0 and <code>Double.POSITIVE_INFINITY</code> for p=1.</p> 111 * 112 * @param p the desired probability 113 * @return x, such that P(X < x) = <code>p</code> 114 * @throws MathException if the inverse cumulative probability can not be 115 * computed due to convergence or other numerical errors. 116 * @throws IllegalArgumentException if <code>p</code> is not a valid 117 * probability. 118 */ 119 @Override 120 public double inverseCumulativeProbability(final double p) 121 throws MathException { 122 if (p == 0) { 123 return 0d; 124 } 125 if (p == 1) { 126 return Double.POSITIVE_INFINITY; 127 } 128 return super.inverseCumulativeProbability(p); 129 } 130 131 /** 132 * Modify the shape parameter, alpha. 133 * @param alpha the new shape parameter. 134 * @throws IllegalArgumentException if <code>alpha</code> is not positive. 135 * @deprecated as of 2.1 (class will become immutable in 3.0) 136 */ 137 @Deprecated 138 public void setAlpha(double alpha) { 139 setAlphaInternal(alpha); 140 } 141 142 /** 143 * Modify the shape parameter, alpha. 144 * @param newAlpha the new shape parameter. 145 * @throws IllegalArgumentException if <code>newAlpha</code> is not positive. 146 */ 147 private void setAlphaInternal(double newAlpha) { 148 if (newAlpha <= 0.0) { 149 throw MathRuntimeException.createIllegalArgumentException( 150 LocalizedFormats.NOT_POSITIVE_ALPHA, 151 newAlpha); 152 } 153 this.alpha = newAlpha; 154 } 155 156 /** 157 * Access the shape parameter, alpha 158 * @return alpha. 159 */ 160 public double getAlpha() { 161 return alpha; 162 } 163 164 /** 165 * Modify the scale parameter, beta. 166 * @param newBeta the new scale parameter. 167 * @throws IllegalArgumentException if <code>newBeta</code> is not positive. 168 * @deprecated as of 2.1 (class will become immutable in 3.0) 169 */ 170 @Deprecated 171 public void setBeta(double newBeta) { 172 setBetaInternal(newBeta); 173 } 174 175 /** 176 * Modify the scale parameter, beta. 177 * @param newBeta the new scale parameter. 178 * @throws IllegalArgumentException if <code>newBeta</code> is not positive. 179 */ 180 private void setBetaInternal(double newBeta) { 181 if (newBeta <= 0.0) { 182 throw MathRuntimeException.createIllegalArgumentException( 183 LocalizedFormats.NOT_POSITIVE_BETA, 184 newBeta); 185 } 186 this.beta = newBeta; 187 } 188 189 /** 190 * Access the scale parameter, beta 191 * @return beta. 192 */ 193 public double getBeta() { 194 return beta; 195 } 196 197 /** 198 * Returns the probability density for a particular point. 199 * 200 * @param x The point at which the density should be computed. 201 * @return The pdf at point x. 202 */ 203 @Override 204 public double density(double x) { 205 if (x < 0) return 0; 206 return FastMath.pow(x / beta, alpha - 1) / beta * FastMath.exp(-x / beta) / FastMath.exp(Gamma.logGamma(alpha)); 207 } 208 209 /** 210 * Return the probability density for a particular point. 211 * 212 * @param x The point at which the density should be computed. 213 * @return The pdf at point x. 214 * @deprecated 215 */ 216 @Deprecated 217 public double density(Double x) { 218 return density(x.doubleValue()); 219 } 220 221 /** 222 * Access the domain value lower bound, based on <code>p</code>, used to 223 * bracket a CDF root. This method is used by 224 * {@link #inverseCumulativeProbability(double)} to find critical values. 225 * 226 * @param p the desired probability for the critical value 227 * @return domain value lower bound, i.e. 228 * P(X < <i>lower bound</i>) < <code>p</code> 229 */ 230 @Override 231 protected double getDomainLowerBound(double p) { 232 // TODO: try to improve on this estimate 233 return Double.MIN_VALUE; 234 } 235 236 /** 237 * Access the domain value upper bound, based on <code>p</code>, used to 238 * bracket a CDF root. This method is used by 239 * {@link #inverseCumulativeProbability(double)} to find critical values. 240 * 241 * @param p the desired probability for the critical value 242 * @return domain value upper bound, i.e. 243 * P(X < <i>upper bound</i>) > <code>p</code> 244 */ 245 @Override 246 protected double getDomainUpperBound(double p) { 247 // TODO: try to improve on this estimate 248 // NOTE: gamma is skewed to the left 249 // NOTE: therefore, P(X < μ) > .5 250 251 double ret; 252 253 if (p < .5) { 254 // use mean 255 ret = alpha * beta; 256 } else { 257 // use max value 258 ret = Double.MAX_VALUE; 259 } 260 261 return ret; 262 } 263 264 /** 265 * Access the initial domain value, based on <code>p</code>, used to 266 * bracket a CDF root. This method is used by 267 * {@link #inverseCumulativeProbability(double)} to find critical values. 268 * 269 * @param p the desired probability for the critical value 270 * @return initial domain value 271 */ 272 @Override 273 protected double getInitialDomain(double p) { 274 // TODO: try to improve on this estimate 275 // Gamma is skewed to the left, therefore, P(X < μ) > .5 276 277 double ret; 278 279 if (p < .5) { 280 // use 1/2 mean 281 ret = alpha * beta * .5; 282 } else { 283 // use mean 284 ret = alpha * beta; 285 } 286 287 return ret; 288 } 289 290 /** 291 * Return the absolute accuracy setting of the solver used to estimate 292 * inverse cumulative probabilities. 293 * 294 * @return the solver absolute accuracy 295 * @since 2.1 296 */ 297 @Override 298 protected double getSolverAbsoluteAccuracy() { 299 return solverAbsoluteAccuracy; 300 } 301 302 /** 303 * Returns the upper bound of the support for the distribution. 304 * 305 * The lower bound of the support is always 0, regardless of the parameters. 306 * 307 * @return lower bound of the support (always 0) 308 * @since 2.2 309 */ 310 public double getSupportLowerBound() { 311 return 0; 312 } 313 314 /** 315 * Returns the upper bound of the support for the distribution. 316 * 317 * The upper bound of the support is always positive infinity, 318 * regardless of the parameters. 319 * 320 * @return upper bound of the support (always Double.POSITIVE_INFINITY) 321 * @since 2.2 322 */ 323 public double getSupportUpperBound() { 324 return Double.POSITIVE_INFINITY; 325 } 326 327 /** 328 * Returns the mean. 329 * 330 * For shape parameter <code>alpha</code> and scale 331 * parameter <code>beta</code>, the mean is 332 * <code>alpha * beta</code> 333 * 334 * @return the mean 335 * @since 2.2 336 */ 337 public double getNumericalMean() { 338 return getAlpha() * getBeta(); 339 } 340 341 /** 342 * Returns the variance. 343 * 344 * For shape parameter <code>alpha</code> and scale 345 * parameter <code>beta</code>, the variance is 346 * <code>alpha * beta^2</code> 347 * 348 * @return the variance 349 * @since 2.2 350 */ 351 public double getNumericalVariance() { 352 final double b = getBeta(); 353 return getAlpha() * b * b; 354 } 355 } 356
