    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    }