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.RegularizedBeta;
020import org.apache.commons.numbers.gamma.LogGamma;
021import org.apache.commons.rng.UniformRandomProvider;
022import org.apache.commons.rng.sampling.distribution.ChengBetaSampler;
023
024/**
025 * Implementation of the <a href="http://en.wikipedia.org/wiki/Beta_distribution">Beta distribution</a>.
026 */
027public class BetaDistribution extends AbstractContinuousDistribution {
028    /** First shape parameter. */
029    private final double alpha;
030    /** Second shape parameter. */
031    private final double beta;
032    /** Normalizing factor used in density computations.*/
033    private final double z;
034
035    /**
036     * Creates a new instance.
037     *
038     * @param alpha First shape parameter (must be positive).
039     * @param beta Second shape parameter (must be positive).
040     */
041    public BetaDistribution(double alpha,
042                            double beta) {
043        this.alpha = alpha;
044        this.beta = beta;
045        z = LogGamma.value(alpha) + LogGamma.value(beta) - LogGamma.value(alpha + beta);
046    }
047
048    /**
049     * Access the first shape parameter, {@code alpha}.
050     *
051     * @return the first shape parameter.
052     */
053    public double getAlpha() {
054        return alpha;
055    }
056
057    /**
058     * Access the second shape parameter, {@code beta}.
059     *
060     * @return the second shape parameter.
061     */
062    public double getBeta() {
063        return beta;
064    }
065
066    /** {@inheritDoc} */
067    @Override
068    public double density(double x) {
069        final double logDensity = logDensity(x);
070        return logDensity == Double.NEGATIVE_INFINITY ? 0 : Math.exp(logDensity);
071    }
072
073    /** {@inheritDoc} **/
074    @Override
075    public double logDensity(double x) {
076        if (x < 0 ||
077            x > 1) {
078            return Double.NEGATIVE_INFINITY;
079        } else if (x == 0) {
080            if (alpha < 1) {
081                throw new DistributionException(DistributionException.TOO_SMALL,
082                                                alpha, 1);
083            }
084            return Double.NEGATIVE_INFINITY;
085        } else if (x == 1) {
086            if (beta < 1) {
087                throw new DistributionException(DistributionException.TOO_SMALL,
088                                                beta, 1);
089            }
090            return Double.NEGATIVE_INFINITY;
091        } else {
092            final double logX = Math.log(x);
093            final double log1mX = Math.log1p(-x);
094            return (alpha - 1) * logX + (beta - 1) * log1mX - z;
095        }
096    }
097
098    /** {@inheritDoc} */
099    @Override
100    public double cumulativeProbability(double x)  {
101        if (x <= 0) {
102            return 0;
103        } else if (x >= 1) {
104            return 1;
105        } else {
106            return RegularizedBeta.value(x, alpha, beta);
107        }
108    }
109
110    /**
111     * {@inheritDoc}
112     *
113     * For first shape parameter {@code alpha} and second shape parameter
114     * {@code beta}, the mean is {@code alpha / (alpha + beta)}.
115     */
116    @Override
117    public double getMean() {
118        final double a = getAlpha();
119        return a / (a + getBeta());
120    }
121
122    /**
123     * {@inheritDoc}
124     *
125     * For first shape parameter {@code alpha} and second shape parameter
126     * {@code beta}, the variance is
127     * {@code (alpha * beta) / [(alpha + beta)^2 * (alpha + beta + 1)]}.
128     */
129    @Override
130    public double getVariance() {
131        final double a = getAlpha();
132        final double b = getBeta();
133        final double alphabetasum = a + b;
134        return (a * b) / ((alphabetasum * alphabetasum) * (alphabetasum + 1));
135    }
136
137    /**
138     * {@inheritDoc}
139     *
140     * The lower bound of the support is always 0 no matter the parameters.
141     *
142     * @return lower bound of the support (always 0)
143     */
144    @Override
145    public double getSupportLowerBound() {
146        return 0;
147    }
148
149    /**
150     * {@inheritDoc}
151     *
152     * The upper bound of the support is always 1 no matter the parameters.
153     *
154     * @return upper bound of the support (always 1)
155     */
156    @Override
157    public double getSupportUpperBound() {
158        return 1;
159    }
160
161    /**
162     * {@inheritDoc}
163     *
164     * The support of this distribution is connected.
165     *
166     * @return {@code true}
167     */
168    @Override
169    public boolean isSupportConnected() {
170        return true;
171    }
172
173    /**
174     * {@inheritDoc}
175     *
176     * Sampling is performed using Cheng's algorithm:
177     * <blockquote>
178     * <pre>
179     * R. C. H. Cheng,
180     * "Generating beta variates with nonintegral shape parameters",
181     * Communications of the ACM, 21, 317-322, 1978.
182     * </pre>
183     * </blockquote>
184     */
185    @Override
186    public ContinuousDistribution.Sampler createSampler(final UniformRandomProvider rng) {
187        // Beta distribution sampler.
188        return new ChengBetaSampler(rng, alpha, beta)::sample;
189    }
190}