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    }