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    }