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.  *      https://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.  */

  18. /*
  19.  * This is not the original file distributed by the Apache Software Foundation
  20.  * It has been modified by the Hipparchus project
  21.  */
  22. package org.hipparchus.distribution.discrete;

  24. import org.hipparchus.exception.LocalizedCoreFormats;
  25. import org.hipparchus.exception.MathIllegalArgumentException;
  26. import org.hipparchus.special.Beta;
  27. import org.hipparchus.util.CombinatoricsUtils;
  28. import org.hipparchus.util.FastMath;
  29. import org.hipparchus.util.MathUtils;

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

  76.     /**
  77.      * Create a Pascal distribution with the given number of successes and
  78.      * probability of success.
  79.      *
  80.      * @param r Number of successes.
  81.      * @param p Probability of success.
  82.      * @throws MathIllegalArgumentException if the number of successes is not positive
  83.      * @throws MathIllegalArgumentException if the probability of success is not in the
  84.      * range {@code [0, 1]}.
  85.      */
  86.     public PascalDistribution(int r, double p)
  87.         throws MathIllegalArgumentException {
  88.         if (r <= 0) {
  89.             throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_OF_SUCCESSES,
  90.                                                    r);
  91.         }

  93.         MathUtils.checkRangeInclusive(p, 0, 1);

  95.         numberOfSuccesses = r;
  96.         probabilityOfSuccess = p;
  97.         logProbabilityOfSuccess = FastMath.log(p);
  98.         log1mProbabilityOfSuccess = FastMath.log1p(-p);
  99.     }

  101.     /**
  102.      * Access the number of successes for this distribution.
  103.      *
  104.      * @return the number of successes.
  105.      */
  106.     public int getNumberOfSuccesses() {
  107.         return numberOfSuccesses;
  108.     }

  110.     /**
  111.      * Access the probability of success for this distribution.
  112.      *
  113.      * @return the probability of success.
  114.      */
  115.     public double getProbabilityOfSuccess() {
  116.         return probabilityOfSuccess;
  117.     }

  119.     /** {@inheritDoc} */
  120.     @Override
  121.     public double probability(int x) {
  122.         double ret;
  123.         if (x < 0) {
  124.             ret = 0.0;
  125.         } else {
  126.             ret = CombinatoricsUtils.binomialCoefficientDouble(x +
  127.                   numberOfSuccesses - 1, numberOfSuccesses - 1) *
  128.                   FastMath.pow(probabilityOfSuccess, numberOfSuccesses) *
  129.                   FastMath.pow(1.0 - probabilityOfSuccess, x);
  130.         }
  131.         return ret;
  132.     }

  134.     /** {@inheritDoc} */
  135.     @Override
  136.     public double logProbability(int x) {
  137.         double ret;
  138.         if (x < 0) {
  139.             ret = Double.NEGATIVE_INFINITY;
  140.         } else {
  141.             ret = CombinatoricsUtils.binomialCoefficientLog(x +
  142.                   numberOfSuccesses - 1, numberOfSuccesses - 1) +
  143.                   logProbabilityOfSuccess * numberOfSuccesses +
  144.                   log1mProbabilityOfSuccess * x;
  145.         }
  146.         return ret;
  147.     }

  149.     /** {@inheritDoc} */
  150.     @Override
  151.     public double cumulativeProbability(int x) {
  152.         double ret;
  153.         if (x < 0) {
  154.             ret = 0.0;
  155.         } else {
  156.             ret = Beta.regularizedBeta(probabilityOfSuccess,
  157.                     numberOfSuccesses, x + 1.0);
  158.         }
  159.         return ret;
  160.     }

  162.     /**
  163.      * {@inheritDoc}
  164.      *
  165.      * For number of successes {@code r} and probability of success {@code p},
  166.      * the mean is {@code r * (1 - p) / p}.
  167.      */
  168.     @Override
  169.     public double getNumericalMean() {
  170.         final double p = getProbabilityOfSuccess();
  171.         final double r = getNumberOfSuccesses();
  172.         return (r * (1 - p)) / p;
  173.     }

  175.     /**
  176.      * {@inheritDoc}
  177.      *
  178.      * For number of successes {@code r} and probability of success {@code p},
  179.      * the variance is {@code r * (1 - p) / p^2}.
  180.      */
  181.     @Override
  182.     public double getNumericalVariance() {
  183.         final double p = getProbabilityOfSuccess();
  184.         final double r = getNumberOfSuccesses();
  185.         return r * (1 - p) / (p * p);
  186.     }

  188.     /**
  189.      * {@inheritDoc}
  190.      *
  191.      * The lower bound of the support is always 0 no matter the parameters.
  192.      *
  193.      * @return lower bound of the support (always 0)
  194.      */
  195.     @Override
  196.     public int getSupportLowerBound() {
  197.         return 0;
  198.     }

  200.     /**
  201.      * {@inheritDoc}
  202.      *
  203.      * The upper bound of the support is always positive infinity no matter the
  204.      * parameters. Positive infinity is symbolized by {@code Integer.MAX_VALUE}.
  205.      *
  206.      * @return upper bound of the support (always {@code Integer.MAX_VALUE}
  207.      * for positive infinity)
  208.      */
  209.     @Override
  210.     public int getSupportUpperBound() {
  211.         return Integer.MAX_VALUE;
  212.     }

  214.     /**
  215.      * {@inheritDoc}
  216.      *
  217.      * The support of this distribution is connected.
  218.      *
  219.      * @return {@code true}
  220.      */
  221.     @Override
  222.     public boolean isSupportConnected() {
  223.         return true;
  224.     }
  225. }
Created with JaCoCo 0.8.12.202403310830