What is the original derivation of the Poisson distribution?

Question

I am learning the Poisson distribution. I understand it, but its probability mass function is not natural to me. I think its probability mass function seems to be derived from somewhere with more contexts (or from one application). Where is the original source?

Answer 1

I think it first appeared in Poisson's "Recherches sur la probabilité des jugements en matière criminelle et en matière civile' (researches into the probabilities of judgements in criminal and civil cases).

He did not define it through the PMF but as a limit of a Binomial distribution (using the law of large numbers which he calls "Bernoulli's theorem")

At some point he obtains the formula

$$ P = \left(1 + \omega + \frac{\omega^2}{1 \cdot 2} + \dots + \frac{\omega^n}{1\cdot 2 \cdots n} \right) e^{-\omega} $$

for the probability of obtaining at most $n$ successess out of $\mu$ trials in which there is a probability $\frac{\omega}{\mu}$ of success (for $\mu$ "very large").

This is the cumulative distribution function of a Poisson distribution with parameter $\omega$ evaluated at $n$.

You can see the formula above here page 206 (it's in French).

Answer 2

Thanks to @Glen-b, and @Xi'an for historical references, as requested in the Question. However, looking at historical references is not always a good way to clarify a now-familiar concept. The archaic language and setting may be confusing, and earlier discussions will not mention clarifying insights that came later.

Perhaps an intuitive insight comes from comparing distributions $\mathsf{Binom}(n=100, p = .2),$ $\mathsf{Binom}(n=1000, p = .02),$ and $\mathsf{Pois}(\lambda = 20).$ The former two distributions may be considered as indications of the limit of a binomial distribution as $n \rightarrow \infty,$ while $\mu = \lambda = np$ remains constant.

k = 10:30
pb.100  = round(dbinom(k, 100, .2) ,3)
pb.1000 = round(dbinom(k, 1000, .02), 3)
p.pois =  round(dpois(k, 20), 3)
cbind(k, pb.100, pb.1000, p.pois)

       k pb.100 pb.1000 p.pois
 [1,] 10  0.003   0.006  0.006
 [2,] 11  0.007   0.010  0.011
 [3,] 12  0.013   0.017  0.018
 [4,] 13  0.022   0.027  0.027
 [5,] 14  0.034   0.038  0.039

 [6,] 15  0.048   0.051  0.052
 [7,] 16  0.064   0.065  0.065
 [8,] 17  0.079   0.076  0.076
 [9,] 18  0.091   0.085  0.084
[10,] 19  0.098   0.090  0.089

[11,] 20  0.099   0.090  0.089
[12,] 21  0.095   0.085  0.085
[13,] 22  0.085   0.078  0.077
[14,] 23  0.072   0.067  0.067
[15,] 24  0.058   0.056  0.056

[16,] 25  0.044   0.045  0.045
[17,] 26  0.032   0.034  0.034
[18,] 27  0.022   0.025  0.025
[19,] 28  0.014   0.018  0.018
[20,] 29  0.009   0.012  0.013
[21,] 30  0.005   0.008  0.008

The vertical resolution of the plot below is about $0.005$ or $0.01.$

By the CLT, all three discrete distributions are well approximated by $\mathsf{Norm}(\mu = 20, \sigma = \sqrt{20}).$ But plotting that density function would clutter the figure.

R code for figure:

k = 0:35
PB.100  = dbinom(k, 100, .2)
PB.1000 = dbinom(k, 1000, .02)
P.pois  = dpois(k, 20)

hdr = "PDFs of BINOM(100,.2) [blue], 
        BINOM(1000,.02) [cyan], and POIS(20)"
plot(k-.2, PB.100,, type="h", col="blue", 
      ylab="PDF", xlab="x", main=hdr)
lines(k, PB.1000, type="h", col="cyan3")
lines(k+.2, P.pois, type="h")
 abline(h=0, col="green2")