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.statistics.distribution;
018
019import org.apache.commons.numbers.gamma.RegularizedGamma;
020import org.apache.commons.rng.UniformRandomProvider;
021import org.apache.commons.rng.sampling.distribution.GaussianSampler;
022import org.apache.commons.rng.sampling.distribution.PoissonSampler;
023import org.apache.commons.rng.sampling.distribution.SharedStateContinuousSampler;
024import org.apache.commons.rng.sampling.distribution.ZigguratSampler;
025
026/**
027 * Implementation of the Poisson distribution.
028 *
029 * <p>The probability mass function of \( X \) is:
030 *
031 * <p>\[ f(k; \lambda) = \frac{\lambda^k e^{-k}}{k!} \]
032 *
033 * <p>for \( \lambda \in (0, \infty) \) the mean and
034 * \( k \in \{0, 1, 2, \dots\} \) the number of events.
035 *
036 * @see <a href="https://en.wikipedia.org/wiki/Poisson_distribution">Poisson distribution (Wikipedia)</a>
037 * @see <a href="https://mathworld.wolfram.com/PoissonDistribution.html">Poisson distribution (MathWorld)</a>
038 */
039public final class PoissonDistribution extends AbstractDiscreteDistribution {
040    /** 0.5 * ln(2 * pi). Computed to 25-digits precision. */
041    private static final double HALF_LOG_TWO_PI = 0.9189385332046727417803297;
042    /** Upper bound on the mean to use the PoissonSampler. */
043    private static final double MAX_MEAN = 0.5 * Integer.MAX_VALUE;
044    /** Mean of the distribution. */
045    private final double mean;
046
047    /**
048     * @param mean Poisson mean.
049     * probabilities.
050     */
051    private PoissonDistribution(double mean) {
052        this.mean = mean;
053    }
054
055    /**
056     * Creates a Poisson distribution.
057     *
058     * @param mean Poisson mean.
059     * @return the distribution
060     * @throws IllegalArgumentException if {@code mean <= 0}.
061     */
062    public static PoissonDistribution of(double mean) {
063        if (mean <= 0) {
064            throw new DistributionException(DistributionException.NOT_STRICTLY_POSITIVE, mean);
065        }
066        return new PoissonDistribution(mean);
067    }
068
069    /** {@inheritDoc} */
070    @Override
071    public double probability(int x) {
072        return Math.exp(logProbability(x));
073    }
074
075    /** {@inheritDoc} */
076    @Override
077    public double logProbability(int x) {
078        if (x < 0) {
079            return Double.NEGATIVE_INFINITY;
080        } else if (x == 0) {
081            return -mean;
082        }
083        return -SaddlePointExpansionUtils.getStirlingError(x) -
084              SaddlePointExpansionUtils.getDeviancePart(x, mean) -
085              HALF_LOG_TWO_PI - 0.5 * Math.log(x);
086    }
087
088    /** {@inheritDoc} */
089    @Override
090    public double cumulativeProbability(int x) {
091        if (x < 0) {
092            return 0;
093        } else if (x == 0) {
094            return Math.exp(-mean);
095        }
096        return RegularizedGamma.Q.value((double) x + 1, mean);
097    }
098
099    /** {@inheritDoc} */
100    @Override
101    public double survivalProbability(int x) {
102        if (x < 0) {
103            return 1;
104        } else if (x == 0) {
105            // 1 - exp(-mean)
106            return -Math.expm1(-mean);
107        }
108        return RegularizedGamma.P.value((double) x + 1, mean);
109    }
110
111    /** {@inheritDoc} */
112    @Override
113    public double getMean() {
114        return mean;
115    }
116
117    /**
118     * {@inheritDoc}
119     *
120     * <p>The variance is equal to the {@link #getMean() mean}.
121     */
122    @Override
123    public double getVariance() {
124        return getMean();
125    }
126
127    /**
128     * {@inheritDoc}
129     *
130     * <p>The lower bound of the support is always 0.
131     *
132     * @return 0.
133     */
134    @Override
135    public int getSupportLowerBound() {
136        return 0;
137    }
138
139    /**
140     * {@inheritDoc}
141     *
142     * <p>The upper bound of the support is always positive infinity.
143     *
144     * @return {@link Integer#MAX_VALUE}
145     */
146    @Override
147    public int getSupportUpperBound() {
148        return Integer.MAX_VALUE;
149    }
150
151    /** {@inheritDoc} */
152    @Override
153    public DiscreteDistribution.Sampler createSampler(final UniformRandomProvider rng) {
154        // Poisson distribution sampler.
155        // Large means are not supported.
156        // See STATISTICS-35.
157        final double mu = getMean();
158        if (mu < MAX_MEAN) {
159            return PoissonSampler.of(rng, mu)::sample;
160        }
161        // Switch to a Gaussian approximation.
162        // Use a 0.5 shift to round samples to the correct integer.
163        final SharedStateContinuousSampler s =
164            GaussianSampler.of(ZigguratSampler.NormalizedGaussian.of(rng),
165                               mu + 0.5, Math.sqrt(mu));
166        return () -> {
167            final double x = s.sample();
168            return Math.max(0, (int) x);
169        };
170    }
171}