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    *      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;
18  
19  import org.apache.commons.math3.exception.NotStrictlyPositiveException;
20  import org.apache.commons.math3.exception.util.LocalizedFormats;
21  import org.apache.commons.math3.special.Gamma;
22  import org.apache.commons.math3.util.CombinatoricsUtils;
23  import org.apache.commons.math3.util.MathUtils;
24  import org.apache.commons.math3.util.FastMath;
25  import org.apache.commons.math3.random.RandomGenerator;
26  import org.apache.commons.math3.random.Well19937c;
27  
28  /**
29   * Implementation of the Poisson distribution.
30   *
31   * @see <a href="http://en.wikipedia.org/wiki/Poisson_distribution">Poisson distribution (Wikipedia)</a>
32   * @see <a href="http://mathworld.wolfram.com/PoissonDistribution.html">Poisson distribution (MathWorld)</a>
33   * @version $Id: PoissonDistribution.java 1540217 2013-11-08 23:27:49Z psteitz $
34   */
35  public class PoissonDistribution extends AbstractIntegerDistribution {
36      /**
37       * Default maximum number of iterations for cumulative probability calculations.
38       * @since 2.1
39       */
40      public static final int DEFAULT_MAX_ITERATIONS = 10000000;
41      /**
42       * Default convergence criterion.
43       * @since 2.1
44       */
45      public static final double DEFAULT_EPSILON = 1e-12;
46      /** Serializable version identifier. */
47      private static final long serialVersionUID = -3349935121172596109L;
48      /** Distribution used to compute normal approximation. */
49      private final NormalDistribution normal;
50      /** Distribution needed for the {@link #sample()} method. */
51      private final ExponentialDistribution exponential;
52      /** Mean of the distribution. */
53      private final double mean;
54  
55      /**
56       * Maximum number of iterations for cumulative probability. Cumulative
57       * probabilities are estimated using either Lanczos series approximation
58       * of {@link Gamma#regularizedGammaP(double, double, double, int)}
59       * or continued fraction approximation of
60       * {@link Gamma#regularizedGammaQ(double, double, double, int)}.
61       */
62      private final int maxIterations;
63  
64      /** Convergence criterion for cumulative probability. */
65      private final double epsilon;
66  
67      /**
68       * Creates a new Poisson distribution with specified mean.
69       *
70       * @param p the Poisson mean
71       * @throws NotStrictlyPositiveException if {@code p <= 0}.
72       */
73      public PoissonDistribution(double p) throws NotStrictlyPositiveException {
74          this(p, DEFAULT_EPSILON, DEFAULT_MAX_ITERATIONS);
75      }
76  
77      /**
78       * Creates a new Poisson distribution with specified mean, convergence
79       * criterion and maximum number of iterations.
80       *
81       * @param p Poisson mean.
82       * @param epsilon Convergence criterion for cumulative probabilities.
83       * @param maxIterations the maximum number of iterations for cumulative
84       * probabilities.
85       * @throws NotStrictlyPositiveException if {@code p <= 0}.
86       * @since 2.1
87       */
88      public PoissonDistribution(double p, double epsilon, int maxIterations)
89      throws NotStrictlyPositiveException {
90          this(new Well19937c(), p, epsilon, maxIterations);
91      }
92  
93      /**
94       * Creates a new Poisson distribution with specified mean, convergence
95       * criterion and maximum number of iterations.
96       *
97       * @param rng Random number generator.
98       * @param p Poisson mean.
99       * @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         final double logProbability = logProbability(x);
165         return logProbability == Double.NEGATIVE_INFINITY ? 0 : FastMath.exp(logProbability);
166     }
167 
168     /** {@inheritDoc} */
169     @Override
170     public double logProbability(int x) {
171         double ret;
172         if (x < 0 || x == Integer.MAX_VALUE) {
173             ret = Double.NEGATIVE_INFINITY;
174         } else if (x == 0) {
175             ret = -mean;
176         } else {
177             ret = -SaddlePointExpansion.getStirlingError(x) -
178                   SaddlePointExpansion.getDeviancePart(x, mean) -
179                   0.5 * FastMath.log(MathUtils.TWO_PI) - 0.5 * FastMath.log(x);
180         }
181         return ret;
182     }
183 
184     /** {@inheritDoc} */
185     public double cumulativeProbability(int x) {
186         if (x < 0) {
187             return 0;
188         }
189         if (x == Integer.MAX_VALUE) {
190             return 1;
191         }
192         return Gamma.regularizedGammaQ((double) x + 1, mean, epsilon,
193                                        maxIterations);
194     }
195 
196     /**
197      * Calculates the Poisson distribution function using a normal
198      * approximation. The {@code N(mean, sqrt(mean))} distribution is used
199      * to approximate the Poisson distribution. The computation uses
200      * "half-correction" (evaluating the normal distribution function at
201      * {@code x + 0.5}).
202      *
203      * @param x Upper bound, inclusive.
204      * @return the distribution function value calculated using a normal
205      * approximation.
206      */
207     public double normalApproximateProbability(int x)  {
208         // calculate the probability using half-correction
209         return normal.cumulativeProbability(x + 0.5);
210     }
211 
212     /**
213      * {@inheritDoc}
214      *
215      * For mean parameter {@code p}, the mean is {@code p}.
216      */
217     public double getNumericalMean() {
218         return getMean();
219     }
220 
221     /**
222      * {@inheritDoc}
223      *
224      * For mean parameter {@code p}, the variance is {@code p}.
225      */
226     public double getNumericalVariance() {
227         return getMean();
228     }
229 
230     /**
231      * {@inheritDoc}
232      *
233      * The lower bound of the support is always 0 no matter the mean parameter.
234      *
235      * @return lower bound of the support (always 0)
236      */
237     public int getSupportLowerBound() {
238         return 0;
239     }
240 
241     /**
242      * {@inheritDoc}
243      *
244      * The upper bound of the support is positive infinity,
245      * regardless of the parameter values. There is no integer infinity,
246      * so this method returns {@code Integer.MAX_VALUE}.
247      *
248      * @return upper bound of the support (always {@code Integer.MAX_VALUE} for
249      * positive infinity)
250      */
251     public int getSupportUpperBound() {
252         return Integer.MAX_VALUE;
253     }
254 
255     /**
256      * {@inheritDoc}
257      *
258      * The support of this distribution is connected.
259      *
260      * @return {@code true}
261      */
262     public boolean isSupportConnected() {
263         return true;
264     }
265 
266     /**
267      * {@inheritDoc}
268      * <p>
269      * <strong>Algorithm Description</strong>:
270      * <ul>
271      *  <li>For small means, uses simulation of a Poisson process
272      *   using Uniform deviates, as described
273      *   <a href="http://irmi.epfl.ch/cmos/Pmmi/interactive/rng7.htm"> here</a>.
274      *   The Poisson process (and hence value returned) is bounded by 1000 * mean.
275      *  </li>
276      *  <li>For large means, uses the rejection algorithm described in
277      *   <quote>
278      *    Devroye, Luc. (1981).<i>The Computer Generation of Poisson Random Variables</i>
279      *    <strong>Computing</strong> vol. 26 pp. 197-207.
280      *   </quote>
281      *  </li>
282      * </ul>
283      * </p>
284      *
285      * @return a random value.
286      * @since 2.2
287      */
288     @Override
289     public int sample() {
290         return (int) FastMath.min(nextPoisson(mean), Integer.MAX_VALUE);
291     }
292 
293     /**
294      * @param meanPoisson Mean of the Poisson distribution.
295      * @return the next sample.
296      */
297     private long nextPoisson(double meanPoisson) {
298         final double pivot = 40.0d;
299         if (meanPoisson < pivot) {
300             double p = FastMath.exp(-meanPoisson);
301             long n = 0;
302             double r = 1.0d;
303             double rnd = 1.0d;
304 
305             while (n < 1000 * meanPoisson) {
306                 rnd = random.nextDouble();
307                 r *= rnd;
308                 if (r >= p) {
309                     n++;
310                 } else {
311                     return n;
312                 }
313             }
314             return n;
315         } else {
316             final double lambda = FastMath.floor(meanPoisson);
317             final double lambdaFractional = meanPoisson - lambda;
318             final double logLambda = FastMath.log(lambda);
319             final double logLambdaFactorial = CombinatoricsUtils.factorialLog((int) lambda);
320             final long y2 = lambdaFractional < Double.MIN_VALUE ? 0 : nextPoisson(lambdaFractional);
321             final double delta = FastMath.sqrt(lambda * FastMath.log(32 * lambda / FastMath.PI + 1));
322             final double halfDelta = delta / 2;
323             final double twolpd = 2 * lambda + delta;
324             final double a1 = FastMath.sqrt(FastMath.PI * twolpd) * FastMath.exp(1 / (8 * lambda));
325             final double a2 = (twolpd / delta) * FastMath.exp(-delta * (1 + delta) / twolpd);
326             final double aSum = a1 + a2 + 1;
327             final double p1 = a1 / aSum;
328             final double p2 = a2 / aSum;
329             final double c1 = 1 / (8 * lambda);
330 
331             double x = 0;
332             double y = 0;
333             double v = 0;
334             int a = 0;
335             double t = 0;
336             double qr = 0;
337             double qa = 0;
338             for (;;) {
339                 final double u = random.nextDouble();
340                 if (u <= p1) {
341                     final double n = random.nextGaussian();
342                     x = n * FastMath.sqrt(lambda + halfDelta) - 0.5d;
343                     if (x > delta || x < -lambda) {
344                         continue;
345                     }
346                     y = x < 0 ? FastMath.floor(x) : FastMath.ceil(x);
347                     final double e = exponential.sample();
348                     v = -e - (n * n / 2) + c1;
349                 } else {
350                     if (u > p1 + p2) {
351                         y = lambda;
352                         break;
353                     } else {
354                         x = delta + (twolpd / delta) * exponential.sample();
355                         y = FastMath.ceil(x);
356                         v = -exponential.sample() - delta * (x + 1) / twolpd;
357                     }
358                 }
359                 a = x < 0 ? 1 : 0;
360                 t = y * (y + 1) / (2 * lambda);
361                 if (v < -t && a == 0) {
362                     y = lambda + y;
363                     break;
364                 }
365                 qr = t * ((2 * y + 1) / (6 * lambda) - 1);
366                 qa = qr - (t * t) / (3 * (lambda + a * (y + 1)));
367                 if (v < qa) {
368                     y = lambda + y;
369                     break;
370                 }
371                 if (v > qr) {
372                     continue;
373                 }
374                 if (v < y * logLambda - CombinatoricsUtils.factorialLog((int) (y + lambda)) + logLambdaFactorial) {
375                     y = lambda + y;
376                     break;
377                 }
378             }
379             return y2 + (long) y;
380         }
381     }
382 }