    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    }