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    /** Upper bound on the mean to use the PoissonSampler. */
041    private static final double MAX_MEAN = 0.5 * Integer.MAX_VALUE;
042    /** Mean of the distribution. */
043    private final double mean;
044
045    /**
046     * @param mean Poisson mean.
047     * probabilities.
048     */
049    private PoissonDistribution(double mean) {
050        this.mean = mean;
051    }
052
053    /**
054     * Creates a Poisson distribution.
055     *
056     * @param mean Poisson mean.
057     * @return the distribution
058     * @throws IllegalArgumentException if {@code mean <= 0}.
059     */
060    public static PoissonDistribution of(double mean) {
061        if (mean <= 0) {
062            throw new DistributionException(DistributionException.NOT_STRICTLY_POSITIVE, mean);
063        }
064        return new PoissonDistribution(mean);
065    }
066
067    /** {@inheritDoc} */
068    @Override
069    public double probability(int x) {
070        return Math.exp(logProbability(x));
071    }
072
073    /** {@inheritDoc} */
074    @Override
075    public double logProbability(int x) {
076        if (x < 0) {
077            return Double.NEGATIVE_INFINITY;
078        } else if (x == 0) {
079            return -mean;
080        }
081        return -SaddlePointExpansionUtils.getStirlingError(x) -
082              SaddlePointExpansionUtils.getDeviancePart(x, mean) -
083              Constants.HALF_LOG_TWO_PI - 0.5 * Math.log(x);
084    }
085
086    /** {@inheritDoc} */
087    @Override
088    public double cumulativeProbability(int x) {
089        if (x < 0) {
090            return 0;
091        } else if (x == 0) {
092            return Math.exp(-mean);
093        }
094        return RegularizedGamma.Q.value((double) x + 1, mean);
095    }
096
097    /** {@inheritDoc} */
098    @Override
099    public double survivalProbability(int x) {
100        if (x < 0) {
101            return 1;
102        } else if (x == 0) {
103            // 1 - exp(-mean)
104            return -Math.expm1(-mean);
105        }
106        return RegularizedGamma.P.value((double) x + 1, mean);
107    }
108
109    /** {@inheritDoc} */
110    @Override
111    public double getMean() {
112        return mean;
113    }
114
115    /**
116     * {@inheritDoc}
117     *
118     * <p>The variance is equal to the {@linkplain #getMean() mean}.
119     */
120    @Override
121    public double getVariance() {
122        return getMean();
123    }
124
125    /**
126     * {@inheritDoc}
127     *
128     * <p>The lower bound of the support is always 0.
129     *
130     * @return 0.
131     */
132    @Override
133    public int getSupportLowerBound() {
134        return 0;
135    }
136
137    /**
138     * {@inheritDoc}
139     *
140     * <p>The upper bound of the support is always positive infinity.
141     *
142     * @return {@link Integer#MAX_VALUE}
143     */
144    @Override
145    public int getSupportUpperBound() {
146        return Integer.MAX_VALUE;
147    }
148
149    /** {@inheritDoc} */
150    @Override
151    public DiscreteDistribution.Sampler createSampler(final UniformRandomProvider rng) {
152        // Poisson distribution sampler.
153        // Large means are not supported.
154        // See STATISTICS-35.
155        final double mu = getMean();
156        if (mu < MAX_MEAN) {
157            return PoissonSampler.of(rng, mu)::sample;
158        }
159        // Switch to a Gaussian approximation.
160        // Use a 0.5 shift to round samples to the correct integer.
161        final SharedStateContinuousSampler s =
162            GaussianSampler.of(ZigguratSampler.NormalizedGaussian.of(rng),
163                               mu + 0.5, Math.sqrt(mu));
164        return () -> {
165            final double x = s.sample();
166            return Math.max(0, (int) x);
167        };
168    }
169}