    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.math3.distribution;
    018
    019    import org.apache.commons.math3.exception.NumberIsTooSmallException;
    020    import org.apache.commons.math3.exception.util.LocalizedFormats;
    021    import org.apache.commons.math3.special.Gamma;
    022    import org.apache.commons.math3.special.Beta;
    023    import org.apache.commons.math3.util.FastMath;
    024    import org.apache.commons.math3.random.RandomGenerator;
    025    import org.apache.commons.math3.random.Well19937c;
    026
    027    /**
    028     * Implements the Beta distribution.
    029     *
    030     * @see <a href="http://en.wikipedia.org/wiki/Beta_distribution">Beta distribution</a>
    031     * @version $Id: BetaDistribution.java 1416643 2012-12-03 19:37:14Z tn $
    032     * @since 2.0 (changed to concrete class in 3.0)
    033     */
    034    public class BetaDistribution extends AbstractRealDistribution {
    035        /**
    036         * Default inverse cumulative probability accuracy.
    037         * @since 2.1
    038         */
    039        public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9;
    040        /** Serializable version identifier. */
    041        private static final long serialVersionUID = -1221965979403477668L;
    042        /** First shape parameter. */
    043        private final double alpha;
    044        /** Second shape parameter. */
    045        private final double beta;
    046        /** Normalizing factor used in density computations.
    047         * updated whenever alpha or beta are changed.
    048         */
    049        private double z;
    050        /** Inverse cumulative probability accuracy. */
    051        private final double solverAbsoluteAccuracy;
    052
    053        /**
    054         * Build a new instance.
    055         *
    056         * @param alpha First shape parameter (must be positive).
    057         * @param beta Second shape parameter (must be positive).
    058         */
    059        public BetaDistribution(double alpha, double beta) {
    060            this(alpha, beta, DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
    061        }
    062
    063        /**
    064         * Build a new instance.
    065         *
    066         * @param alpha First shape parameter (must be positive).
    067         * @param beta Second shape parameter (must be positive).
    068         * @param inverseCumAccuracy Maximum absolute error in inverse
    069         * cumulative probability estimates (defaults to
    070         * {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}).
    071         * @since 2.1
    072         */
    073        public BetaDistribution(double alpha, double beta, double inverseCumAccuracy) {
    074            this(new Well19937c(), alpha, beta, inverseCumAccuracy);
    075        }
    076
    077        /**
    078         * Creates a &beta; distribution.
    079         *
    080         * @param rng Random number generator.
    081         * @param alpha First shape parameter (must be positive).
    082         * @param beta Second shape parameter (must be positive).
    083         * @param inverseCumAccuracy Maximum absolute error in inverse
    084         * cumulative probability estimates (defaults to
    085         * {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}).
    086         * @since 3.1
    087         */
    088        public BetaDistribution(RandomGenerator rng,
    089                                double alpha,
    090                                double beta,
    091                                double inverseCumAccuracy) {
    092            super(rng);
    093
    094            this.alpha = alpha;
    095            this.beta = beta;
    096            z = Double.NaN;
    097            solverAbsoluteAccuracy = inverseCumAccuracy;
    098        }
    099
    100        /**
    101         * Access the first shape parameter, {@code alpha}.
    102         *
    103         * @return the first shape parameter.
    104         */
    105        public double getAlpha() {
    106            return alpha;
    107        }
    108
    109        /**
    110         * Access the second shape parameter, {@code beta}.
    111         *
    112         * @return the second shape parameter.
    113         */
    114        public double getBeta() {
    115            return beta;
    116        }
    117
    118        /** Recompute the normalization factor. */
    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        /** {@inheritDoc} */
    126        public double density(double x) {
    127            recomputeZ();
    128            if (x < 0 || x > 1) {
    129                return 0;
    130            } else if (x == 0) {
    131                if (alpha < 1) {
    132                    throw new NumberIsTooSmallException(LocalizedFormats.CANNOT_COMPUTE_BETA_DENSITY_AT_0_FOR_SOME_ALPHA, alpha, 1, false);
    133                }
    134                return 0;
    135            } else if (x == 1) {
    136                if (beta < 1) {
    137                    throw new NumberIsTooSmallException(LocalizedFormats.CANNOT_COMPUTE_BETA_DENSITY_AT_1_FOR_SOME_BETA, beta, 1, false);
    138                }
    139                return 0;
    140            } else {
    141                double logX = FastMath.log(x);
    142                double log1mX = FastMath.log1p(-x);
    143                return FastMath.exp((alpha - 1) * logX + (beta - 1) * log1mX - z);
    144            }
    145        }
    146
    147        /** {@inheritDoc} */
    148        public double cumulativeProbability(double x)  {
    149            if (x <= 0) {
    150                return 0;
    151            } else if (x >= 1) {
    152                return 1;
    153            } else {
    154                return Beta.regularizedBeta(x, alpha, beta);
    155            }
    156        }
    157
    158        /**
    159         * Return the absolute accuracy setting of the solver used to estimate
    160         * inverse cumulative probabilities.
    161         *
    162         * @return the solver absolute accuracy.
    163         * @since 2.1
    164         */
    165        @Override
    166        protected double getSolverAbsoluteAccuracy() {
    167            return solverAbsoluteAccuracy;
    168        }
    169
    170        /**
    171         * {@inheritDoc}
    172         *
    173         * For first shape parameter {@code alpha} and second shape parameter
    174         * {@code beta}, the mean is {@code alpha / (alpha + beta)}.
    175         */
    176        public double getNumericalMean() {
    177            final double a = getAlpha();
    178            return a / (a + getBeta());
    179        }
    180
    181        /**
    182         * {@inheritDoc}
    183         *
    184         * For first shape parameter {@code alpha} and second shape parameter
    185         * {@code beta}, the variance is
    186         * {@code (alpha * beta) / [(alpha + beta)^2 * (alpha + beta + 1)]}.
    187         */
    188        public double getNumericalVariance() {
    189            final double a = getAlpha();
    190            final double b = getBeta();
    191            final double alphabetasum = a + b;
    192            return (a * b) / ((alphabetasum * alphabetasum) * (alphabetasum + 1));
    193        }
    194
    195        /**
    196         * {@inheritDoc}
    197         *
    198         * The lower bound of the support is always 0 no matter the parameters.
    199         *
    200         * @return lower bound of the support (always 0)
    201         */
    202        public double getSupportLowerBound() {
    203            return 0;
    204        }
    205
    206        /**
    207         * {@inheritDoc}
    208         *
    209         * The upper bound of the support is always 1 no matter the parameters.
    210         *
    211         * @return upper bound of the support (always 1)
    212         */
    213        public double getSupportUpperBound() {
    214            return 1;
    215        }
    216
    217        /** {@inheritDoc} */
    218        public boolean isSupportLowerBoundInclusive() {
    219            return false;
    220        }
    221
    222        /** {@inheritDoc} */
    223        public boolean isSupportUpperBoundInclusive() {
    224            return false;
    225        }
    226
    227        /**
    228         * {@inheritDoc}
    229         *
    230         * The support of this distribution is connected.
    231         *
    232         * @return {@code true}
    233         */
    234        public boolean isSupportConnected() {
    235            return true;
    236        }
    237    }