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
19 import org.apache.commons.numbers.gamma.RegularizedGamma;
20 import org.apache.commons.rng.UniformRandomProvider;
21 import org.apache.commons.rng.sampling.distribution.GaussianSampler;
22 import org.apache.commons.rng.sampling.distribution.PoissonSampler;
23 import org.apache.commons.rng.sampling.distribution.SharedStateContinuousSampler;
24 import org.apache.commons.rng.sampling.distribution.ZigguratSampler;
25
26 /**
27 * Implementation of the Poisson distribution.
28 *
29 * <p>The probability mass function of \( X \) is:
30 *
31 * <p>\[ f(k; \lambda) = \frac{\lambda^k e^{-k}}{k!} \]
32 *
33 * <p>for \( \lambda \in (0, \infty) \) the mean and
34 * \( k \in \{0, 1, 2, \dots\} \) the number of events.
35 *
36 * @see <a href="https://en.wikipedia.org/wiki/Poisson_distribution">Poisson distribution (Wikipedia)</a>
37 * @see <a href="https://mathworld.wolfram.com/PoissonDistribution.html">Poisson distribution (MathWorld)</a>
38 */
39 public final class PoissonDistribution extends AbstractDiscreteDistribution {
40 /** Upper bound on the mean to use the PoissonSampler. */
41 private static final double MAX_MEAN = 0.5 * Integer.MAX_VALUE;
42 /** Mean of the distribution. */
43 private final double mean;
44
45 /**
46 * @param mean Poisson mean.
47 * probabilities.
48 */
49 private PoissonDistribution(double mean) {
50 this.mean = mean;
51 }
52
53 /**
54 * Creates a Poisson distribution.
55 *
56 * @param mean Poisson mean.
57 * @return the distribution
58 * @throws IllegalArgumentException if {@code mean <= 0}.
59 */
60 public static PoissonDistribution of(double mean) {
61 if (mean <= 0) {
62 throw new DistributionException(DistributionException.NOT_STRICTLY_POSITIVE, mean);
63 }
64 return new PoissonDistribution(mean);
65 }
66
67 /** {@inheritDoc} */
68 @Override
69 public double probability(int x) {
70 return Math.exp(logProbability(x));
71 }
72
73 /** {@inheritDoc} */
74 @Override
75 public double logProbability(int x) {
76 if (x < 0) {
77 return Double.NEGATIVE_INFINITY;
78 } else if (x == 0) {
79 return -mean;
80 }
81 return -SaddlePointExpansionUtils.getStirlingError(x) -
82 SaddlePointExpansionUtils.getDeviancePart(x, mean) -
83 Constants.HALF_LOG_TWO_PI - 0.5 * Math.log(x);
84 }
85
86 /** {@inheritDoc} */
87 @Override
88 public double cumulativeProbability(int x) {
89 if (x < 0) {
90 return 0;
91 } else if (x == 0) {
92 return Math.exp(-mean);
93 }
94 return RegularizedGamma.Q.value((double) x + 1, mean);
95 }
96
97 /** {@inheritDoc} */
98 @Override
99 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 }