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}