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.math3.distribution;
018
019import org.apache.commons.math3.exception.NotStrictlyPositiveException;
020import org.apache.commons.math3.exception.OutOfRangeException;
021import org.apache.commons.math3.exception.util.LocalizedFormats;
022import org.apache.commons.math3.random.RandomGenerator;
023import org.apache.commons.math3.random.Well19937c;
024import org.apache.commons.math3.special.Beta;
025import org.apache.commons.math3.util.CombinatoricsUtils;
026import org.apache.commons.math3.util.FastMath;
027
028/**
029 * <p>
030 * Implementation of the Pascal distribution. The Pascal distribution is a
031 * special case of the Negative Binomial distribution where the number of
032 * successes parameter is an integer.
033 * </p>
034 * <p>
035 * There are various ways to express the probability mass and distribution
036 * functions for the Pascal distribution. The present implementation represents
037 * the distribution of the number of failures before {@code r} successes occur.
038 * This is the convention adopted in e.g.
039 * <a href="http://mathworld.wolfram.com/NegativeBinomialDistribution.html">MathWorld</a>,
040 * but <em>not</em> in
041 * <a href="http://en.wikipedia.org/wiki/Negative_binomial_distribution">Wikipedia</a>.
042 * </p>
043 * <p>
044 * For a random variable {@code X} whose values are distributed according to this
045 * distribution, the probability mass function is given by<br/>
046 * {@code P(X = k) = C(k + r - 1, r - 1) * p^r * (1 - p)^k,}<br/>
047 * where {@code r} is the number of successes, {@code p} is the probability of
048 * success, and {@code X} is the total number of failures. {@code C(n, k)} is
049 * the binomial coefficient ({@code n} choose {@code k}). The mean and variance
050 * of {@code X} are<br/>
051 * {@code E(X) = (1 - p) * r / p, var(X) = (1 - p) * r / p^2.}<br/>
052 * Finally, the cumulative distribution function is given by<br/>
053 * {@code P(X <= k) = I(p, r, k + 1)},
054 * where I is the regularized incomplete Beta function.
055 * </p>
056 *
057 * @see <a href="http://en.wikipedia.org/wiki/Negative_binomial_distribution">
058 * Negative binomial distribution (Wikipedia)</a>
059 * @see <a href="http://mathworld.wolfram.com/NegativeBinomialDistribution.html">
060 * Negative binomial distribution (MathWorld)</a>
061 * @since 1.2 (changed to concrete class in 3.0)
062 */
063public class PascalDistribution extends AbstractIntegerDistribution {
064    /** Serializable version identifier. */
065    private static final long serialVersionUID = 6751309484392813623L;
066    /** The number of successes. */
067    private final int numberOfSuccesses;
068    /** The probability of success. */
069    private final double probabilityOfSuccess;
070    /** The value of {@code log(p)}, where {@code p} is the probability of success,
071     * stored for faster computation. */
072    private final double logProbabilityOfSuccess;
073    /** The value of {@code log(1-p)}, where {@code p} is the probability of success,
074     * stored for faster computation. */
075    private final double log1mProbabilityOfSuccess;
076
077    /**
078     * Create a Pascal distribution with the given number of successes and
079     * probability of success.
080     * <p>
081     * <b>Note:</b> this constructor will implicitly create an instance of
082     * {@link Well19937c} as random generator to be used for sampling only (see
083     * {@link #sample()} and {@link #sample(int)}). In case no sampling is
084     * needed for the created distribution, it is advised to pass {@code null}
085     * as random generator via the appropriate constructors to avoid the
086     * additional initialisation overhead.
087     *
088     * @param r Number of successes.
089     * @param p Probability of success.
090     * @throws NotStrictlyPositiveException if the number of successes is not positive
091     * @throws OutOfRangeException if the probability of success is not in the
092     * range {@code [0, 1]}.
093     */
094    public PascalDistribution(int r, double p)
095        throws NotStrictlyPositiveException, OutOfRangeException {
096        this(new Well19937c(), r, p);
097    }
098
099    /**
100     * Create a Pascal distribution with the given number of successes and
101     * probability of success.
102     *
103     * @param rng Random number generator.
104     * @param r Number of successes.
105     * @param p Probability of success.
106     * @throws NotStrictlyPositiveException if the number of successes is not positive
107     * @throws OutOfRangeException if the probability of success is not in the
108     * range {@code [0, 1]}.
109     * @since 3.1
110     */
111    public PascalDistribution(RandomGenerator rng,
112                              int r,
113                              double p)
114        throws NotStrictlyPositiveException, OutOfRangeException {
115        super(rng);
116
117        if (r <= 0) {
118            throw new NotStrictlyPositiveException(LocalizedFormats.NUMBER_OF_SUCCESSES,
119                                                   r);
120        }
121        if (p < 0 || p > 1) {
122            throw new OutOfRangeException(p, 0, 1);
123        }
124
125        numberOfSuccesses = r;
126        probabilityOfSuccess = p;
127        logProbabilityOfSuccess = FastMath.log(p);
128        log1mProbabilityOfSuccess = FastMath.log1p(-p);
129    }
130
131    /**
132     * Access the number of successes for this distribution.
133     *
134     * @return the number of successes.
135     */
136    public int getNumberOfSuccesses() {
137        return numberOfSuccesses;
138    }
139
140    /**
141     * Access the probability of success for this distribution.
142     *
143     * @return the probability of success.
144     */
145    public double getProbabilityOfSuccess() {
146        return probabilityOfSuccess;
147    }
148
149    /** {@inheritDoc} */
150    public double probability(int x) {
151        double ret;
152        if (x < 0) {
153            ret = 0.0;
154        } else {
155            ret = CombinatoricsUtils.binomialCoefficientDouble(x +
156                  numberOfSuccesses - 1, numberOfSuccesses - 1) *
157                  FastMath.pow(probabilityOfSuccess, numberOfSuccesses) *
158                  FastMath.pow(1.0 - probabilityOfSuccess, x);
159        }
160        return ret;
161    }
162
163    /** {@inheritDoc} */
164    @Override
165    public double logProbability(int x) {
166        double ret;
167        if (x < 0) {
168            ret = Double.NEGATIVE_INFINITY;
169        } else {
170            ret = CombinatoricsUtils.binomialCoefficientLog(x +
171                  numberOfSuccesses - 1, numberOfSuccesses - 1) +
172                  logProbabilityOfSuccess * numberOfSuccesses +
173                  log1mProbabilityOfSuccess * x;
174        }
175        return ret;
176    }
177
178    /** {@inheritDoc} */
179    public double cumulativeProbability(int x) {
180        double ret;
181        if (x < 0) {
182            ret = 0.0;
183        } else {
184            ret = Beta.regularizedBeta(probabilityOfSuccess,
185                    numberOfSuccesses, x + 1.0);
186        }
187        return ret;
188    }
189
190    /**
191     * {@inheritDoc}
192     *
193     * For number of successes {@code r} and probability of success {@code p},
194     * the mean is {@code r * (1 - p) / p}.
195     */
196    public double getNumericalMean() {
197        final double p = getProbabilityOfSuccess();
198        final double r = getNumberOfSuccesses();
199        return (r * (1 - p)) / p;
200    }
201
202    /**
203     * {@inheritDoc}
204     *
205     * For number of successes {@code r} and probability of success {@code p},
206     * the variance is {@code r * (1 - p) / p^2}.
207     */
208    public double getNumericalVariance() {
209        final double p = getProbabilityOfSuccess();
210        final double r = getNumberOfSuccesses();
211        return r * (1 - p) / (p * p);
212    }
213
214    /**
215     * {@inheritDoc}
216     *
217     * The lower bound of the support is always 0 no matter the parameters.
218     *
219     * @return lower bound of the support (always 0)
220     */
221    public int getSupportLowerBound() {
222        return 0;
223    }
224
225    /**
226     * {@inheritDoc}
227     *
228     * The upper bound of the support is always positive infinity no matter the
229     * parameters. Positive infinity is symbolized by {@code Integer.MAX_VALUE}.
230     *
231     * @return upper bound of the support (always {@code Integer.MAX_VALUE}
232     * for positive infinity)
233     */
234    public int getSupportUpperBound() {
235        return Integer.MAX_VALUE;
236    }
237
238    /**
239     * {@inheritDoc}
240     *
241     * The support of this distribution is connected.
242     *
243     * @return {@code true}
244     */
245    public boolean isSupportConnected() {
246        return true;
247    }
248}