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.random.RandomGenerator;
022import org.apache.commons.math3.random.Well19937c;
023import org.apache.commons.math3.special.Gamma;
024import org.apache.commons.math3.util.CombinatoricsUtils;
025import org.apache.commons.math3.util.FastMath;
026import org.apache.commons.math3.util.MathUtils;
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 */
034public class PoissonDistribution extends AbstractIntegerDistribution {
035    /**
036     * Default maximum number of iterations for cumulative probability calculations.
037     * @since 2.1
038     */
039    public static final int DEFAULT_MAX_ITERATIONS = 10000000;
040    /**
041     * Default convergence criterion.
042     * @since 2.1
043     */
044    public static final double DEFAULT_EPSILON = 1e-12;
045    /** Serializable version identifier. */
046    private static final long serialVersionUID = -3349935121172596109L;
047    /** Distribution used to compute normal approximation. */
048    private final NormalDistribution normal;
049    /** Distribution needed for the {@link #sample()} method. */
050    private final ExponentialDistribution exponential;
051    /** Mean of the distribution. */
052    private final double mean;
053
054    /**
055     * Maximum number of iterations for cumulative probability. Cumulative
056     * probabilities are estimated using either Lanczos series approximation
057     * of {@link Gamma#regularizedGammaP(double, double, double, int)}
058     * or continued fraction approximation of
059     * {@link Gamma#regularizedGammaQ(double, double, double, int)}.
060     */
061    private final int maxIterations;
062
063    /** Convergence criterion for cumulative probability. */
064    private final double epsilon;
065
066    /**
067     * Creates a new Poisson distribution with specified mean.
068     * <p>
069     * <b>Note:</b> this constructor will implicitly create an instance of
070     * {@link Well19937c} as random generator to be used for sampling only (see
071     * {@link #sample()} and {@link #sample(int)}). In case no sampling is
072     * needed for the created distribution, it is advised to pass {@code null}
073     * as random generator via the appropriate constructors to avoid the
074     * additional initialisation overhead.
075     *
076     * @param p the Poisson mean
077     * @throws NotStrictlyPositiveException if {@code p <= 0}.
078     */
079    public PoissonDistribution(double p) throws NotStrictlyPositiveException {
080        this(p, DEFAULT_EPSILON, DEFAULT_MAX_ITERATIONS);
081    }
082
083    /**
084     * Creates a new Poisson distribution with specified mean, convergence
085     * criterion and maximum number of iterations.
086     * <p>
087     * <b>Note:</b> this constructor will implicitly create an instance of
088     * {@link Well19937c} as random generator to be used for sampling only (see
089     * {@link #sample()} and {@link #sample(int)}). In case no sampling is
090     * needed for the created distribution, it is advised to pass {@code null}
091     * as random generator via the appropriate constructors to avoid the
092     * additional initialisation overhead.
093     *
094     * @param p Poisson mean.
095     * @param epsilon Convergence criterion for cumulative probabilities.
096     * @param maxIterations the maximum number of iterations for cumulative
097     * probabilities.
098     * @throws NotStrictlyPositiveException if {@code p <= 0}.
099     * @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}