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.statistics.distribution;

  18. import org.apache.commons.numbers.gamma.RegularizedGamma;
  19. import org.apache.commons.rng.UniformRandomProvider;
  20. import org.apache.commons.rng.sampling.distribution.GaussianSampler;
  21. import org.apache.commons.rng.sampling.distribution.PoissonSampler;
  22. import org.apache.commons.rng.sampling.distribution.SharedStateContinuousSampler;
  23. import org.apache.commons.rng.sampling.distribution.ZigguratSampler;

  24. /**
  25.  * Implementation of the Poisson distribution.
  26.  *
  27.  * <p>The probability mass function of \( X \) is:
  28.  *
  29.  * <p>\[ f(k; \lambda) = \frac{\lambda^k e^{-k}}{k!} \]
  30.  *
  31.  * <p>for \( \lambda \in (0, \infty) \) the mean and
  32.  * \( k \in \{0, 1, 2, \dots\} \) the number of events.
  33.  *
  34.  * @see <a href="https://en.wikipedia.org/wiki/Poisson_distribution">Poisson distribution (Wikipedia)</a>
  35.  * @see <a href="https://mathworld.wolfram.com/PoissonDistribution.html">Poisson distribution (MathWorld)</a>
  36.  */
  37. public final class PoissonDistribution extends AbstractDiscreteDistribution {
  38.     /** Upper bound on the mean to use the PoissonSampler. */
  39.     private static final double MAX_MEAN = 0.5 * Integer.MAX_VALUE;
  40.     /** Mean of the distribution. */
  41.     private final double mean;

  42.     /**
  43.      * @param mean Poisson mean.
  44.      * probabilities.
  45.      */
  46.     private PoissonDistribution(double mean) {
  47.         this.mean = mean;
  48.     }

  49.     /**
  50.      * Creates a Poisson distribution.
  51.      *
  52.      * @param mean Poisson mean.
  53.      * @return the distribution
  54.      * @throws IllegalArgumentException if {@code mean <= 0}.
  55.      */
  56.     public static PoissonDistribution of(double mean) {
  57.         if (mean <= 0) {
  58.             throw new DistributionException(DistributionException.NOT_STRICTLY_POSITIVE, mean);
  59.         }
  60.         return new PoissonDistribution(mean);
  61.     }

  62.     /** {@inheritDoc} */
  63.     @Override
  64.     public double probability(int x) {
  65.         return Math.exp(logProbability(x));
  66.     }

  67.     /** {@inheritDoc} */
  68.     @Override
  69.     public double logProbability(int x) {
  70.         if (x < 0) {
  71.             return Double.NEGATIVE_INFINITY;
  72.         } else if (x == 0) {
  73.             return -mean;
  74.         }
  75.         return -SaddlePointExpansionUtils.getStirlingError(x) -
  76.               SaddlePointExpansionUtils.getDeviancePart(x, mean) -
  77.               Constants.HALF_LOG_TWO_PI - 0.5 * Math.log(x);
  78.     }

  79.     /** {@inheritDoc} */
  80.     @Override
  81.     public double cumulativeProbability(int x) {
  82.         if (x < 0) {
  83.             return 0;
  84.         } else if (x == 0) {
  85.             return Math.exp(-mean);
  86.         }
  87.         return RegularizedGamma.Q.value((double) x + 1, mean);
  88.     }

  89.     /** {@inheritDoc} */
  90.     @Override
  91.     public double survivalProbability(int x) {
  92.         if (x < 0) {
  93.             return 1;
  94.         } else if (x == 0) {
  95.             // 1 - exp(-mean)
  96.             return -Math.expm1(-mean);
  97.         }
  98.         return RegularizedGamma.P.value((double) x + 1, mean);
  99.     }

  100.     /** {@inheritDoc} */
  101.     @Override
  102.     public double getMean() {
  103.         return mean;
  104.     }

  105.     /**
  106.      * {@inheritDoc}
  107.      *
  108.      * <p>The variance is equal to the {@linkplain #getMean() mean}.
  109.      */
  110.     @Override
  111.     public double getVariance() {
  112.         return getMean();
  113.     }

  114.     /**
  115.      * {@inheritDoc}
  116.      *
  117.      * <p>The lower bound of the support is always 0.
  118.      *
  119.      * @return 0.
  120.      */
  121.     @Override
  122.     public int getSupportLowerBound() {
  123.         return 0;
  124.     }

  125.     /**
  126.      * {@inheritDoc}
  127.      *
  128.      * <p>The upper bound of the support is always positive infinity.
  129.      *
  130.      * @return {@link Integer#MAX_VALUE}
  131.      */
  132.     @Override
  133.     public int getSupportUpperBound() {
  134.         return Integer.MAX_VALUE;
  135.     }

  136.     /** {@inheritDoc} */
  137.     @Override
  138.     public DiscreteDistribution.Sampler createSampler(final UniformRandomProvider rng) {
  139.         // Poisson distribution sampler.
  140.         // Large means are not supported.
  141.         // See STATISTICS-35.
  142.         final double mu = getMean();
  143.         if (mu < MAX_MEAN) {
  144.             return PoissonSampler.of(rng, mu)::sample;
  145.         }
  146.         // Switch to a Gaussian approximation.
  147.         // Use a 0.5 shift to round samples to the correct integer.
  148.         final SharedStateContinuousSampler s =
  149.             GaussianSampler.of(ZigguratSampler.NormalizedGaussian.of(rng),
  150.                                mu + 0.5, Math.sqrt(mu));
  151.         return () -> {
  152.             final double x = s.sample();
  153.             return Math.max(0, (int) x);
  154.         };
  155.     }
  156. }