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 org.apache.commons.math.MathException; 020 import org.apache.commons.math.MathRuntimeException; 021 import org.apache.commons.math.exception.util.LocalizedFormats; 022 import org.apache.commons.math.special.Gamma; 023 import org.apache.commons.math.special.Beta; 024 import org.apache.commons.math.util.FastMath; 025 026 /** 027 * Implements the Beta distribution. 028 * <p> 029 * References: 030 * <ul> 031 * <li><a href="http://en.wikipedia.org/wiki/Beta_distribution"> 032 * Beta distribution</a></li> 033 * </ul> 034 * </p> 035 * @version $Revision: 1054524 $ $Date: 2011-01-03 05:59:18 +0100 (lun. 03 janv. 2011) $ 036 * @since 2.0 037 */ 038 public class BetaDistributionImpl 039 extends AbstractContinuousDistribution implements BetaDistribution { 040 041 /** 042 * Default inverse cumulative probability accuracy 043 * @since 2.1 044 */ 045 public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9; 046 047 /** Serializable version identifier. */ 048 private static final long serialVersionUID = -1221965979403477668L; 049 050 /** First shape parameter. */ 051 private double alpha; 052 053 /** Second shape parameter. */ 054 private double beta; 055 056 /** Normalizing factor used in density computations. 057 * updated whenever alpha or beta are changed. 058 */ 059 private double z; 060 061 /** Inverse cumulative probability accuracy */ 062 private final double solverAbsoluteAccuracy; 063 064 /** 065 * Build a new instance. 066 * @param alpha first shape parameter (must be positive) 067 * @param beta second shape parameter (must be positive) 068 * @param inverseCumAccuracy the maximum absolute error in inverse cumulative probability estimates 069 * (defaults to {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}) 070 * @since 2.1 071 */ 072 public BetaDistributionImpl(double alpha, double beta, double inverseCumAccuracy) { 073 this.alpha = alpha; 074 this.beta = beta; 075 z = Double.NaN; 076 solverAbsoluteAccuracy = inverseCumAccuracy; 077 } 078 079 /** 080 * Build a new instance. 081 * @param alpha first shape parameter (must be positive) 082 * @param beta second shape parameter (must be positive) 083 */ 084 public BetaDistributionImpl(double alpha, double beta) { 085 this(alpha, beta, DEFAULT_INVERSE_ABSOLUTE_ACCURACY); 086 } 087 088 /** {@inheritDoc} 089 * @deprecated as of 2.1 (class will become immutable in 3.0) 090 */ 091 @Deprecated 092 public void setAlpha(double alpha) { 093 this.alpha = alpha; 094 z = Double.NaN; 095 } 096 097 /** {@inheritDoc} */ 098 public double getAlpha() { 099 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 }