PascalDistribution.java

  1. /*
  2.  * Licensed to the Apache Software Foundation (ASF) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * The ASF licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *      http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */
  17. package org.apache.commons.math3.distribution;

  19. import org.apache.commons.math3.exception.NotStrictlyPositiveException;
  20. import org.apache.commons.math3.exception.OutOfRangeException;
  21. import org.apache.commons.math3.exception.util.LocalizedFormats;
  22. import org.apache.commons.math3.random.RandomGenerator;
  23. import org.apache.commons.math3.random.Well19937c;
  24. import org.apache.commons.math3.special.Beta;
  25. import org.apache.commons.math3.util.CombinatoricsUtils;
  26. import org.apache.commons.math3.util.FastMath;

  28. /**
  29.  * <p>
  30.  * Implementation of the Pascal distribution. The Pascal distribution is a
  31.  * special case of the Negative Binomial distribution where the number of
  32.  * successes parameter is an integer.
  33.  * </p>
  34.  * <p>
  35.  * There are various ways to express the probability mass and distribution
  36.  * functions for the Pascal distribution. The present implementation represents
  37.  * the distribution of the number of failures before {@code r} successes occur.
  38.  * This is the convention adopted in e.g.
  39.  * <a href="http://mathworld.wolfram.com/NegativeBinomialDistribution.html">MathWorld</a>,
  40.  * but <em>not</em> in
  41.  * <a href="http://en.wikipedia.org/wiki/Negative_binomial_distribution">Wikipedia</a>.
  42.  * </p>
  43.  * <p>
  44.  * For a random variable {@code X} whose values are distributed according to this
  45.  * distribution, the probability mass function is given by<br/>
  46.  * {@code P(X = k) = C(k + r - 1, r - 1) * p^r * (1 - p)^k,}<br/>
  47.  * where {@code r} is the number of successes, {@code p} is the probability of
  48.  * success, and {@code X} is the total number of failures. {@code C(n, k)} is
  49.  * the binomial coefficient ({@code n} choose {@code k}). The mean and variance
  50.  * of {@code X} are<br/>
  51.  * {@code E(X) = (1 - p) * r / p, var(X) = (1 - p) * r / p^2.}<br/>
  52.  * Finally, the cumulative distribution function is given by<br/>
  53.  * {@code P(X <= k) = I(p, r, k + 1)},
  54.  * where I is the regularized incomplete Beta function.
  55.  * </p>
  56.  *
  57.  * @see <a href="http://en.wikipedia.org/wiki/Negative_binomial_distribution">
  58.  * Negative binomial distribution (Wikipedia)</a>
  59.  * @see <a href="http://mathworld.wolfram.com/NegativeBinomialDistribution.html">
  60.  * Negative binomial distribution (MathWorld)</a>
  61.  * @since 1.2 (changed to concrete class in 3.0)
  62.  */
  63. public class PascalDistribution extends AbstractIntegerDistribution {
  64.     /** Serializable version identifier. */
  65.     private static final long serialVersionUID = 6751309484392813623L;
  66.     /** The number of successes. */
  67.     private final int numberOfSuccesses;
  68.     /** The probability of success. */
  69.     private final double probabilityOfSuccess;
  70.     /** The value of {@code log(p)}, where {@code p} is the probability of success,
  71.      * stored for faster computation. */
  72.     private final double logProbabilityOfSuccess;
  73.     /** The value of {@code log(1-p)}, where {@code p} is the probability of success,
  74.      * stored for faster computation. */
  75.     private final double log1mProbabilityOfSuccess;

  77.     /**
  78.      * Create a Pascal distribution with the given number of successes and
  79.      * probability of success.
  80.      * <p>
  81.      * <b>Note:</b> this constructor will implicitly create an instance of
  82.      * {@link Well19937c} as random generator to be used for sampling only (see
  83.      * {@link #sample()} and {@link #sample(int)}). In case no sampling is
  84.      * needed for the created distribution, it is advised to pass {@code null}
  85.      * as random generator via the appropriate constructors to avoid the
  86.      * additional initialisation overhead.
  87.      *
  88.      * @param r Number of successes.
  89.      * @param p Probability of success.
  90.      * @throws NotStrictlyPositiveException if the number of successes is not positive
  91.      * @throws OutOfRangeException if the probability of success is not in the
  92.      * range {@code [0, 1]}.
  93.      */
  94.     public PascalDistribution(int r, double p)
  95.         throws NotStrictlyPositiveException, OutOfRangeException {
  96.         this(new Well19937c(), r, p);
  97.     }

  99.     /**
  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);

  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.         }

  125.         numberOfSuccesses = r;
  126.         probabilityOfSuccess = p;
  127.         logProbabilityOfSuccess = FastMath.log(p);
  128.         log1mProbabilityOfSuccess = FastMath.log1p(-p);
  129.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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. }
Created with JaCoCo 0.7.2.201409121644