PoissonDistribution.java

  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. import org.apache.commons.math3.exception.NotStrictlyPositiveException;
  19. import org.apache.commons.math3.exception.util.LocalizedFormats;
  20. import org.apache.commons.math3.random.RandomGenerator;
  21. import org.apache.commons.math3.random.Well19937c;
  22. import org.apache.commons.math3.special.Gamma;
  23. import org.apache.commons.math3.util.CombinatoricsUtils;
  24. import org.apache.commons.math3.util.FastMath;
  25. import org.apache.commons.math3.util.MathUtils;

  26. /**
  27.  * Implementation of the Poisson distribution.
  28.  *
  29.  * @see <a href="http://en.wikipedia.org/wiki/Poisson_distribution">Poisson distribution (Wikipedia)</a>
  30.  * @see <a href="http://mathworld.wolfram.com/PoissonDistribution.html">Poisson distribution (MathWorld)</a>
  31.  */
  32. public class PoissonDistribution extends AbstractIntegerDistribution {
  33.     /**
  34.      * Default maximum number of iterations for cumulative probability calculations.
  35.      * @since 2.1
  36.      */
  37.     public static final int DEFAULT_MAX_ITERATIONS = 10000000;
  38.     /**
  39.      * Default convergence criterion.
  40.      * @since 2.1
  41.      */
  42.     public static final double DEFAULT_EPSILON = 1e-12;
  43.     /** Serializable version identifier. */
  44.     private static final long serialVersionUID = -3349935121172596109L;
  45.     /** Distribution used to compute normal approximation. */
  46.     private final NormalDistribution normal;
  47.     /** Distribution needed for the {@link #sample()} method. */
  48.     private final ExponentialDistribution exponential;
  49.     /** Mean of the distribution. */
  50.     private final double mean;

  51.     /**
  52.      * Maximum number of iterations for cumulative probability. Cumulative
  53.      * probabilities are estimated using either Lanczos series approximation
  54.      * of {@link Gamma#regularizedGammaP(double, double, double, int)}
  55.      * or continued fraction approximation of
  56.      * {@link Gamma#regularizedGammaQ(double, double, double, int)}.
  57.      */
  58.     private final int maxIterations;

  59.     /** Convergence criterion for cumulative probability. */
  60.     private final double epsilon;

  61.     /**
  62.      * Creates a new Poisson distribution with specified mean.
  63.      * <p>
  64.      * <b>Note:</b> this constructor will implicitly create an instance of
  65.      * {@link Well19937c} as random generator to be used for sampling only (see
  66.      * {@link #sample()} and {@link #sample(int)}). In case no sampling is
  67.      * needed for the created distribution, it is advised to pass {@code null}
  68.      * as random generator via the appropriate constructors to avoid the
  69.      * additional initialisation overhead.
  70.      *
  71.      * @param p the Poisson mean
  72.      * @throws NotStrictlyPositiveException if {@code p <= 0}.
  73.      */
  74.     public PoissonDistribution(double p) throws NotStrictlyPositiveException {
  75.         this(p, DEFAULT_EPSILON, DEFAULT_MAX_ITERATIONS);
  76.     }

  77.     /**
  78.      * Creates a new Poisson distribution with specified mean, convergence
  79.      * criterion and maximum number of iterations.
  80.      * <p>
  81.      * <b>Note:</b> this constructor will implicitly create an instance of
  82.      * {@link Well19937c} as random generator to be used for sampling only (see
  83.      * {@link #sample()} and {@link #sample(int)}). In case no sampling is
  84.      * needed for the created distribution, it is advised to pass {@code null}
  85.      * as random generator via the appropriate constructors to avoid the
  86.      * additional initialisation overhead.
  87.      *
  88.      * @param p Poisson mean.
  89.      * @param epsilon Convergence criterion for cumulative probabilities.
  90.      * @param maxIterations the maximum number of iterations for cumulative
  91.      * probabilities.
  92.      * @throws NotStrictlyPositiveException if {@code p <= 0}.
  93.      * @since 2.1
  94.      */
  95.     public PoissonDistribution(double p, double epsilon, int maxIterations)
  96.     throws NotStrictlyPositiveException {
  97.         this(new Well19937c(), p, epsilon, maxIterations);
  98.     }

  99.     /**
  100.      * Creates a new Poisson distribution with specified mean, convergence
  101.      * criterion and maximum number of iterations.
  102.      *
  103.      * @param rng Random number generator.
  104.      * @param p Poisson mean.
  105.      * @param epsilon Convergence criterion for cumulative probabilities.
  106.      * @param maxIterations the maximum number of iterations for cumulative
  107.      * probabilities.
  108.      * @throws NotStrictlyPositiveException if {@code p <= 0}.
  109.      * @since 3.1
  110.      */
  111.     public PoissonDistribution(RandomGenerator rng,
  112.                                double p,
  113.                                double epsilon,
  114.                                int maxIterations)
  115.     throws NotStrictlyPositiveException {
  116.         super(rng);

  117.         if (p <= 0) {
  118.             throw new NotStrictlyPositiveException(LocalizedFormats.MEAN, p);
  119.         }
  120.         mean = p;
  121.         this.epsilon = epsilon;
  122.         this.maxIterations = maxIterations;

  123.         // Use the same RNG instance as the parent class.
  124.         normal = new NormalDistribution(rng, p, FastMath.sqrt(p),
  125.                                         NormalDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
  126.         exponential = new ExponentialDistribution(rng, 1,
  127.                                                   ExponentialDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
  128.     }

  129.     /**
  130.      * Creates a new Poisson distribution with the specified mean and
  131.      * convergence criterion.
  132.      *
  133.      * @param p Poisson mean.
  134.      * @param epsilon Convergence criterion for cumulative probabilities.
  135.      * @throws NotStrictlyPositiveException if {@code p <= 0}.
  136.      * @since 2.1
  137.      */
  138.     public PoissonDistribution(double p, double epsilon)
  139.     throws NotStrictlyPositiveException {
  140.         this(p, epsilon, DEFAULT_MAX_ITERATIONS);
  141.     }

  142.     /**
  143.      * Creates a new Poisson distribution with the specified mean and maximum
  144.      * number of iterations.
  145.      *
  146.      * @param p Poisson mean.
  147.      * @param maxIterations Maximum number of iterations for cumulative
  148.      * probabilities.
  149.      * @since 2.1
  150.      */
  151.     public PoissonDistribution(double p, int maxIterations) {
  152.         this(p, DEFAULT_EPSILON, maxIterations);
  153.     }

  154.     /**
  155.      * Get the mean for the distribution.
  156.      *
  157.      * @return the mean for the distribution.
  158.      */
  159.     public double getMean() {
  160.         return mean;
  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.     /** {@inheritDoc} */
  168.     @Override
  169.     public double logProbability(int x) {
  170.         double ret;
  171.         if (x < 0 || x == Integer.MAX_VALUE) {
  172.             ret = Double.NEGATIVE_INFINITY;
  173.         } else if (x == 0) {
  174.             ret = -mean;
  175.         } else {
  176.             ret = -SaddlePointExpansion.getStirlingError(x) -
  177.                   SaddlePointExpansion.getDeviancePart(x, mean) -
  178.                   0.5 * FastMath.log(MathUtils.TWO_PI) - 0.5 * FastMath.log(x);
  179.         }
  180.         return ret;
  181.     }

  182.     /** {@inheritDoc} */
  183.     public double cumulativeProbability(int x) {
  184.         if (x < 0) {
  185.             return 0;
  186.         }
  187.         if (x == Integer.MAX_VALUE) {
  188.             return 1;
  189.         }
  190.         return Gamma.regularizedGammaQ((double) x + 1, mean, epsilon,
  191.                                        maxIterations);
  192.     }

  193.     /**
  194.      * Calculates the Poisson distribution function using a normal
  195.      * approximation. The {@code N(mean, sqrt(mean))} distribution is used
  196.      * to approximate the Poisson distribution. The computation uses
  197.      * "half-correction" (evaluating the normal distribution function at
  198.      * {@code x + 0.5}).
  199.      *
  200.      * @param x Upper bound, inclusive.
  201.      * @return the distribution function value calculated using a normal
  202.      * approximation.
  203.      */
  204.     public double normalApproximateProbability(int x)  {
  205.         // calculate the probability using half-correction
  206.         return normal.cumulativeProbability(x + 0.5);
  207.     }

  208.     /**
  209.      * {@inheritDoc}
  210.      *
  211.      * For mean parameter {@code p}, the mean is {@code p}.
  212.      */
  213.     public double getNumericalMean() {
  214.         return getMean();
  215.     }

  216.     /**
  217.      * {@inheritDoc}
  218.      *
  219.      * For mean parameter {@code p}, the variance is {@code p}.
  220.      */
  221.     public double getNumericalVariance() {
  222.         return getMean();
  223.     }

  224.     /**
  225.      * {@inheritDoc}
  226.      *
  227.      * The lower bound of the support is always 0 no matter the mean parameter.
  228.      *
  229.      * @return lower bound of the support (always 0)
  230.      */
  231.     public int getSupportLowerBound() {
  232.         return 0;
  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.      * {@inheritDoc}
  249.      *
  250.      * The support of this distribution is connected.
  251.      *
  252.      * @return {@code true}
  253.      */
  254.     public boolean isSupportConnected() {
  255.         return true;
  256.     }

  257.     /**
  258.      * {@inheritDoc}
  259.      * <p>
  260.      * <strong>Algorithm Description</strong>:
  261.      * <ul>
  262.      *  <li>For small means, uses simulation of a Poisson process
  263.      *   using Uniform deviates, as described
  264.      *   <a href="http://mathaa.epfl.ch/cours/PMMI2001/interactive/rng7.htm"> here</a>.
  265.      *   The Poisson process (and hence value returned) is bounded by 1000 * mean.
  266.      *  </li>
  267.      *  <li>For large means, uses the rejection algorithm described in
  268.      *   <blockquote>
  269.      *    Devroye, Luc. (1981).<i>The Computer Generation of Poisson Random Variables</i><br>
  270.      *    <strong>Computing</strong> vol. 26 pp. 197-207.<br>
  271.      *   </blockquote>
  272.      *  </li>
  273.      * </ul>
  274.      * </p>
  275.      *
  276.      * @return a random value.
  277.      * @since 2.2
  278.      */
  279.     @Override
  280.     public int sample() {
  281.         return (int) FastMath.min(nextPoisson(mean), Integer.MAX_VALUE);
  282.     }

  283.     /**
  284.      * @param meanPoisson Mean of the Poisson distribution.
  285.      * @return the next sample.
  286.      */
  287.     private long nextPoisson(double meanPoisson) {
  288.         final double pivot = 40.0d;
  289.         if (meanPoisson < pivot) {
  290.             double p = FastMath.exp(-meanPoisson);
  291.             long n = 0;
  292.             double r = 1.0d;
  293.             double rnd = 1.0d;

  294.             while (n < 1000 * meanPoisson) {
  295.                 rnd = random.nextDouble();
  296.                 r *= rnd;
  297.                 if (r >= p) {
  298.                     n++;
  299.                 } else {
  300.                     return n;
  301.                 }
  302.             }
  303.             return n;
  304.         } else {
  305.             final double lambda = FastMath.floor(meanPoisson);
  306.             final double lambdaFractional = meanPoisson - lambda;
  307.             final double logLambda = FastMath.log(lambda);
  308.             final double logLambdaFactorial = CombinatoricsUtils.factorialLog((int) lambda);
  309.             final long y2 = lambdaFractional < Double.MIN_VALUE ? 0 : nextPoisson(lambdaFractional);
  310.             final double delta = FastMath.sqrt(lambda * FastMath.log(32 * lambda / FastMath.PI + 1));
  311.             final double halfDelta = delta / 2;
  312.             final double twolpd = 2 * lambda + delta;
  313.             final double a1 = FastMath.sqrt(FastMath.PI * twolpd) * FastMath.exp(1 / (8 * lambda));
  314.             final double a2 = (twolpd / delta) * FastMath.exp(-delta * (1 + delta) / twolpd);
  315.             final double aSum = a1 + a2 + 1;
  316.             final double p1 = a1 / aSum;
  317.             final double p2 = a2 / aSum;
  318.             final double c1 = 1 / (8 * lambda);

  319.             double x = 0;
  320.             double y = 0;
  321.             double v = 0;
  322.             int a = 0;
  323.             double t = 0;
  324.             double qr = 0;
  325.             double qa = 0;
  326.             for (;;) {
  327.                 final double u = random.nextDouble();
  328.                 if (u <= p1) {
  329.                     final double n = random.nextGaussian();
  330.                     x = n * FastMath.sqrt(lambda + halfDelta) - 0.5d;
  331.                     if (x > delta || x < -lambda) {
  332.                         continue;
  333.                     }
  334.                     y = x < 0 ? FastMath.floor(x) : FastMath.ceil(x);
  335.                     final double e = exponential.sample();
  336.                     v = -e - (n * n / 2) + c1;
  337.                 } else {
  338.                     if (u > p1 + p2) {
  339.                         y = lambda;
  340.                         break;
  341.                     } else {
  342.                         x = delta + (twolpd / delta) * exponential.sample();
  343.                         y = FastMath.ceil(x);
  344.                         v = -exponential.sample() - delta * (x + 1) / twolpd;
  345.                     }
  346.                 }
  347.                 a = x < 0 ? 1 : 0;
  348.                 t = y * (y + 1) / (2 * lambda);
  349.                 if (v < -t && a == 0) {
  350.                     y = lambda + y;
  351.                     break;
  352.                 }
  353.                 qr = t * ((2 * y + 1) / (6 * lambda) - 1);
  354.                 qa = qr - (t * t) / (3 * (lambda + a * (y + 1)));
  355.                 if (v < qa) {
  356.                     y = lambda + y;
  357.                     break;
  358.                 }
  359.                 if (v > qr) {
  360.                     continue;
  361.                 }
  362.                 if (v < y * logLambda - CombinatoricsUtils.factorialLog((int) (y + lambda)) + logLambdaFactorial) {
  363.                     y = lambda + y;
  364.                     break;
  365.                 }
  366.             }
  367.             return y2 + (long) y;
  368.         }
  369.     }
  370. }