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 */
017package org.apache.commons.statistics.distribution;
018
019import org.apache.commons.numbers.gamma.LogBeta;
020import org.apache.commons.numbers.gamma.RegularizedBeta;
021import org.apache.commons.rng.UniformRandomProvider;
022import org.apache.commons.rng.sampling.distribution.ChengBetaSampler;
023
024/**
025 * Implementation of the beta distribution.
026 *
027 * <p>The probability density function of \( X \) is:
028 *
029 * <p>\[ f(x; \alpha, \beta) = \frac{1}{ B(\alpha, \beta)} x^{\alpha-1} (1-x)^{\beta-1} \]
030 *
031 * <p>for \( \alpha &gt; 0 \),
032 * \( \beta &gt; 0 \), \( x \in [0, 1] \), and
033 * the beta function, \( B \), is a normalization constant:
034 *
035 * <p>\[ B(\alpha, \beta) = \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha) \Gamma(\beta)} \]
036 *
037 * <p>where \( \Gamma \) is the gamma function.
038 *
039 * <p>\( \alpha \) and \( \beta \) are <em>shape</em> parameters.
040 *
041 * @see <a href="https://en.wikipedia.org/wiki/Beta_distribution">Beta distribution (Wikipedia)</a>
042 * @see <a href="https://mathworld.wolfram.com/BetaDistribution.html">Beta distribution (MathWorld)</a>
043 */
044public final class BetaDistribution extends AbstractContinuousDistribution {
045    /** First shape parameter. */
046    private final double alpha;
047    /** Second shape parameter. */
048    private final double beta;
049    /** Normalizing factor used in log density computations. log(beta(a, b)). */
050    private final double logBeta;
051    /** Cached value for inverse probability function. */
052    private final double mean;
053    /** Cached value for inverse probability function. */
054    private final double variance;
055
056    /**
057     * @param alpha First shape parameter (must be positive).
058     * @param beta Second shape parameter (must be positive).
059     */
060    private BetaDistribution(double alpha,
061                             double beta) {
062        this.alpha = alpha;
063        this.beta = beta;
064        logBeta = LogBeta.value(alpha, beta);
065        final double alphabetasum = alpha + beta;
066        mean = alpha / alphabetasum;
067        variance = (alpha * beta) / ((alphabetasum * alphabetasum) * (alphabetasum + 1));
068    }
069
070    /**
071     * Creates a beta distribution.
072     *
073     * @param alpha First shape parameter (must be positive).
074     * @param beta Second shape parameter (must be positive).
075     * @return the distribution
076     * @throws IllegalArgumentException if {@code alpha <= 0} or {@code beta <= 0}.
077     */
078    public static BetaDistribution of(double alpha,
079                                      double beta) {
080        if (alpha <= 0) {
081            throw new DistributionException(DistributionException.NOT_STRICTLY_POSITIVE, alpha);
082        }
083        if (beta <= 0) {
084            throw new DistributionException(DistributionException.NOT_STRICTLY_POSITIVE, beta);
085        }
086        return new BetaDistribution(alpha, beta);
087    }
088
089    /**
090     * Gets the first shape parameter of this distribution.
091     *
092     * @return the first shape parameter.
093     */
094    public double getAlpha() {
095        return alpha;
096    }
097
098    /**
099     * Gets the second shape parameter of this distribution.
100     *
101     * @return the second shape parameter.
102     */
103    public double getBeta() {
104        return beta;
105    }
106
107    /** {@inheritDoc}
108     *
109     * <p>The density is not defined when {@code x = 0, alpha < 1}, or {@code x = 1, beta < 1}.
110     * In this case the limit of infinity is returned.
111     */
112    @Override
113    public double density(double x) {
114        if (x < 0 || x > 1) {
115            return 0;
116        }
117        return RegularizedBeta.derivative(x, alpha, beta);
118    }
119
120    /** {@inheritDoc}
121     *
122     * <p>The density is not defined when {@code x = 0, alpha < 1}, or {@code x = 1, beta < 1}.
123     * In this case the limit of infinity is returned.
124     */
125    @Override
126    public double logDensity(double x) {
127        if (x < 0 || x > 1) {
128            return Double.NEGATIVE_INFINITY;
129        } else if (x == 0) {
130            if (alpha < 1) {
131                // Distribution is not valid when x=0, alpha<1
132                // due to a divide by zero error.
133                // Do not raise an exception and return the limit.
134                return Double.POSITIVE_INFINITY;
135            }
136            // Special case of cancellation: x^(a-1) (1-x)^(b-1) / B(a, b) = 1 / B(a, b)
137            if (alpha == 1) {
138                return -logBeta;
139            }
140            return Double.NEGATIVE_INFINITY;
141        } else if (x == 1) {
142            if (beta < 1) {
143                // Distribution is not valid when x=1, beta<1
144                // due to a divide by zero error.
145                // Do not raise an exception and return the limit.
146                return Double.POSITIVE_INFINITY;
147            }
148            // Special case of cancellation: x^(a-1) (1-x)^(b-1) / B(a, b) = 1 / B(a, b)
149            if (beta == 1) {
150                return -logBeta;
151            }
152            return Double.NEGATIVE_INFINITY;
153        }
154
155        // Log computation
156        final double logX = Math.log(x);
157        final double log1mX = Math.log1p(-x);
158        return (alpha - 1) * logX + (beta - 1) * log1mX - logBeta;
159    }
160
161    /** {@inheritDoc} */
162    @Override
163    public double cumulativeProbability(double x)  {
164        if (x <= 0) {
165            return 0;
166        } else if (x >= 1) {
167            return 1;
168        } else {
169            return RegularizedBeta.value(x, alpha, beta);
170        }
171    }
172
173    /** {@inheritDoc} */
174    @Override
175    public double survivalProbability(double x) {
176        if (x <= 0) {
177            return 1;
178        } else if (x >= 1) {
179            return 0;
180        } else {
181            return RegularizedBeta.complement(x, alpha, beta);
182        }
183    }
184
185    /**
186     * {@inheritDoc}
187     *
188     * <p>For first shape parameter \( \alpha \) and second shape parameter
189     * \( \beta \), the mean is:
190     *
191     * <p>\[ \frac{\alpha}{\alpha + \beta} \]
192     */
193    @Override
194    public double getMean() {
195        return mean;
196    }
197
198    /**
199     * {@inheritDoc}
200     *
201     * <p>For first shape parameter \( \alpha \) and second shape parameter
202     * \( \beta \), the variance is:
203     *
204     * <p>\[ \frac{\alpha \beta}{(\alpha + \beta)^2 (\alpha + \beta + 1)} \]
205     */
206    @Override
207    public double getVariance() {
208        return variance;
209    }
210
211    /**
212     * {@inheritDoc}
213     *
214     * <p>The lower bound of the support is always 0.
215     *
216     * @return 0.
217     */
218    @Override
219    public double getSupportLowerBound() {
220        return 0;
221    }
222
223    /**
224     * {@inheritDoc}
225     *
226     * <p>The upper bound of the support is always 1.
227     *
228     * @return 1.
229     */
230    @Override
231    public double getSupportUpperBound() {
232        return 1;
233    }
234
235    /** {@inheritDoc} */
236    @Override
237    public ContinuousDistribution.Sampler createSampler(final UniformRandomProvider rng) {
238        // Beta distribution sampler.
239        return ChengBetaSampler.of(rng, alpha, beta)::sample;
240    }
241}