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    }