diff options
Diffstat (limited to 'src/main/java/org/apache/commons/math3/distribution/PascalDistribution.java')
-rw-r--r-- | src/main/java/org/apache/commons/math3/distribution/PascalDistribution.java | 240 |
1 files changed, 240 insertions, 0 deletions
diff --git a/src/main/java/org/apache/commons/math3/distribution/PascalDistribution.java b/src/main/java/org/apache/commons/math3/distribution/PascalDistribution.java new file mode 100644 index 0000000..c850f8f --- /dev/null +++ b/src/main/java/org/apache/commons/math3/distribution/PascalDistribution.java @@ -0,0 +1,240 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.math3.distribution; + +import org.apache.commons.math3.exception.NotStrictlyPositiveException; +import org.apache.commons.math3.exception.OutOfRangeException; +import org.apache.commons.math3.exception.util.LocalizedFormats; +import org.apache.commons.math3.random.RandomGenerator; +import org.apache.commons.math3.random.Well19937c; +import org.apache.commons.math3.special.Beta; +import org.apache.commons.math3.util.CombinatoricsUtils; +import org.apache.commons.math3.util.FastMath; + +/** + * Implementation of the Pascal distribution. The Pascal distribution is a special case of the + * Negative Binomial distribution where the number of successes parameter is an integer. + * + * <p>There are various ways to express the probability mass and distribution functions for the + * Pascal distribution. The present implementation represents the distribution of the number of + * failures before {@code r} successes occur. This is the convention adopted in e.g. <a + * href="http://mathworld.wolfram.com/NegativeBinomialDistribution.html">MathWorld</a>, but + * <em>not</em> in <a + * href="http://en.wikipedia.org/wiki/Negative_binomial_distribution">Wikipedia</a>. + * + * <p>For a random variable {@code X} whose values are distributed according to this distribution, + * the probability mass function is given by<br> + * {@code P(X = k) = C(k + r - 1, r - 1) * p^r * (1 - p)^k,}<br> + * where {@code r} is the number of successes, {@code p} is the probability of success, and {@code + * X} is the total number of failures. {@code C(n, k)} is the binomial coefficient ({@code n} choose + * {@code k}). The mean and variance of {@code X} are<br> + * {@code E(X) = (1 - p) * r / p, var(X) = (1 - p) * r / p^2.}<br> + * Finally, the cumulative distribution function is given by<br> + * {@code P(X <= k) = I(p, r, k + 1)}, where I is the regularized incomplete Beta function. + * + * @see <a href="http://en.wikipedia.org/wiki/Negative_binomial_distribution">Negative binomial + * distribution (Wikipedia)</a> + * @see <a href="http://mathworld.wolfram.com/NegativeBinomialDistribution.html">Negative binomial + * distribution (MathWorld)</a> + * @since 1.2 (changed to concrete class in 3.0) + */ +public class PascalDistribution extends AbstractIntegerDistribution { + /** Serializable version identifier. */ + private static final long serialVersionUID = 6751309484392813623L; + + /** The number of successes. */ + private final int numberOfSuccesses; + + /** The probability of success. */ + private final double probabilityOfSuccess; + + /** + * The value of {@code log(p)}, where {@code p} is the probability of success, stored for faster + * computation. + */ + private final double logProbabilityOfSuccess; + + /** + * The value of {@code log(1-p)}, where {@code p} is the probability of success, stored for + * faster computation. + */ + private final double log1mProbabilityOfSuccess; + + /** + * Create a Pascal distribution with the given number of successes and probability of success. + * + * <p><b>Note:</b> this constructor will implicitly create an instance of {@link Well19937c} as + * random generator to be used for sampling only (see {@link #sample()} and {@link + * #sample(int)}). In case no sampling is needed for the created distribution, it is advised to + * pass {@code null} as random generator via the appropriate constructors to avoid the + * additional initialisation overhead. + * + * @param r Number of successes. + * @param p Probability of success. + * @throws NotStrictlyPositiveException if the number of successes is not positive + * @throws OutOfRangeException if the probability of success is not in the range {@code [0, 1]}. + */ + public PascalDistribution(int r, double p) + throws NotStrictlyPositiveException, OutOfRangeException { + this(new Well19937c(), r, p); + } + + /** + * Create a Pascal distribution with the given number of successes and probability of success. + * + * @param rng Random number generator. + * @param r Number of successes. + * @param p Probability of success. + * @throws NotStrictlyPositiveException if the number of successes is not positive + * @throws OutOfRangeException if the probability of success is not in the range {@code [0, 1]}. + * @since 3.1 + */ + public PascalDistribution(RandomGenerator rng, int r, double p) + throws NotStrictlyPositiveException, OutOfRangeException { + super(rng); + + if (r <= 0) { + throw new NotStrictlyPositiveException(LocalizedFormats.NUMBER_OF_SUCCESSES, r); + } + if (p < 0 || p > 1) { + throw new OutOfRangeException(p, 0, 1); + } + + numberOfSuccesses = r; + probabilityOfSuccess = p; + logProbabilityOfSuccess = FastMath.log(p); + log1mProbabilityOfSuccess = FastMath.log1p(-p); + } + + /** + * Access the number of successes for this distribution. + * + * @return the number of successes. + */ + public int getNumberOfSuccesses() { + return numberOfSuccesses; + } + + /** + * Access the probability of success for this distribution. + * + * @return the probability of success. + */ + public double getProbabilityOfSuccess() { + return probabilityOfSuccess; + } + + /** {@inheritDoc} */ + public double probability(int x) { + double ret; + if (x < 0) { + ret = 0.0; + } else { + ret = + CombinatoricsUtils.binomialCoefficientDouble( + x + numberOfSuccesses - 1, numberOfSuccesses - 1) + * FastMath.pow(probabilityOfSuccess, numberOfSuccesses) + * FastMath.pow(1.0 - probabilityOfSuccess, x); + } + return ret; + } + + /** {@inheritDoc} */ + @Override + public double logProbability(int x) { + double ret; + if (x < 0) { + ret = Double.NEGATIVE_INFINITY; + } else { + ret = + CombinatoricsUtils.binomialCoefficientLog( + x + numberOfSuccesses - 1, numberOfSuccesses - 1) + + logProbabilityOfSuccess * numberOfSuccesses + + log1mProbabilityOfSuccess * x; + } + return ret; + } + + /** {@inheritDoc} */ + public double cumulativeProbability(int x) { + double ret; + if (x < 0) { + ret = 0.0; + } else { + ret = Beta.regularizedBeta(probabilityOfSuccess, numberOfSuccesses, x + 1.0); + } + return ret; + } + + /** + * {@inheritDoc} + * + * <p>For number of successes {@code r} and probability of success {@code p}, the mean is {@code + * r * (1 - p) / p}. + */ + public double getNumericalMean() { + final double p = getProbabilityOfSuccess(); + final double r = getNumberOfSuccesses(); + return (r * (1 - p)) / p; + } + + /** + * {@inheritDoc} + * + * <p>For number of successes {@code r} and probability of success {@code p}, the variance is + * {@code r * (1 - p) / p^2}. + */ + public double getNumericalVariance() { + final double p = getProbabilityOfSuccess(); + final double r = getNumberOfSuccesses(); + return r * (1 - p) / (p * p); + } + + /** + * {@inheritDoc} + * + * <p>The lower bound of the support is always 0 no matter the parameters. + * + * @return lower bound of the support (always 0) + */ + public int getSupportLowerBound() { + return 0; + } + + /** + * {@inheritDoc} + * + * <p>The upper bound of the support is always positive infinity no matter the parameters. + * Positive infinity is symbolized by {@code Integer.MAX_VALUE}. + * + * @return upper bound of the support (always {@code Integer.MAX_VALUE} for positive infinity) + */ + public int getSupportUpperBound() { + return Integer.MAX_VALUE; + } + + /** + * {@inheritDoc} + * + * <p>The support of this distribution is connected. + * + * @return {@code true} + */ + public boolean isSupportConnected() { + return true; + } +} |