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;
020
021/**
022 * Implementation of the binomial distribution.
023 *
024 * <p>The probability mass function of \( X \) is:
025 *
026 * <p>\[ f(k; n, p) = \binom{n}{k} p^k (1-p)^{n-k} \]
027 *
028 * <p>for \( n \in \{0, 1, 2, \dots\} \) the number of trials,
029 * \( p \in [0, 1] \) the probability of success,
030 * \( k \in \{0, 1, \dots, n\} \) the number of successes, and
031 *
032 * <p>\[ \binom{n}{k} = \frac{n!}{k! \, (n-k)!} \]
033 *
034 * <p>is the binomial coefficient.
035 *
036 * @see <a href="https://en.wikipedia.org/wiki/Binomial_distribution">Binomial distribution (Wikipedia)</a>
037 * @see <a href="https://mathworld.wolfram.com/BinomialDistribution.html">Binomial distribution (MathWorld)</a>
038 */
039public final class BinomialDistribution extends AbstractDiscreteDistribution {
040    /** 1/2. */
041    private static final float HALF = 0.5f;
042
043    /** The number of trials. */
044    private final int numberOfTrials;
045    /** The probability of success. */
046    private final double probabilityOfSuccess;
047    /** Cached value for pmf(x=0). */
048    private final double pmf0;
049    /** Cached value for pmf(x=n). */
050    private final double pmfn;
051
052    /**
053     * @param trials Number of trials.
054     * @param p Probability of success.
055     */
056    private BinomialDistribution(int trials,
057                                 double p) {
058        probabilityOfSuccess = p;
059        numberOfTrials = trials;
060        // Special pmf cases where the power function is more accurate:
061        //   (n choose k) == 1 for k=0, k=n
062        //   pmf = p^k (1-p)^(n-k)
063        // Note: This handles the edge case of n=0: pmf(k=0) = 1, else 0
064        if (probabilityOfSuccess >= HALF) {
065            pmf0 = Math.pow(1 - probabilityOfSuccess, numberOfTrials);
066        } else {
067            pmf0 = Math.exp(numberOfTrials * Math.log1p(-probabilityOfSuccess));
068        }
069        pmfn = Math.pow(probabilityOfSuccess, numberOfTrials);
070    }
071
072    /**
073     * Creates a binomial distribution.
074     *
075     * @param trials Number of trials.
076     * @param p Probability of success.
077     * @return the distribution
078     * @throws IllegalArgumentException if {@code trials < 0}, or if {@code p < 0}
079     * or {@code p > 1}.
080     */
081    public static BinomialDistribution of(int trials,
082                                          double p) {
083        if (trials < 0) {
084            throw new DistributionException(DistributionException.NEGATIVE,
085                                            trials);
086        }
087        ArgumentUtils.checkProbability(p);
088        // Avoid p = -0.0 to avoid returning -0.0 for some probability computations.
089        return new BinomialDistribution(trials, Math.abs(p));
090    }
091
092    /**
093     * Gets the number of trials parameter of this distribution.
094     *
095     * @return the number of trials.
096     */
097    public int getNumberOfTrials() {
098        return numberOfTrials;
099    }
100
101    /**
102     * Gets the probability of success parameter of this distribution.
103     *
104     * @return the probability of success.
105     */
106    public double getProbabilityOfSuccess() {
107        return probabilityOfSuccess;
108    }
109
110    /** {@inheritDoc} */
111    @Override
112    public double probability(int x) {
113        if (x < 0 || x > numberOfTrials) {
114            return 0;
115        } else if (x == 0) {
116            return pmf0;
117        } else if (x == numberOfTrials) {
118            return pmfn;
119        }
120        return Math.exp(SaddlePointExpansionUtils.logBinomialProbability(x,
121                        numberOfTrials, probabilityOfSuccess,
122                        1.0 - probabilityOfSuccess));
123    }
124
125    /** {@inheritDoc} **/
126    @Override
127    public double logProbability(int x) {
128        if (numberOfTrials == 0) {
129            return (x == 0) ? 0.0 : Double.NEGATIVE_INFINITY;
130        } else if (x < 0 || x > numberOfTrials) {
131            return Double.NEGATIVE_INFINITY;
132        }
133        // Special cases for x=0, x=n
134        // are handled in the saddle point expansion
135        return SaddlePointExpansionUtils.logBinomialProbability(x,
136                numberOfTrials, probabilityOfSuccess,
137                1.0 - probabilityOfSuccess);
138    }
139
140    /** {@inheritDoc} */
141    @Override
142    public double cumulativeProbability(int x) {
143        if (x < 0) {
144            return 0.0;
145        } else if (x >= numberOfTrials) {
146            return 1.0;
147        } else if (x == 0) {
148            return pmf0;
149        }
150        return RegularizedBeta.complement(probabilityOfSuccess,
151                                          x + 1.0, (double) numberOfTrials - x);
152    }
153
154    /** {@inheritDoc} */
155    @Override
156    public double survivalProbability(int x) {
157        if (x < 0) {
158            return 1.0;
159        } else if (x >= numberOfTrials) {
160            return 0.0;
161        } else if (x == numberOfTrials - 1) {
162            return pmfn;
163        }
164        return RegularizedBeta.value(probabilityOfSuccess,
165                                     x + 1.0, (double) numberOfTrials - x);
166    }
167
168    /**
169     * {@inheritDoc}
170     *
171     * <p>For number of trials \( n \) and probability of success \( p \), the mean is \( np \).
172     */
173    @Override
174    public double getMean() {
175        return numberOfTrials * probabilityOfSuccess;
176    }
177
178    /**
179     * {@inheritDoc}
180     *
181     * <p>For number of trials \( n \) and probability of success \( p \), the variance is \( np (1 - p) \).
182     */
183    @Override
184    public double getVariance() {
185        final double p = probabilityOfSuccess;
186        return numberOfTrials * p * (1 - p);
187    }
188
189    /**
190     * {@inheritDoc}
191     *
192     * <p>The lower bound of the support is always 0 except for the probability
193     * parameter {@code p = 1}.
194     *
195     * @return 0 or the number of trials.
196     */
197    @Override
198    public int getSupportLowerBound() {
199        return probabilityOfSuccess < 1.0 ? 0 : numberOfTrials;
200    }
201
202    /**
203     * {@inheritDoc}
204     *
205     * <p>The upper bound of the support is the number of trials except for the
206     * probability parameter {@code p = 0}.
207     *
208     * @return number of trials or 0.
209     */
210    @Override
211    public int getSupportUpperBound() {
212        return probabilityOfSuccess > 0.0 ? numberOfTrials : 0;
213    }
214
215    /** {@inheritDoc} */
216    @Override
217    int getMedian() {
218        // Overridden for the probability(int, int) method.
219        // This is intentionally not a public method.
220        // Can be floor or ceiling of np. For the probability in a range use the floor
221        // as this only used for values >= median+1.
222        return (int) (numberOfTrials * probabilityOfSuccess);
223    }
224}