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     */
017    package org.apache.commons.math3.distribution;
018    
019    import org.apache.commons.math3.exception.NotStrictlyPositiveException;
020    import org.apache.commons.math3.exception.OutOfRangeException;
021    import org.apache.commons.math3.exception.util.LocalizedFormats;
022    import org.apache.commons.math3.special.Beta;
023    import org.apache.commons.math3.util.ArithmeticUtils;
024    import org.apache.commons.math3.util.FastMath;
025    import org.apache.commons.math3.random.RandomGenerator;
026    import org.apache.commons.math3.random.Well19937c;
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     * @version $Id: PascalDistribution.java 1416643 2012-12-03 19:37:14Z tn $
062     * @since 1.2 (changed to concrete class in 3.0)
063     */
064    public class PascalDistribution extends AbstractIntegerDistribution {
065        /** Serializable version identifier. */
066        private static final long serialVersionUID = 6751309484392813623L;
067        /** The number of successes. */
068        private final int numberOfSuccesses;
069        /** The probability of success. */
070        private final double probabilityOfSuccess;
071    
072        /**
073         * Create a Pascal distribution with the given number of successes and
074         * probability of success.
075         *
076         * @param r Number of successes.
077         * @param p Probability of success.
078         * @throws NotStrictlyPositiveException if the number of successes is not positive
079         * @throws OutOfRangeException if the probability of success is not in the
080         * range {@code [0, 1]}.
081         */
082        public PascalDistribution(int r, double p)
083            throws NotStrictlyPositiveException, OutOfRangeException {
084            this(new Well19937c(), r, p);
085        }
086    
087        /**
088         * Create a Pascal distribution with the given number of successes and
089         * probability of success.
090         *
091         * @param rng Random number generator.
092         * @param r Number of successes.
093         * @param p Probability of success.
094         * @throws NotStrictlyPositiveException if the number of successes is not positive
095         * @throws OutOfRangeException if the probability of success is not in the
096         * range {@code [0, 1]}.
097         * @since 3.1
098         */
099        public PascalDistribution(RandomGenerator rng,
100                                  int r,
101                                  double p)
102            throws NotStrictlyPositiveException, OutOfRangeException {
103            super(rng);
104    
105            if (r <= 0) {
106                throw new NotStrictlyPositiveException(LocalizedFormats.NUMBER_OF_SUCCESSES,
107                                                       r);
108            }
109            if (p < 0 || p > 1) {
110                throw new OutOfRangeException(p, 0, 1);
111            }
112    
113            numberOfSuccesses = r;
114            probabilityOfSuccess = p;
115        }
116    
117        /**
118         * Access the number of successes for this distribution.
119         *
120         * @return the number of successes.
121         */
122        public int getNumberOfSuccesses() {
123            return numberOfSuccesses;
124        }
125    
126        /**
127         * Access the probability of success for this distribution.
128         *
129         * @return the probability of success.
130         */
131        public double getProbabilityOfSuccess() {
132            return probabilityOfSuccess;
133        }
134    
135        /** {@inheritDoc} */
136        public double probability(int x) {
137            double ret;
138            if (x < 0) {
139                ret = 0.0;
140            } else {
141                ret = ArithmeticUtils.binomialCoefficientDouble(x +
142                      numberOfSuccesses - 1, numberOfSuccesses - 1) *
143                      FastMath.pow(probabilityOfSuccess, numberOfSuccesses) *
144                      FastMath.pow(1.0 - probabilityOfSuccess, x);
145            }
146            return ret;
147        }
148    
149        /** {@inheritDoc} */
150        public double cumulativeProbability(int x) {
151            double ret;
152            if (x < 0) {
153                ret = 0.0;
154            } else {
155                ret = Beta.regularizedBeta(probabilityOfSuccess,
156                        numberOfSuccesses, x + 1.0);
157            }
158            return ret;
159        }
160    
161        /**
162         * {@inheritDoc}
163         *
164         * For number of successes {@code r} and probability of success {@code p},
165         * the mean is {@code r * (1 - p) / p}.
166         */
167        public double getNumericalMean() {
168            final double p = getProbabilityOfSuccess();
169            final double r = getNumberOfSuccesses();
170            return (r * (1 - p)) / p;
171        }
172    
173        /**
174         * {@inheritDoc}
175         *
176         * For number of successes {@code r} and probability of success {@code p},
177         * the variance is {@code r * (1 - p) / p^2}.
178         */
179        public double getNumericalVariance() {
180            final double p = getProbabilityOfSuccess();
181            final double r = getNumberOfSuccesses();
182            return r * (1 - p) / (p * p);
183        }
184    
185        /**
186         * {@inheritDoc}
187         *
188         * The lower bound of the support is always 0 no matter the parameters.
189         *
190         * @return lower bound of the support (always 0)
191         */
192        public int getSupportLowerBound() {
193            return 0;
194        }
195    
196        /**
197         * {@inheritDoc}
198         *
199         * The upper bound of the support is always positive infinity no matter the
200         * parameters. Positive infinity is symbolized by {@code Integer.MAX_VALUE}.
201         *
202         * @return upper bound of the support (always {@code Integer.MAX_VALUE}
203         * for positive infinity)
204         */
205        public int getSupportUpperBound() {
206            return Integer.MAX_VALUE;
207        }
208    
209        /**
210         * {@inheritDoc}
211         *
212         * The support of this distribution is connected.
213         *
214         * @return {@code true}
215         */
216        public boolean isSupportConnected() {
217            return true;
218        }
219    }