    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    }