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.special.Gamma; 022 import org.apache.commons.math.special.Beta; 023 024 /** 025 * Implements the Beta distribution. 026 * <p> 027 * References: 028 * <ul> 029 * <li><a href="http://en.wikipedia.org/wiki/Beta_distribution"> 030 * Beta distribution</a></li> 031 * </ul> 032 * </p> 033 * @version $Revision: 925900 $ $Date: 2010-03-21 17:10:07 -0400 (Sun, 21 Mar 2010) $ 034 * @since 2.0 035 */ 036 public class BetaDistributionImpl 037 extends AbstractContinuousDistribution implements BetaDistribution { 038 039 /** 040 * Default inverse cumulative probability accurac 041 * @since 2.1 042 */ 043 public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9; 044 045 /** Serializable version identifier. */ 046 private static final long serialVersionUID = -1221965979403477668L; 047 048 /** First shape parameter. */ 049 private double alpha; 050 051 /** Second shape parameter. */ 052 private double beta; 053 054 /** Normalizing factor used in density computations. 055 * updated whenever alpha or beta are changed. 056 */ 057 private double z; 058 059 /** Inverse cumulative probability accuracy */ 060 private final double solverAbsoluteAccuracy; 061 062 /** 063 * Build a new instance. 064 * @param alpha first shape parameter (must be positive) 065 * @param beta second shape parameter (must be positive) 066 * @param inverseCumAccuracy the maximum absolute error in inverse cumulative probability estimates 067 * (defaults to {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}) 068 * @since 2.1 069 */ 070 public BetaDistributionImpl(double alpha, double beta, double inverseCumAccuracy) { 071 this.alpha = alpha; 072 this.beta = beta; 073 z = Double.NaN; 074 solverAbsoluteAccuracy = inverseCumAccuracy; 075 } 076 077 /** 078 * Build a new instance. 079 * @param alpha first shape parameter (must be positive) 080 * @param beta second shape parameter (must be positive) 081 */ 082 public BetaDistributionImpl(double alpha, double beta) { 083 this(alpha, beta, DEFAULT_INVERSE_ABSOLUTE_ACCURACY); 084 } 085 086 /** {@inheritDoc} 087 * @deprecated as of 2.1 (class will become immutable in 3.0) 088 */ 089 @Deprecated 090 public void setAlpha(double alpha) { 091 this.alpha = alpha; 092 z = Double.NaN; 093 } 094 095 /** {@inheritDoc} */ 096 public double getAlpha() { 097 return alpha; 098 } 099 100 /** {@inheritDoc} 101 * @deprecated as of 2.1 (class will become immutable in 3.0) 102 */ 103 @Deprecated 104 public void setBeta(double beta) { 105 this.beta = beta; 106 z = Double.NaN; 107 } 108 109 /** {@inheritDoc} */ 110 public double getBeta() { 111 return beta; 112 } 113 114 /** 115 * Recompute the normalization factor. 116 */ 117 private void recomputeZ() { 118 if (Double.isNaN(z)) { 119 z = Gamma.logGamma(alpha) + Gamma.logGamma(beta) - Gamma.logGamma(alpha + beta); 120 } 121 } 122 123 /** 124 * Return the probability density for a particular point. 125 * 126 * @param x The point at which the density should be computed. 127 * @return The pdf at point x. 128 * @deprecated 129 */ 130 public double density(Double x) { 131 return density(x.doubleValue()); 132 } 133 134 /** 135 * Return the probability density for a particular point. 136 * 137 * @param x The point at which the density should be computed. 138 * @return The pdf at point x. 139 * @since 2.1 140 */ 141 public double density(double x) { 142 recomputeZ(); 143 if (x < 0 || x > 1) { 144 return 0; 145 } else if (x == 0) { 146 if (alpha < 1) { 147 throw MathRuntimeException.createIllegalArgumentException( 148 "Cannot compute beta density at 0 when alpha = {0,number}", alpha); 149 } 150 return 0; 151 } else if (x == 1) { 152 if (beta < 1) { 153 throw MathRuntimeException.createIllegalArgumentException( 154 "Cannot compute beta density at 1 when beta = %.3g", beta); 155 } 156 return 0; 157 } else { 158 double logX = Math.log(x); 159 double log1mX = Math.log1p(-x); 160 return Math.exp((alpha - 1) * logX + (beta - 1) * log1mX - z); 161 } 162 } 163 164 /** {@inheritDoc} */ 165 @Override 166 public double inverseCumulativeProbability(double p) throws MathException { 167 if (p == 0) { 168 return 0; 169 } else if (p == 1) { 170 return 1; 171 } else { 172 return super.inverseCumulativeProbability(p); 173 } 174 } 175 176 /** {@inheritDoc} */ 177 @Override 178 protected double getInitialDomain(double p) { 179 return p; 180 } 181 182 /** {@inheritDoc} */ 183 @Override 184 protected double getDomainLowerBound(double p) { 185 return 0; 186 } 187 188 /** {@inheritDoc} */ 189 @Override 190 protected double getDomainUpperBound(double p) { 191 return 1; 192 } 193 194 /** {@inheritDoc} */ 195 public double cumulativeProbability(double x) throws MathException { 196 if (x <= 0) { 197 return 0; 198 } else if (x >= 1) { 199 return 1; 200 } else { 201 return Beta.regularizedBeta(x, alpha, beta); 202 } 203 } 204 205 /** {@inheritDoc} */ 206 @Override 207 public double cumulativeProbability(double x0, double x1) throws MathException { 208 return cumulativeProbability(x1) - cumulativeProbability(x0); 209 } 210 211 /** 212 * Return the absolute accuracy setting of the solver used to estimate 213 * inverse cumulative probabilities. 214 * 215 * @return the solver absolute accuracy 216 * @since 2.1 217 */ 218 @Override 219 protected double getSolverAbsoluteAccuracy() { 220 return solverAbsoluteAccuracy; 221 } 222 }