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.random.RandomGenerator;
22  import org.apache.commons.math3.random.Well19937c;
23  import org.apache.commons.math3.special.Gamma;
24  import org.apache.commons.math3.util.CombinatoricsUtils;
25  import org.apache.commons.math3.util.FastMath;
26  import org.apache.commons.math3.util.MathUtils;
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   */
34  public class PoissonDistribution extends AbstractIntegerDistribution {
35      /**
36       * Default maximum number of iterations for cumulative probability calculations.
37       * @since 2.1
38       */
39      public static final int DEFAULT_MAX_ITERATIONS = 10000000;
40      /**
41       * Default convergence criterion.
42       * @since 2.1
43       */
44      public static final double DEFAULT_EPSILON = 1e-12;
45      /** Serializable version identifier. */
46      private static final long serialVersionUID = -3349935121172596109L;
47      /** Distribution used to compute normal approximation. */
48      private final NormalDistribution normal;
49      /** Distribution needed for the {@link #sample()} method. */
50      private final ExponentialDistribution exponential;
51      /** Mean of the distribution. */
52      private final double mean;
53  
54      /**
55       * Maximum number of iterations for cumulative probability. Cumulative
56       * probabilities are estimated using either Lanczos series approximation
57       * of {@link Gamma#regularizedGammaP(double, double, double, int)}
58       * or continued fraction approximation of
59       * {@link Gamma#regularizedGammaQ(double, double, double, int)}.
60       */
61      private final int maxIterations;
62  
63      /** Convergence criterion for cumulative probability. */
64      private final double epsilon;
65  
66      /**
67       * Creates a new Poisson distribution with specified mean.
68       * <p>
69       * <b>Note:</b> this constructor will implicitly create an instance of
70       * {@link Well19937c} as random generator to be used for sampling only (see
71       * {@link #sample()} and {@link #sample(int)}). In case no sampling is
72       * needed for the created distribution, it is advised to pass {@code null}
73       * as random generator via the appropriate constructors to avoid the
74       * additional initialisation overhead.
75       *
76       * @param p the Poisson mean
77       * @throws NotStrictlyPositiveException if {@code p <= 0}.
78       */
79      public PoissonDistribution(double p) throws NotStrictlyPositiveException {
80          this(p, DEFAULT_EPSILON, DEFAULT_MAX_ITERATIONS);
81      }
82  
83      /**
84       * Creates a new Poisson distribution with specified mean, convergence
85       * criterion and maximum number of iterations.
86       * <p>
87       * <b>Note:</b> this constructor will implicitly create an instance of
88       * {@link Well19937c} as random generator to be used for sampling only (see
89       * {@link #sample()} and {@link #sample(int)}). In case no sampling is
90       * needed for the created distribution, it is advised to pass {@code null}
91       * as random generator via the appropriate constructors to avoid the
92       * additional initialisation overhead.
93       *
94       * @param p Poisson mean.
95       * @param epsilon Convergence criterion for cumulative probabilities.
96       * @param maxIterations the maximum number of iterations for cumulative
97       * probabilities.
98       * @throws NotStrictlyPositiveException if {@code p <= 0}.
99       * @since 2.1
100      */
101     public PoissonDistribution(double p, double epsilon, int maxIterations)
102     throws NotStrictlyPositiveException {
103         this(new Well19937c(), p, epsilon, maxIterations);
104     }
105 
106     /**
107      * Creates a new Poisson distribution with specified mean, convergence
108      * criterion and maximum number of iterations.
109      *
110      * @param rng Random number generator.
111      * @param p Poisson mean.
112      * @param epsilon Convergence criterion for cumulative probabilities.
113      * @param maxIterations the maximum number of iterations for cumulative
114      * probabilities.
115      * @throws NotStrictlyPositiveException if {@code p <= 0}.
116      * @since 3.1
117      */
118     public PoissonDistribution(RandomGenerator rng,
119                                double p,
120                                double epsilon,
121                                int maxIterations)
122     throws NotStrictlyPositiveException {
123         super(rng);
124 
125         if (p <= 0) {
126             throw new NotStrictlyPositiveException(LocalizedFormats.MEAN, p);
127         }
128         mean = p;
129         this.epsilon = epsilon;
130         this.maxIterations = maxIterations;
131 
132         // Use the same RNG instance as the parent class.
133         normal = new NormalDistribution(rng, p, FastMath.sqrt(p),
134                                         NormalDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
135         exponential = new ExponentialDistribution(rng, 1,
136                                                   ExponentialDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
137     }
138 
139     /**
140      * Creates a new Poisson distribution with the specified mean and
141      * convergence criterion.
142      *
143      * @param p Poisson mean.
144      * @param epsilon Convergence criterion for cumulative probabilities.
145      * @throws NotStrictlyPositiveException if {@code p <= 0}.
146      * @since 2.1
147      */
148     public PoissonDistribution(double p, double epsilon)
149     throws NotStrictlyPositiveException {
150         this(p, epsilon, DEFAULT_MAX_ITERATIONS);
151     }
152 
153     /**
154      * Creates a new Poisson distribution with the specified mean and maximum
155      * number of iterations.
156      *
157      * @param p Poisson mean.
158      * @param maxIterations Maximum number of iterations for cumulative
159      * probabilities.
160      * @since 2.1
161      */
162     public PoissonDistribution(double p, int maxIterations) {
163         this(p, DEFAULT_EPSILON, maxIterations);
164     }
165 
166     /**
167      * Get the mean for the distribution.
168      *
169      * @return the mean for the distribution.
170      */
171     public double getMean() {
172         return mean;
173     }
174 
175     /** {@inheritDoc} */
176     public double probability(int x) {
177         final double logProbability = logProbability(x);
178         return logProbability == Double.NEGATIVE_INFINITY ? 0 : FastMath.exp(logProbability);
179     }
180 
181     /** {@inheritDoc} */
182     @Override
183     public double logProbability(int x) {
184         double ret;
185         if (x < 0 || x == Integer.MAX_VALUE) {
186             ret = Double.NEGATIVE_INFINITY;
187         } else if (x == 0) {
188             ret = -mean;
189         } else {
190             ret = -SaddlePointExpansion.getStirlingError(x) -
191                   SaddlePointExpansion.getDeviancePart(x, mean) -
192                   0.5 * FastMath.log(MathUtils.TWO_PI) - 0.5 * FastMath.log(x);
193         }
194         return ret;
195     }
196 
197     /** {@inheritDoc} */
198     public double cumulativeProbability(int x) {
199         if (x < 0) {
200             return 0;
201         }
202         if (x == Integer.MAX_VALUE) {
203             return 1;
204         }
205         return Gamma.regularizedGammaQ((double) x + 1, mean, epsilon,
206                                        maxIterations);
207     }
208 
209     /**
210      * Calculates the Poisson distribution function using a normal
211      * approximation. The {@code N(mean, sqrt(mean))} distribution is used
212      * to approximate the Poisson distribution. The computation uses
213      * "half-correction" (evaluating the normal distribution function at
214      * {@code x + 0.5}).
215      *
216      * @param x Upper bound, inclusive.
217      * @return the distribution function value calculated using a normal
218      * approximation.
219      */
220     public double normalApproximateProbability(int x)  {
221         // calculate the probability using half-correction
222         return normal.cumulativeProbability(x + 0.5);
223     }
224 
225     /**
226      * {@inheritDoc}
227      *
228      * For mean parameter {@code p}, the mean is {@code p}.
229      */
230     public double getNumericalMean() {
231         return getMean();
232     }
233 
234     /**
235      * {@inheritDoc}
236      *
237      * For mean parameter {@code p}, the variance is {@code p}.
238      */
239     public double getNumericalVariance() {
240         return getMean();
241     }
242 
243     /**
244      * {@inheritDoc}
245      *
246      * The lower bound of the support is always 0 no matter the mean parameter.
247      *
248      * @return lower bound of the support (always 0)
249      */
250     public int getSupportLowerBound() {
251         return 0;
252     }
253 
254     /**
255      * {@inheritDoc}
256      *
257      * The upper bound of the support is positive infinity,
258      * regardless of the parameter values. There is no integer infinity,
259      * so this method returns {@code Integer.MAX_VALUE}.
260      *
261      * @return upper bound of the support (always {@code Integer.MAX_VALUE} for
262      * positive infinity)
263      */
264     public int getSupportUpperBound() {
265         return Integer.MAX_VALUE;
266     }
267 
268     /**
269      * {@inheritDoc}
270      *
271      * The support of this distribution is connected.
272      *
273      * @return {@code true}
274      */
275     public boolean isSupportConnected() {
276         return true;
277     }
278 
279     /**
280      * {@inheritDoc}
281      * <p>
282      * <strong>Algorithm Description</strong>:
283      * <ul>
284      *  <li>For small means, uses simulation of a Poisson process
285      *   using Uniform deviates, as described
286      *   <a href="http://irmi.epfl.ch/cmos/Pmmi/interactive/rng7.htm"> here</a>.
287      *   The Poisson process (and hence value returned) is bounded by 1000 * mean.
288      *  </li>
289      *  <li>For large means, uses the rejection algorithm described in
290      *   <quote>
291      *    Devroye, Luc. (1981).<i>The Computer Generation of Poisson Random Variables</i>
292      *    <strong>Computing</strong> vol. 26 pp. 197-207.
293      *   </quote>
294      *  </li>
295      * </ul>
296      * </p>
297      *
298      * @return a random value.
299      * @since 2.2
300      */
301     @Override
302     public int sample() {
303         return (int) FastMath.min(nextPoisson(mean), Integer.MAX_VALUE);
304     }
305 
306     /**
307      * @param meanPoisson Mean of the Poisson distribution.
308      * @return the next sample.
309      */
310     private long nextPoisson(double meanPoisson) {
311         final double pivot = 40.0d;
312         if (meanPoisson < pivot) {
313             double p = FastMath.exp(-meanPoisson);
314             long n = 0;
315             double r = 1.0d;
316             double rnd = 1.0d;
317 
318             while (n < 1000 * meanPoisson) {
319                 rnd = random.nextDouble();
320                 r *= rnd;
321                 if (r >= p) {
322                     n++;
323                 } else {
324                     return n;
325                 }
326             }
327             return n;
328         } else {
329             final double lambda = FastMath.floor(meanPoisson);
330             final double lambdaFractional = meanPoisson - lambda;
331             final double logLambda = FastMath.log(lambda);
332             final double logLambdaFactorial = CombinatoricsUtils.factorialLog((int) lambda);
333             final long y2 = lambdaFractional < Double.MIN_VALUE ? 0 : nextPoisson(lambdaFractional);
334             final double delta = FastMath.sqrt(lambda * FastMath.log(32 * lambda / FastMath.PI + 1));
335             final double halfDelta = delta / 2;
336             final double twolpd = 2 * lambda + delta;
337             final double a1 = FastMath.sqrt(FastMath.PI * twolpd) * FastMath.exp(1 / (8 * lambda));
338             final double a2 = (twolpd / delta) * FastMath.exp(-delta * (1 + delta) / twolpd);
339             final double aSum = a1 + a2 + 1;
340             final double p1 = a1 / aSum;
341             final double p2 = a2 / aSum;
342             final double c1 = 1 / (8 * lambda);
343 
344             double x = 0;
345             double y = 0;
346             double v = 0;
347             int a = 0;
348             double t = 0;
349             double qr = 0;
350             double qa = 0;
351             for (;;) {
352                 final double u = random.nextDouble();
353                 if (u <= p1) {
354                     final double n = random.nextGaussian();
355                     x = n * FastMath.sqrt(lambda + halfDelta) - 0.5d;
356                     if (x > delta || x < -lambda) {
357                         continue;
358                     }
359                     y = x < 0 ? FastMath.floor(x) : FastMath.ceil(x);
360                     final double e = exponential.sample();
361                     v = -e - (n * n / 2) + c1;
362                 } else {
363                     if (u > p1 + p2) {
364                         y = lambda;
365                         break;
366                     } else {
367                         x = delta + (twolpd / delta) * exponential.sample();
368                         y = FastMath.ceil(x);
369                         v = -exponential.sample() - delta * (x + 1) / twolpd;
370                     }
371                 }
372                 a = x < 0 ? 1 : 0;
373                 t = y * (y + 1) / (2 * lambda);
374                 if (v < -t && a == 0) {
375                     y = lambda + y;
376                     break;
377                 }
378                 qr = t * ((2 * y + 1) / (6 * lambda) - 1);
379                 qa = qr - (t * t) / (3 * (lambda + a * (y + 1)));
380                 if (v < qa) {
381                     y = lambda + y;
382                     break;
383                 }
384                 if (v > qr) {
385                     continue;
386                 }
387                 if (v < y * logLambda - CombinatoricsUtils.factorialLog((int) (y + lambda)) + logLambdaFactorial) {
388                     y = lambda + y;
389                     break;
390                 }
391             }
392             return y2 + (long) y;
393         }
394     }
395 }