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 }