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 }