View Javadoc
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   */
17  
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;
23  
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;
30  
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;
75  
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          }
92  
93          MathUtils.checkRangeInclusive(p, 0, 1);
94  
95          numberOfSuccesses = r;
96          probabilityOfSuccess = p;
97          logProbabilityOfSuccess = FastMath.log(p);
98          log1mProbabilityOfSuccess = FastMath.log1p(-p);
99      }
100 
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     }
109 
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     }
118 
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     }
133 
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     }
148 
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     }
161 
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     }
174 
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     }
187 
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     }
199 
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     }
213 
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 }