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.util.LocalizedFormats;
021 import org.apache.commons.math3.special.Gamma;
022 import org.apache.commons.math3.util.MathUtils;
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 * Implementation of the Poisson distribution.
030 *
031 * @see <a href="http://en.wikipedia.org/wiki/Poisson_distribution">Poisson distribution (Wikipedia)</a>
032 * @see <a href="http://mathworld.wolfram.com/PoissonDistribution.html">Poisson distribution (MathWorld)</a>
033 * @version $Id: PoissonDistribution.java 1416643 2012-12-03 19:37:14Z tn $
034 */
035 public class PoissonDistribution extends AbstractIntegerDistribution {
036 /**
037 * Default maximum number of iterations for cumulative probability calculations.
038 * @since 2.1
039 */
040 public static final int DEFAULT_MAX_ITERATIONS = 10000000;
041 /**
042 * Default convergence criterion.
043 * @since 2.1
044 */
045 public static final double DEFAULT_EPSILON = 1e-12;
046 /** Serializable version identifier. */
047 private static final long serialVersionUID = -3349935121172596109L;
048 /** Distribution used to compute normal approximation. */
049 private final NormalDistribution normal;
050 /** Distribution needed for the {@link #sample()} method. */
051 private final ExponentialDistribution exponential;
052 /** Mean of the distribution. */
053 private final double mean;
054
055 /**
056 * Maximum number of iterations for cumulative probability. Cumulative
057 * probabilities are estimated using either Lanczos series approximation
058 * of {@link Gamma#regularizedGammaP(double, double, double, int)}
059 * or continued fraction approximation of
060 * {@link Gamma#regularizedGammaQ(double, double, double, int)}.
061 */
062 private final int maxIterations;
063
064 /** Convergence criterion for cumulative probability. */
065 private final double epsilon;
066
067 /**
068 * Creates a new Poisson distribution with specified mean.
069 *
070 * @param p the Poisson mean
071 * @throws NotStrictlyPositiveException if {@code p <= 0}.
072 */
073 public PoissonDistribution(double p) throws NotStrictlyPositiveException {
074 this(p, DEFAULT_EPSILON, DEFAULT_MAX_ITERATIONS);
075 }
076
077 /**
078 * Creates a new Poisson distribution with specified mean, convergence
079 * criterion and maximum number of iterations.
080 *
081 * @param p Poisson mean.
082 * @param epsilon Convergence criterion for cumulative probabilities.
083 * @param maxIterations the maximum number of iterations for cumulative
084 * probabilities.
085 * @throws NotStrictlyPositiveException if {@code p <= 0}.
086 * @since 2.1
087 */
088 public PoissonDistribution(double p, double epsilon, int maxIterations)
089 throws NotStrictlyPositiveException {
090 this(new Well19937c(), p, epsilon, maxIterations);
091 }
092
093 /**
094 * Creates a new Poisson distribution with specified mean, convergence
095 * criterion and maximum number of iterations.
096 *
097 * @param rng Random number generator.
098 * @param p Poisson mean.
099 * @param epsilon Convergence criterion for cumulative probabilities.
100 * @param maxIterations the maximum number of iterations for cumulative
101 * probabilities.
102 * @throws NotStrictlyPositiveException if {@code p <= 0}.
103 * @since 3.1
104 */
105 public PoissonDistribution(RandomGenerator rng,
106 double p,
107 double epsilon,
108 int maxIterations)
109 throws NotStrictlyPositiveException {
110 super(rng);
111
112 if (p <= 0) {
113 throw new NotStrictlyPositiveException(LocalizedFormats.MEAN, p);
114 }
115 mean = p;
116 this.epsilon = epsilon;
117 this.maxIterations = maxIterations;
118
119 // Use the same RNG instance as the parent class.
120 normal = new NormalDistribution(rng, p, FastMath.sqrt(p),
121 NormalDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
122 exponential = new ExponentialDistribution(rng, 1,
123 ExponentialDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
124 }
125
126 /**
127 * Creates a new Poisson distribution with the specified mean and
128 * convergence criterion.
129 *
130 * @param p Poisson mean.
131 * @param epsilon Convergence criterion for cumulative probabilities.
132 * @throws NotStrictlyPositiveException if {@code p <= 0}.
133 * @since 2.1
134 */
135 public PoissonDistribution(double p, double epsilon)
136 throws NotStrictlyPositiveException {
137 this(p, epsilon, DEFAULT_MAX_ITERATIONS);
138 }
139
140 /**
141 * Creates a new Poisson distribution with the specified mean and maximum
142 * number of iterations.
143 *
144 * @param p Poisson mean.
145 * @param maxIterations Maximum number of iterations for cumulative
146 * probabilities.
147 * @since 2.1
148 */
149 public PoissonDistribution(double p, int maxIterations) {
150 this(p, DEFAULT_EPSILON, maxIterations);
151 }
152
153 /**
154 * Get the mean for the distribution.
155 *
156 * @return the mean for the distribution.
157 */
158 public double getMean() {
159 return mean;
160 }
161
162 /** {@inheritDoc} */
163 public double probability(int x) {
164 double ret;
165 if (x < 0 || x == Integer.MAX_VALUE) {
166 ret = 0.0;
167 } else if (x == 0) {
168 ret = FastMath.exp(-mean);
169 } else {
170 ret = FastMath.exp(-SaddlePointExpansion.getStirlingError(x) -
171 SaddlePointExpansion.getDeviancePart(x, mean)) /
172 FastMath.sqrt(MathUtils.TWO_PI * x);
173 }
174 return ret;
175 }
176
177 /** {@inheritDoc} */
178 public double cumulativeProbability(int x) {
179 if (x < 0) {
180 return 0;
181 }
182 if (x == Integer.MAX_VALUE) {
183 return 1;
184 }
185 return Gamma.regularizedGammaQ((double) x + 1, mean, epsilon,
186 maxIterations);
187 }
188
189 /**
190 * Calculates the Poisson distribution function using a normal
191 * approximation. The {@code N(mean, sqrt(mean))} distribution is used
192 * to approximate the Poisson distribution. The computation uses
193 * "half-correction" (evaluating the normal distribution function at
194 * {@code x + 0.5}).
195 *
196 * @param x Upper bound, inclusive.
197 * @return the distribution function value calculated using a normal
198 * approximation.
199 */
200 public double normalApproximateProbability(int x) {
201 // calculate the probability using half-correction
202 return normal.cumulativeProbability(x + 0.5);
203 }
204
205 /**
206 * {@inheritDoc}
207 *
208 * For mean parameter {@code p}, the mean is {@code p}.
209 */
210 public double getNumericalMean() {
211 return getMean();
212 }
213
214 /**
215 * {@inheritDoc}
216 *
217 * For mean parameter {@code p}, the variance is {@code p}.
218 */
219 public double getNumericalVariance() {
220 return getMean();
221 }
222
223 /**
224 * {@inheritDoc}
225 *
226 * The lower bound of the support is always 0 no matter the mean parameter.
227 *
228 * @return lower bound of the support (always 0)
229 */
230 public int getSupportLowerBound() {
231 return 0;
232 }
233
234 /**
235 * {@inheritDoc}
236 *
237 * The upper bound of the support is positive infinity,
238 * regardless of the parameter values. There is no integer infinity,
239 * so this method returns {@code Integer.MAX_VALUE}.
240 *
241 * @return upper bound of the support (always {@code Integer.MAX_VALUE} for
242 * positive infinity)
243 */
244 public int getSupportUpperBound() {
245 return Integer.MAX_VALUE;
246 }
247
248 /**
249 * {@inheritDoc}
250 *
251 * The support of this distribution is connected.
252 *
253 * @return {@code true}
254 */
255 public boolean isSupportConnected() {
256 return true;
257 }
258
259 /**
260 * {@inheritDoc}
261 * <p>
262 * <strong>Algorithm Description</strong>:
263 * <ul>
264 * <li>For small means, uses simulation of a Poisson process
265 * using Uniform deviates, as described
266 * <a href="http://irmi.epfl.ch/cmos/Pmmi/interactive/rng7.htm"> here</a>.
267 * The Poisson process (and hence value returned) is bounded by 1000 * mean.
268 * </li>
269 * <li>For large means, uses the rejection algorithm described in
270 * <quote>
271 * Devroye, Luc. (1981).<i>The Computer Generation of Poisson Random Variables</i>
272 * <strong>Computing</strong> vol. 26 pp. 197-207.
273 * </quote>
274 * </li>
275 * </ul>
276 * </p>
277 *
278 * @return a random value.
279 * @since 2.2
280 */
281 @Override
282 public int sample() {
283 return (int) FastMath.min(nextPoisson(mean), Integer.MAX_VALUE);
284 }
285
286 /**
287 * @param meanPoisson Mean of the Poisson distribution.
288 * @return the next sample.
289 */
290 private long nextPoisson(double meanPoisson) {
291 final double pivot = 40.0d;
292 if (meanPoisson < pivot) {
293 double p = FastMath.exp(-meanPoisson);
294 long n = 0;
295 double r = 1.0d;
296 double rnd = 1.0d;
297
298 while (n < 1000 * meanPoisson) {
299 rnd = random.nextDouble();
300 r = r * rnd;
301 if (r >= p) {
302 n++;
303 } else {
304 return n;
305 }
306 }
307 return n;
308 } else {
309 final double lambda = FastMath.floor(meanPoisson);
310 final double lambdaFractional = meanPoisson - lambda;
311 final double logLambda = FastMath.log(lambda);
312 final double logLambdaFactorial = ArithmeticUtils.factorialLog((int) lambda);
313 final long y2 = lambdaFractional < Double.MIN_VALUE ? 0 : nextPoisson(lambdaFractional);
314 final double delta = FastMath.sqrt(lambda * FastMath.log(32 * lambda / FastMath.PI + 1));
315 final double halfDelta = delta / 2;
316 final double twolpd = 2 * lambda + delta;
317 final double a1 = FastMath.sqrt(FastMath.PI * twolpd) * FastMath.exp(1 / 8 * lambda);
318 final double a2 = (twolpd / delta) * FastMath.exp(-delta * (1 + delta) / twolpd);
319 final double aSum = a1 + a2 + 1;
320 final double p1 = a1 / aSum;
321 final double p2 = a2 / aSum;
322 final double c1 = 1 / (8 * lambda);
323
324 double x = 0;
325 double y = 0;
326 double v = 0;
327 int a = 0;
328 double t = 0;
329 double qr = 0;
330 double qa = 0;
331 for (;;) {
332 final double u = random.nextDouble();
333 if (u <= p1) {
334 final double n = random.nextGaussian();
335 x = n * FastMath.sqrt(lambda + halfDelta) - 0.5d;
336 if (x > delta || x < -lambda) {
337 continue;
338 }
339 y = x < 0 ? FastMath.floor(x) : FastMath.ceil(x);
340 final double e = exponential.sample();
341 v = -e - (n * n / 2) + c1;
342 } else {
343 if (u > p1 + p2) {
344 y = lambda;
345 break;
346 } else {
347 x = delta + (twolpd / delta) * exponential.sample();
348 y = FastMath.ceil(x);
349 v = -exponential.sample() - delta * (x + 1) / twolpd;
350 }
351 }
352 a = x < 0 ? 1 : 0;
353 t = y * (y + 1) / (2 * lambda);
354 if (v < -t && a == 0) {
355 y = lambda + y;
356 break;
357 }
358 qr = t * ((2 * y + 1) / (6 * lambda) - 1);
359 qa = qr - (t * t) / (3 * (lambda + a * (y + 1)));
360 if (v < qa) {
361 y = lambda + y;
362 break;
363 }
364 if (v > qr) {
365 continue;
366 }
367 if (v < y * logLambda - ArithmeticUtils.factorialLog((int) (y + lambda)) + logLambdaFactorial) {
368 y = lambda + y;
369 break;
370 }
371 }
372 return y2 + (long) y;
373 }
374 }
375 }