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 org.apache.commons.math.MathException; 20 import org.apache.commons.math.MathRuntimeException; 21 import org.apache.commons.math.exception.util.LocalizedFormats; 22 import org.apache.commons.math.special.Gamma; 23 import org.apache.commons.math.special.Beta; 24 import org.apache.commons.math.util.FastMath; 25 26 /** 27 * Implements the Beta distribution. 28 * <p> 29 * References: 30 * <ul> 31 * <li><a href="http://en.wikipedia.org/wiki/Beta_distribution"> 32 * Beta distribution</a></li> 33 * </ul> 34 * </p> 35 * @version $Revision: 1054524 $ $Date: 2011-01-03 05:59:18 +0100 (lun. 03 janv. 2011) $ 36 * @since 2.0 37 */ 38 public class BetaDistributionImpl 39 extends AbstractContinuousDistribution implements BetaDistribution { 40 41 /** 42 * Default inverse cumulative probability accuracy 43 * @since 2.1 44 */ 45 public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9; 46 47 /** Serializable version identifier. */ 48 private static final long serialVersionUID = -1221965979403477668L; 49 50 /** First shape parameter. */ 51 private double alpha; 52 53 /** Second shape parameter. */ 54 private double beta; 55 56 /** Normalizing factor used in density computations. 57 * updated whenever alpha or beta are changed. 58 */ 59 private double z; 60 61 /** Inverse cumulative probability accuracy */ 62 private final double solverAbsoluteAccuracy; 63 64 /** 65 * Build a new instance. 66 * @param alpha first shape parameter (must be positive) 67 * @param beta second shape parameter (must be positive) 68 * @param inverseCumAccuracy the maximum absolute error in inverse cumulative probability estimates 69 * (defaults to {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}) 70 * @since 2.1 71 */ 72 public BetaDistributionImpl(double alpha, double beta, double inverseCumAccuracy) { 73 this.alpha = alpha; 74 this.beta = beta; 75 z = Double.NaN; 76 solverAbsoluteAccuracy = inverseCumAccuracy; 77 } 78 79 /** 80 * Build a new instance. 81 * @param alpha first shape parameter (must be positive) 82 * @param beta second shape parameter (must be positive) 83 */ 84 public BetaDistributionImpl(double alpha, double beta) { 85 this(alpha, beta, DEFAULT_INVERSE_ABSOLUTE_ACCURACY); 86 } 87 88 /** {@inheritDoc} 89 * @deprecated as of 2.1 (class will become immutable in 3.0) 90 */ 91 @Deprecated 92 public void setAlpha(double alpha) { 93 this.alpha = alpha; 94 z = Double.NaN; 95 } 96 97 /** {@inheritDoc} */ 98 public double getAlpha() { 99 return alpha; 100 } 101 102 /** {@inheritDoc} 103 * @deprecated as of 2.1 (class will become immutable in 3.0) 104 */ 105 @Deprecated 106 public void setBeta(double beta) { 107 this.beta = beta; 108 z = Double.NaN; 109 } 110 111 /** {@inheritDoc} */ 112 public double getBeta() { 113 return beta; 114 } 115 116 /** 117 * Recompute the normalization factor. 118 */ 119 private void recomputeZ() { 120 if (Double.isNaN(z)) { 121 z = Gamma.logGamma(alpha) + Gamma.logGamma(beta) - Gamma.logGamma(alpha + beta); 122 } 123 } 124 125 /** 126 * Return the probability density for a particular point. 127 * 128 * @param x The point at which the density should be computed. 129 * @return The pdf at point x. 130 * @deprecated 131 */ 132 @Deprecated 133 public double density(Double x) { 134 return density(x.doubleValue()); 135 } 136 137 /** 138 * Return the probability density for a particular point. 139 * 140 * @param x The point at which the density should be computed. 141 * @return The pdf at point x. 142 * @since 2.1 143 */ 144 @Override 145 public double density(double x) { 146 recomputeZ(); 147 if (x < 0 || x > 1) { 148 return 0; 149 } else if (x == 0) { 150 if (alpha < 1) { 151 throw MathRuntimeException.createIllegalArgumentException( 152 LocalizedFormats.CANNOT_COMPUTE_BETA_DENSITY_AT_0_FOR_SOME_ALPHA, alpha); 153 } 154 return 0; 155 } else if (x == 1) { 156 if (beta < 1) { 157 throw MathRuntimeException.createIllegalArgumentException( 158 LocalizedFormats.CANNOT_COMPUTE_BETA_DENSITY_AT_1_FOR_SOME_BETA, beta); 159 } 160 return 0; 161 } else { 162 double logX = FastMath.log(x); 163 double log1mX = FastMath.log1p(-x); 164 return FastMath.exp((alpha - 1) * logX + (beta - 1) * log1mX - z); 165 } 166 } 167 168 /** {@inheritDoc} */ 169 @Override 170 public double inverseCumulativeProbability(double p) throws MathException { 171 if (p == 0) { 172 return 0; 173 } else if (p == 1) { 174 return 1; 175 } else { 176 return super.inverseCumulativeProbability(p); 177 } 178 } 179 180 /** {@inheritDoc} */ 181 @Override 182 protected double getInitialDomain(double p) { 183 return p; 184 } 185 186 /** {@inheritDoc} */ 187 @Override 188 protected double getDomainLowerBound(double p) { 189 return 0; 190 } 191 192 /** {@inheritDoc} */ 193 @Override 194 protected double getDomainUpperBound(double p) { 195 return 1; 196 } 197 198 /** {@inheritDoc} */ 199 public double cumulativeProbability(double x) throws MathException { 200 if (x <= 0) { 201 return 0; 202 } else if (x >= 1) { 203 return 1; 204 } else { 205 return Beta.regularizedBeta(x, alpha, beta); 206 } 207 } 208 209 /** {@inheritDoc} */ 210 @Override 211 public double cumulativeProbability(double x0, double x1) throws MathException { 212 return cumulativeProbability(x1) - cumulativeProbability(x0); 213 } 214 215 /** 216 * Return the absolute accuracy setting of the solver used to estimate 217 * inverse cumulative probabilities. 218 * 219 * @return the solver absolute accuracy 220 * @since 2.1 221 */ 222 @Override 223 protected double getSolverAbsoluteAccuracy() { 224 return solverAbsoluteAccuracy; 225 } 226 227 /** 228 * Returns the lower bound of the support for this distribution. 229 * The support of the Beta distribution is always [0, 1], regardless 230 * of the parameters, so this method always returns 0. 231 * 232 * @return lower bound of the support (always 0) 233 * @since 2.2 234 */ 235 public double getSupportLowerBound() { 236 return 0; 237 } 238 239 /** 240 * Returns the upper bound of the support for this distribution. 241 * The support of the Beta distribution is always [0, 1], regardless 242 * of the parameters, so this method always returns 1. 243 * 244 * @return lower bound of the support (always 1) 245 * @since 2.2 246 */ 247 public double getSupportUpperBound() { 248 return 1; 249 } 250 251 /** 252 * Returns the mean. 253 * 254 * For first shape parameter <code>s1</code> and 255 * second shape parameter <code>s2</code>, the mean is 256 * <code>s1 / (s1 + s2)</code> 257 * 258 * @return the mean 259 * @since 2.2 260 */ 261 public double getNumericalMean() { 262 final double a = getAlpha(); 263 return a / (a + getBeta()); 264 } 265 266 /** 267 * Returns the variance. 268 * 269 * For first shape parameter <code>s1</code> and 270 * second shape parameter <code>s2</code>, 271 * the variance is 272 * <code>[ s1 * s2 ] / [ (s1 + s2)^2 * (s1 + s2 + 1) ]</code> 273 * 274 * @return the variance 275 * @since 2.2 276 */ 277 public double getNumericalVariance() { 278 final double a = getAlpha(); 279 final double b = getBeta(); 280 final double alphabetasum = a + b; 281 return (a * b) / ((alphabetasum * alphabetasum) * (alphabetasum + 1)); 282 } 283 284 } 285
