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 */
017package org.apache.commons.math3.distribution;
018
019import org.apache.commons.math3.exception.NotStrictlyPositiveException;
020import org.apache.commons.math3.exception.util.LocalizedFormats;
021import org.apache.commons.math3.special.Gamma;
022import org.apache.commons.math3.util.CombinatoricsUtils;
023import org.apache.commons.math3.util.MathUtils;
024import org.apache.commons.math3.util.FastMath;
025import org.apache.commons.math3.random.RandomGenerator;
026import 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 1540217 2013-11-08 23:27:49Z psteitz $
034 */
035public 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        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}