Poisson Distribution | |||||||||||||||||||||||||
Type: | mass | ||||||||||||||||||||||||
Pdf Caption: | The horizontal axis is the index, the number of occurrences. is the expected rate of occurrences. The vertical axis is the probability of occurrences given . The function is defined only at integer values of ; the connecting lines are only guides for the eye. | ||||||||||||||||||||||||
Cdf Caption: | The horizontal axis is the index, the number of occurrences. The CDF is discontinuous at the integers of and flat everywhere else because a variable that is Poisson distributed takes on only integer values. | ||||||||||||||||||||||||
Notation: | \operatorname{Pois}(λ) | ||||||||||||||||||||||||
Parameters: | λ\in(0,infty) | ||||||||||||||||||||||||
Support: | k\inN | ||||||||||||||||||||||||
Pdf: |
| ||||||||||||||||||||||||
Cdf: |
, e-λ
, Q(\lfloork+1\rfloor,λ) k\ge0, \Gamma(x,y) \lfloork\rfloor Q | ||||||||||||||||||||||||
Mean: | λ | ||||||||||||||||||||||||
Median: | ≈ \left\lfloorλ+
-
\right\rfloor | ||||||||||||||||||||||||
Mode: | \left\lceilλ\right\rceil-1,\left\lfloorλ\right\rfloor | ||||||||||||||||||||||||
Variance: | λ | ||||||||||||||||||||||||
Skewness: |
| ||||||||||||||||||||||||
Kurtosis: |
| ||||||||||||||||||||||||
Entropy: | λl[1-log(λ)r]+e-λ
λ
log\left(2\pieλ\right)-
-
\\-
+l{O}\left(
\right)\end{align} | ||||||||||||||||||||||||
Pgf: | \exp\left[λ\left(z-1\right)\right] | ||||||||||||||||||||||||
Mgf: | \exp\left[λ\left(et-1\right)\right] | ||||||||||||||||||||||||
Char: | \exp\left[λ\left(eit-1\right)\right] | ||||||||||||||||||||||||
Fisher: |
|
In probability theory and statistics, the Poisson distribution is a discrete probability distribution that expresses the probability of a given number of events occurring in a fixed interval of time if these events occur with a known constant mean rate and independently of the time since the last event. It can also be used for the number of events in other types of intervals than time, and in dimension greater than 1 (e.g., number of events in a given area or volume).
The Poisson distribution is named after French mathematician Siméon Denis Poisson (; in French pronounced as /pwasɔ̃/). It plays an important role for discrete-stable distributions.
Under a Poisson distribution with the expectation of λ events in a given interval, the probability of k events in the same interval is:
λke-λ | |
k! |
.
A classic example used to motivate the Poisson distribution is the number of radioactive decay events during a fixed observation period.[1]
The distribution was first introduced by Siméon Denis Poisson (1781–1840) and published together with his probability theory in his work Recherches sur la probabilité des jugements en matière criminelle et en matière civile (1837). The work theorized about the number of wrongful convictions in a given country by focusing on certain random variables that count, among other things, the number of discrete occurrences (sometimes called "events" or "arrivals") that take place during a time-interval of given length. The result had already been given in 1711 by Abraham de Moivre in De Mensura Sortis seu; de Probabilitate Eventuum in Ludis a Casu Fortuito Pendentibus . This makes it an example of Stigler's law and it has prompted some authors to argue that the Poisson distribution should bear the name of de Moivre.
In 1860, Simon Newcomb fitted the Poisson distribution to the number of stars found in a unit of space.A further practical application was made by Ladislaus Bortkiewicz in 1898. Bortkiewicz showed that the frequency with which soldiers in the Prussian army were accidentally killed by horse kicks could be well modeled by a Poisson distribution..
A discrete random variable is said to have a Poisson distribution, with parameter
λ>0,
f(k;λ)=\Pr(X{=}k)=
λke-λ | |
k! |
,
k=0,1,2,\ldots
e=2.71828\ldots
The positive real number is equal to the expected value of and also to its variance.[2]
λ=\operatorname{E}(X)=\operatorname{Var}(X).
The Poisson distribution can be applied to systems with a large number of possible events, each of which is rare. The number of such events that occur during a fixed time interval is, under the right circumstances, a random number with a Poisson distribution.
The equation can be adapted if, instead of the average number of events
λ,
r
λ=rt,
P(keventsinintervalt)=
(rt)ke-rt | |
k! |
.
The Poisson distribution may be useful to model events such as:
Examples of the occurrence of random points in space are: the locations of asteroid impacts with earth (2-dimensional), the locations of imperfections in a material (3-dimensional), and the locations of trees in a forest (2-dimensional).[4]
The Poisson distribution is an appropriate model if the following assumptions are true:
If these conditions are true, then is a Poisson random variable, and the distribution of is a Poisson distribution.
The Poisson distribution is also the limit of a binomial distribution, for which the probability of success for each trial equals divided by the number of trials, as the number of trials approaches infinity (see Related distributions).
On a particular river, overflow floods occur once every 100 years on average. Calculate the probability of = 0, 1, 2, 3, 4, 5, or 6 overflow floods in a 100-year interval, assuming the Poisson model is appropriate.
Because the average event rate is one overflow flood per 100 years, = 1
P(koverflowfloodsin100years)=
λke-λ | |
k! |
=
1ke-1 | |
k! |
P(k=0overflowfloodsin100years)=
10e-1 | |
0! |
=
e-1 | |
1 |
≈ 0.368
P(k=1overflowfloodin100years)=
11e-1 | |
1! |
=
e-1 | |
1 |
≈ 0.368
P(k=2overflowfloodsin100years)=
12e-1 | |
2! |
=
e-1 | |
2 |
≈ 0.184
(overflow floods in 100 years) | ||
---|---|---|
0 | 0.368 | |
1 | 0.368 | |
2 | 0.184 | |
3 | 0.061 | |
4 | 0.015 | |
5 | 0.003 | |
6 | 0.0005 |
It has been reported that the average number of goals in a World Cup soccer match is approximately 2.5 and the Poisson model is appropriate.Because the average event rate is 2.5 goals per match, = 2.5 .
P(kgoalsinamatch)=
2.5ke-2.5 | |
k! |
P(k=0goalsinamatch)=
2.50e-2.5 | |
0! |
=
e-2.5 | |
1 |
≈ 0.082
P(k=1goalinamatch)=
2.51e-2.5 | |
1! |
=
2.5e-2.5 | |
1 |
≈ 0.205
P(k=2goalsinamatch)=
2.52e-2.5 | |
2! |
=
6.25e-2.5 | |
2 |
≈ 0.257
(goals in a World Cup soccer match) | ||
---|---|---|
0 | 0.082 | |
1 | 0.205 | |
2 | 0.257 | |
3 | 0.213 | |
4 | 0.133 | |
5 | 0.067 | |
6 | 0.028 | |
7 | 0.010 |
Suppose that astronomers estimate that large meteorites (above a certain size) hit the earth on average once every 100 years (event per 100 years), and that the number of meteorite hits follows a Poisson distribution. What is the probability of meteorite hits in the next 100 years?
P(k=0meteoriteshitinnext100years)=
10e-1 | |
0! |
=
1 | |
e |
≈ 0.37.
Under these assumptions, the probability that no large meteorites hit the earth in the next 100 years is roughly 0.37. The remaining is the probability of 1, 2, 3, or more large meteorite hits in the next 100 years.In an example above, an overflow flood occurred once every 100 years The probability of no overflow floods in 100 years was roughly 0.37, by the same calculation.
In general, if an event occurs on average once per interval ( = 1), and the events follow a Poisson distribution, then In addition, as shown in the table for overflow floods.
The number of students who arrive at the student union per minute will likely not follow a Poisson distribution, because the rate is not constant (low rate during class time, high rate between class times) and the arrivals of individual students are not independent (students tend to come in groups). The non-constant arrival rate may be modeled as a mixed Poisson distribution, and the arrival of groups rather than individual students as a compound Poisson process.
The number of magnitude 5 earthquakes per year in a country may not follow a Poisson distribution, if one large earthquake increases the probability of aftershocks of similar magnitude.
Examples in which at least one event is guaranteed are not Poisson distributed; but may be modeled using a zero-truncated Poisson distribution.
Count distributions in which the number of intervals with zero events is higher than predicted by a Poisson model may be modeled using a zero-inflated model.
\lfloorλ\rfloor,
Bounds for the median (
\nu
The higher non-centered moments, of the Poisson distribution, are Touchard polynomials in : where the braces denote Stirling numbers of the second kind. In other words,When the expected value is set to λ = 1, Dobinski's formula implies that the ‑th moment is equal to the number of partitions of a set of size .
A simple upper bound is:[5]
If
Xi\sim\operatorname{Pois}(λi)
i=1,...c,n
It is a maximum-entropy distribution among the set of generalized binomial distributions
Bn(λ)
λ
n → infty
P=\operatorname{Pois}(λ)
P0=\operatorname{Pois}(λ0)
λ\geq1
Y\sim\operatorname{Pois}(λ)
\Pr(Y\geqE[Y])\geq
1 | |
2 |
\Pr(Y\leqE[Y])\geq
1 | |
2 |
.
X\sim\operatorname{Pois}(λ)
where
\operatorname{D}KL(Q\parallelP)
Q=\operatorname{Pois}(x)
P=\operatorname{Pois}(λ)
X\sim\operatorname{Pois}(λ)
\Phi(x)
\operatorname{D}KL(Q-\parallelP)
Q-=\operatorname{Pois}(k)
P=\operatorname{Pois}(λ)
\operatorname{D}KL(Q+\parallelP)
Q+=\operatorname{Pois}(k+1)
P
Let
X\sim\operatorname{Pois}(λ)
Y\sim\operatorname{Pois}(\mu)
λ<\mu,
The upper bound is proved using a standard Chernoff bound.
The lower bound can be proved by noting that
P(X-Y\geq0\midX+Y=i)
D
X+Y\sim\operatorname{Pois}(λ+\mu),
The Poisson distribution can be derived as a limiting case to the binomial distribution as the number of trials goes to infinity and the expected number of successes remains fixed — see law of rare events below. Therefore, it can be used as an approximation of the binomial distribution if is sufficiently large and p is sufficiently small. The Poisson distribution is a good approximation of the binomial distribution if is at least 20 and p is smaller than or equal to 0.05, and an excellent approximation if ≥ 100 and ≤ 10. Letting
FB
FP
λ\leq1
\tfrac{λ}{n}
whose generating function is:,
(n) p k =\binomnk\left(
λ n \right)k\left(1{-}
λ n \right)
Taking the limit as n increases to infinity (with x fixed) and applying the product limit definition of the exponential function, this reduces to the generating function of the Poisson distribution:P(n)
n (x)=\sum k=0
(n) p k xk=\left(1-
λ + n
λ n x\right)n.
\limn\toinftyP(n)(x)=\limn\toinfty\left(1{+}\tfrac{λ(x-1)}{n}\right)n =eλ(x-1)=
infty \sum k=0 e-λ
λk k! xk.
X1\simPois(λ1)
X2\simPois(λ2)
Y=X1-X2
X1\simPois(λ1)
X2\simPois(λ2)
X1
X1+X2
X1+X2=k,
X1|X1+X2=k\simBinom(k,λ1/(λ1+λ2)).
given
n | |
\sum | |
j=1 |
Xj=k,
Xi|\sum
n | |
j=1 |
Xj=k\simBinom\left(k,
λi | |||||||||
|
\right).
\{Xi\}\simMultinom\left(k,\left\{
λi | ||||||
|
\right\}\right).
X\simPois(λ)
Y
Y\mid(X=k)\simBinom(k,p),
Y\simPois(λ ⋅ p).
\{X=k\},
\{Yi\}
\{Yi\}\mid(X=k)\simMultinom\left(k,pi\right),
Yi
Yi\simPois(λ ⋅ pi),\rho(Yi,Yj)=0.
\sqrt{λ}
X\simPois(λ),
λ
Assume
X1\sim\operatorname{Pois}(λ1),X2\sim\operatorname{Pois}(λ2),...,Xn\sim\operatorname{Pois}(λn)
λ1+λ2+...+λn=1,
(X1,X2,...,Xn)
(X1,X2,...,Xn)\sim\operatorname{Mult}(N,λ1,λ2,...,λn)
N=X1+X2+...Xn.
This means, among other things, that for any nonnegative function
f(x1,x2,...,xn),
(Y1,Y2,...,Yn)\sim\operatorname{Mult}(m,p)
(X1,X2,...,Xn)\sim\operatorname{Pois}(p).
The factor of
e\sqrt{m}
f
This distribution has been extended to the bivariate case. The generating function for this distribution is
with
The marginal distributions are Poisson(θ1) and Poisson(θ2) and the correlation coefficient is limited to the range
A simple way to generate a bivariate Poisson distribution
X1,X2
Y1,Y2,Y3
λ1,λ2,λ3
X1=Y1+Y3,X2=Y2+Y3.
The free Poisson distribution[12] with jump size
\alpha
λ
In other words, let
XN
XN
\alpha
X1,X2,\ldots
N\toinfty
X1+ … +XN
λ,\alpha.
This definition is analogous to one of the ways in which the classical Poisson distribution is obtained from a (classical) Poisson process.
The measure associated to the free Poisson law is given by[13] where and has support
[\alpha(1-\sqrt{λ})2,\alpha(1+\sqrt{λ})2].
This law also arises in random matrix theory as the Marchenko–Pastur law. Its free cumulants are equal to
n. | |
\kappa | |
n=λ\alpha |
We give values of some important transforms of the free Poisson law; the computation can be found in e.g. in the book Lectures on the Combinatorics of Free Probability by A. Nica and R. Speicher[14]
The R-transform of the free Poisson law is given by
The Cauchy transform (which is the negative of the Stieltjes transformation) is given by
The S-transform is given byin the case that
\alpha=1.
Poisson's probability mass function
f(k;λ)
(k+1)
ak{N}\alpha(\nu)
\alpha=1/\left(k+1\right),
Wk+1(x)
k+1.
See also: Poisson regression.
Given a sample of measured values
ki\in\{0,1,...\},
\widehat{λ} | ||||
|
n | |
\sum | |
i=1 |
ki .
Since each observation has expectation so does the sample mean. Therefore, the maximum likelihood estimate is an unbiased estimator of . It is also an efficient estimator since its variance achieves the Cramér–Rao lower bound (CRLB).[16] Hence it is minimum-variance unbiased. Also it can be proven that the sum (and hence the sample mean as it is a one-to-one function of the sum) is a complete and sufficient statistic for .
To prove sufficiency we may use the factorization theorem. Consider partitioning the probability mass function of the joint Poisson distribution for the sample into two parts: one that depends solely on the sample
x
h(x)
λ
x
T(x).
T(x)
λ.
| |||||||||||
P(x)=\prod | = | ||||||||||
i=1 |
1 | |||||||||
|
x
| ||||||||||
λ |
e-nλ
The first term
h(x)
x
g(T(x)|λ)
T(x)
To find the parameter that maximizes the probability function for the Poisson population, we can use the logarithm of the likelihood function:
\begin{align} \ell(λ)&=ln
n | |
\prod | |
i=1 |
f(ki\midλ)\\ &=
n | ||
\sum | ln\left( | |
i=1 |
| |||||||||||
ki! |
\right)\\ &=-nλ+
n | |
\left(\sum | |
i=1 |
ki\right)ln(λ)-
n | |
\sum | |
i=1 |
ln(ki!). \end{align}
We take the derivative of
\ell
d | |
dλ |
\ell(λ)=0\iff-n+
n | |
\left(\sum | |
i=1 |
ki\right)
1 | |
λ |
=0.
Solving for gives a stationary point.
λ=
| ||||||||||
n |
So is the average of the i values. Obtaining the sign of the second derivative of L at the stationary point will determine what kind of extreme value is.
\partial2\ell | |
\partialλ2 |
=-λ-2
n | |
\sum | |
i=1 |
ki
Evaluating the second derivative at the stationary point gives:
\partial2\ell | |
\partialλ2 |
=-
n2 | |||||||||
|
which is the negative of times the reciprocal of the average of the ki. This expression is negative when the average is positive. If this is satisfied, then the stationary point maximizes the probability function.
For completeness, a family of distributions is said to be complete if and only if
E(g(T))=0
Pλ(g(T)=0)=1
λ.
Xi
Po(λ),
infty | ||
E(g(T))=\sum | g(t) | |
t=0 |
(nλ)te-nλ | |
t! |
=0
For this equality to hold,
g(t)
t
λ.
E(g(T))=0
λ
Pλ(g(T)=0)=1,
The confidence interval for the mean of a Poisson distribution can be expressed using the relationship between the cumulative distribution functions of the Poisson and chi-squared distributions. The chi-squared distribution is itself closely related to the gamma distribution, and this leads to an alternative expression. Given an observation from a Poisson distribution with mean μ, a confidence interval for μ with confidence level is
\tfrac{1}{2}\chi2(\alpha/2;2k)\le\mu\le\tfrac{1}{2}\chi2(1-\alpha/2;2k+2),
or equivalently,
F-1(\alpha/2;k,1)\le\mu\leF-1(1-\alpha/2;k+1,1),
where
\chi2(p;n)
F-1(p;n,1)
When quantiles of the gamma distribution are not available, an accurate approximation to this exact interval has been proposed (based on the Wilson–Hilferty transformation):
k\left(1-
1 | |
9k |
-
z\alpha/2 | |
3\sqrt{k |
z\alpha/2
For application of these formulae in the same context as above (given a sample of measured values i each drawn from a Poisson distribution with mean), one would set
n | |
k=\sum | |
i=1 |
ki,
In Bayesian inference, the conjugate prior for the rate parameter of the Poisson distribution is the gamma distribution. Let
λ\simGamma(\alpha,\beta)
denote that is distributed according to the gamma density g parameterized in terms of a shape parameter α and an inverse scale parameter β:
g(λ\mid\alpha,\beta)=
\beta\alpha | |
\Gamma(\alpha) |
λ\alpha-1 e-\betaλ forλ>0.
Then, given the same sample of measured values i as before, and a prior of Gamma(α, β), the posterior distribution is
λ\simGamma\left(\alpha+
n | |
\sum | |
i=1 |
ki,\beta+n\right).
Note that the posterior mean is linear and is given by
E[λ|k1,\ldots,kn]=
| |||||||||
\beta+n |
.
L2
The posterior mean E[{{mvar|λ}}] approaches the maximum likelihood estimate
\widehat{λ}MLE
\alpha\to0,\beta\to0,
The posterior predictive distribution for a single additional observation is a negative binomial distribution, sometimes called a gamma–Poisson distribution.
Suppose
X1,X2,...,Xp
p
λi,
i=1,...,p,
p>1,
{\hatλ}i=Xi
In this case, a family of minimax estimators is given for any
0<c\leq2(p-1)
b\geq(p-2+p-1)
{\hatλ}i=\left(1-
c | ||||||||
|
\right)Xi, i=1,...,p.
Some applications of the Poisson distribution to count data (number of events):
More examples of counting events that may be modelled as Poisson processes include:
In probabilistic number theory, Gallagher showed in 1976 that, if a certain version of the unproved prime r-tuple conjecture holds, then the counts of prime numbers in short intervals would obey a Poisson distribution.
See main article: Poisson limit theorem.
The rate of an event is related to the probability of an event occurring in some small subinterval (of time, space or otherwise). In the case of the Poisson distribution, one assumes that there exists a small enough subinterval for which the probability of an event occurring twice is "negligible". With this assumption one can derive the Poisson distribution from the Binomial one, given only the information of expected number of total events in the whole interval.
Let the total number of events in the whole interval be denoted by
λ.
n
I1,...,In
n>λ
λ/n.
Now we assume that the occurrence of an event in the whole interval can be seen as a sequence of Bernoulli trials, where the
i
Ii
λ/n.
n
λ,
rm{B}(n,λ/n).
n
In this case the binomial distribution converges to what is known as the Poisson distribution by the Poisson limit theorem.
In several of the above examples — such as, the number of mutations in a given sequence of DNA—the events being counted are actually the outcomes of discrete trials, and would more precisely be modelled using the binomial distribution, that is
In such cases is very large and is very small (and so the expectation is of intermediate magnitude). Then the distribution may be approximated by the less cumbersome Poisson distribution
This approximation is sometimes known as the law of rare events, since each of the individual Bernoulli events rarely occurs.
The name "law of rare events" may be misleading because the total count of success events in a Poisson process need not be rare if the parameter is not small. For example, the number of telephone calls to a busy switchboard in one hour follows a Poisson distribution with the events appearing frequent to the operator, but they are rare from the point of view of the average member of the population who is very unlikely to make a call to that switchboard in that hour.
The variance of the binomial distribution is 1 − p times that of the Poisson distribution, so almost equal when p is very small.
The word law is sometimes used as a synonym of probability distribution, and convergence in law means convergence in distribution. Accordingly, the Poisson distribution is sometimes called the "law of small numbers" because it is the probability distribution of the number of occurrences of an event that happens rarely but has very many opportunities to happen. The Law of Small Numbers is a book by Ladislaus Bortkiewicz about the Poisson distribution, published in 1898.
See main article: Poisson point process.
The Poisson distribution arises as the number of points of a Poisson point process located in some finite region. More specifically, if D is some region space, for example Euclidean space Rd, for which |D|, the area, volume or, more generally, the Lebesgue measure of the region is finite, and if denotes the number of points in D, then
P(N(D)=k)= | (λ|D|)ke-λ|D| |
k! |
.
Poisson regression and negative binomial regression are useful for analyses where the dependent (response) variable is the count of the number of events or occurrences in an interval.
The Luria–Delbrück experiment tested against the hypothesis of Lamarckian evolution, which should result in a Poisson distribution.
Katz and Miledi measured the membrane potential with and without the presence of acetylcholine (ACh).[18] When ACh is present, ion channels on the membrane would be open randomly at a small fraction of the time. As there are a large number of ion channels each open for a small fraction of the time, the total number of ion channels open at any moment is Poisson distributed. When ACh is not present, effectively no ion channels are open. The membrane potential is
V=NopenVion+V0+Vnoise
8.5 x 10-3 V,(29.2 x 10-6 V)2
Vion=10-7 V
During each cellular replication event, the number of mutations is roughly Poisson distributed. For example, the HIV virus has 10,000 base pairs, and has a mutation rate of about 1 per 30,000 base pairs, meaning the number of mutations per replication event is distributed as
Pois(1/3)
\sigmak=\sqrt{λ}.
The correlation of the mean and standard deviation in counting independent discrete occurrences is useful scientifically. By monitoring how the fluctuations vary with the mean signal, one can estimate the contribution of a single occurrence, even if that contribution is too small to be detected directly. For example, the charge e on an electron can be estimated by correlating the magnitude of an electric current with its shot noise. If N electrons pass a point in a given time t on the average, the mean current is
I=eN/t
\sigmaI=e\sqrt{N}/t
e
2/I. | |
t\sigma | |
I |
An everyday example is the graininess that appears as photographs are enlarged; the graininess is due to Poisson fluctuations in the number of reduced silver grains, not to the individual grains themselves. By correlating the graininess with the degree of enlargement, one can estimate the contribution of an individual grain (which is otherwise too small to be seen unaided).
In causal set theory the discrete elements of spacetime follow a Poisson distribution in the volume.
The Poisson distribution also appears in quantum mechanics, especially quantum optics. Namely, for a quantum harmonic oscillator system in a coherent state, the probability of measuring a particular energy level has a Poisson distribution.
The Poisson distribution poses two different tasks for dedicated software libraries: evaluating the distribution
P(k;λ)
Computing
P(k;λ)
k
λ
P(k;λ)
f(k;λ)=\exp\left[klnλ-λ-ln\Gamma(k+1)\right],
lgamma
function in the C standard library (C99 version) or R, the gammaln
function in MATLAB or SciPy, or the log_gamma
function in Fortran 2008 and later.Some computing languages provide built-in functions to evaluate the Poisson distribution, namely
dpois(x, lambda)
;POISSON(x, mean, cumulative)
, with a flag to specify the cumulative distribution;PoissonDistribution[<math>\lambda</math>]
,[20] bivariate Poisson distribution as MultivariatePoissonDistribution[<math>\theta_{12},</math>{ <math>\theta_1 - \theta_{12},</math> <math>\theta_2 - \theta_{12}</math>}]
,.[21]The less trivial task is to draw integer random variate from the Poisson distribution with given
λ.
Solutions are provided by:
rpois(n, lambda)
;A simple algorithm to generate random Poisson-distributed numbers (pseudo-random number sampling) has been given by Knuth:
algorithm poisson random number (Knuth): init: Let L ← e−λ, k ← 0 and p ← 1. do: k ← k + 1. Generate uniform random number u in [0,1] and let p ← p × u. while p > L. return k − 1.
The complexity is linear in the returned value, which is on average. There are many other algorithms to improve this. Some are given in Ahrens & Dieter, see below.
For large values of, the value of = e− may be so small that it is hard to represent. This can be solved by a change to the algorithm which uses an additional parameter STEP such that e−STEP does not underflow:
algorithm poisson random number (Junhao, based on Knuth): init: Let Left ←, k ← 0 and p ← 1. do: k ← k + 1. Generate uniform random number u in (0,1) and let p ← p × u. while p < 1 and Left > 0: if Left > STEP: p ← p × eSTEP Left ← Left − STEP else: p ← p × eLeft Left ← 0 while p > 1. return k − 1.
The choice of STEP depends on the threshold of overflow. For double precision floating point format the threshold is near e700, so 500 should be a safe STEP.
Other solutions for large values of include rejection sampling and using Gaussian approximation.
Inverse transform sampling is simple and efficient for small values of, and requires only one uniform random number u per sample. Cumulative probabilities are examined in turn until one exceeds u.
algorithm Poisson generator based upon the inversion by sequential search: init: Let x ← 0, p ← e−λ, s ← p. Generate uniform random number u in [0,1]. while u > s do: x ← x + 1. p ← p × / x. s ← s + p. return x.