The Number of Exponential Summands in a Fixed Interval is Poisson

Question

What is the most clever way to prove that the number of independent exponential summands in a fixed interval is distributed as a Poisson random variable? I can do it one way, but I would like to know if there is another way that gets more style points.

Let $S_1, S_2, \ldots \overset{iid}{\sim} \text{Exponential}(\mu)$. Just to be clear, the densities are each $f_S(s) = \mu e^{-\mu s}$. Now, for $t > 0$, define $K_t = \{j : S_1 + \cdots S_j < t < S_1 + \cdots + S_j + S_{j+1}\}$.

The way I can do it is by directly integrating the joint density of the exponentials. For example, $$P(K_t = 1) = P(S_1 < t < S_1 + S_2) = \int_0^t \int_{t-s_1}^{\infty} \mu e^{-\mu s_1}\ \mu e^{-\mu s_2}\ ds_2\ ds_1 = (\mu t) e^{\mu t}.$$ But is there another way?

Edit: This is a well-known result of basic Poisson processes. This fact was motivation for the question. So, ideally, answers won't rely on any results from Poisson or any other stochastic processes.

Answer 1

I propose taking, as a point of departure, the concept of a homogeneous Poisson process. This is a point process on the line (often thought of, and referred to, as a "time" line). The realizations are sets of points. Almost surely, any bounded set of real numbers will contain only finitely many points.

The fundamental properties enjoyed by this process, the ones I will use repeatedly in the analysis, are

1. (Independence) The outcomes in any two disjoint sets are independent.

2. (Homogeneity) The expected number of points in any measurable set $\mathcal{A}$ with finite measure $|\mathcal A|$ is directly proportional to $|\mathcal{A}|$. The constant of proportionality, $\lambda$, is nonzero.

Everything flows from these properties, as we will see.

Waiting times

Let's study the "waiting times" of this process. Given a start time $s$ and elapsed duration $t \ge 0$, let $S(s, t)$ be the chance that no points occur between $s$ and $s+t$: that is, within the interval $(s, s+t]$. Consider two adjacent intervals, one from $r$ to $r+s$ and another from $r+s$ to $r+s+t$. By independence, that chance of no point being in their union is the product of the chance that no point is in the first and the chance that no point is in the second:

$$S(r, s+t) = S(r,s) S(r+s,t).\tag{1}$$

By homogeneity, these chances remain the same when we slide the intervals around. That is, for any $s$, $S(r, t) = S(r+s, t)$. In particular, we may always take $r=-s$ to obtain

$$S(r, t) = S(0,t) = S(t)$$

for all $r$, allowing us to drop the explicit dependence on $r$ in the notation. Plugging this into $(1)$ gives

$$S(s+t) = S(s)S(t).\tag{2}$$

Homogeneity makes it obvious $S$ must be continuous (actually, differentiable). It is well known that the only solutions to $(2)$ are exponential. One simple way to see this is to consider that the logarithm of $S$ is linear and, since $S(0)=1$, there consequently must be some number $\kappa$ for which

$$\log(S(t)) = \kappa t.$$

Since $S$ must decrease as time goes on, $\kappa \lt 0$. Ergo, all solutions are of the form

$$S(t) = e^{-\kappa t}.$$

It is the chance that no points occur within any specified interval of length $t$.

Exponential summands

Fix an interval; by virtue of homogeneity, we may assume it starts at $0$ and ends (say) at $b$. Almost surely there are only finitely many points of the process in the interval $(0,b]$, allowing us to order them $0 \lt t_1 \lt t_2 \lt \cdots \lt t_n$. $t_1$ is a realization of a random variable $T_1$ governed by $S$: that is,

$$\Pr(T_1 \le t) = 1 - \Pr(T_1 \gt t) = 1 - S(t) = 1 - e^{-\kappa t}.$$

$T_2$ similarly is an independent random variable, where $T_2 - T_1$ also is governed by $S$, whence

$$\Pr(T_2 \le T_1 + t) = 1 - S(t) = 1 - e^{-\kappa t},$$

and likewise for the remaining $T_i$. This demonstrates that $n$ is exactly as described in the question: it is the largest number of "exponential summands" that fit within the interval $(0, b]$. It is a realization of the random variable

$$N(0,b) = \max\{i\,|\,T_1+T_2+\cdots+T_i \le b\}.$$

Poisson distribution

Let $k\ge 0$ be an integer. What is the chance, $p_k$, that there are exactly $k$ points in the interval $[0,1]$? I am going to deduce the answer from the ergodic property of the Poisson process, which I take to be intuitive: because the process within the unit interval $[0,1]$ is the same as (and independent of) the process within any unit interval $[t-1, t]$, we may deduce properties of the process by varying $t$ for a single realization and studying the point patterns that show up. In particular, $p_k$ must equal the limiting proportion of time that the number of points $N(t)=N(t-1,t)$ in the interval $[t-1,t]$ equals $k$. We may express this formally using the indicator function $I$ which is equal to $1$ when its argument is true and is $0$ otherwise:

$$p_k = \lim_{x\to\infty}\frac{1}{x-1}\int_1^x I(N(t)=k) dt.$$

The integral is the total duration between $1$ and $x$ when interval $[t-1,t]$ contains exactly $k$ points, while the denominator of the fraction of course is the total elapsed time from $1$ to $x$.

This is a plot of $N(t)$ for one realization of a Poisson process with rate $\lambda=2.5$. White horizontal lines mark the values $k=1,2,3$ on the vertical axis, which extends from $0$ to $4$. The gray vertical lines show the times at which points occur in this realization. The dotted blue vertical lines show the same points shifted one unit to the right. The solid red curve plots $N(t)$ beginning at $t=1$. (It extends infinitely far to the right but not all of it could be drawn!)

$N(1)=4$ because four points (gray lines) appear in the first unit of time, $[0,1]$. The plot of $N(t)$ then rises by one unit every time $t$ another gray line is encountered moving left to right, because this is when the interval $[t-1,t]$ picks up that point. It falls by one unit every time $t$ a blue line is encountered, because this is the same thing as losing the value $t-1$ from the interval $[t-1,t]$.

The proportion of time spent at each height $k$ estimates $p_k$, the chance that any unit interval contains exactly $k$ points.

Suppose the interval $[t-1, t]$ contains $k$ points. As we slide it to the right, let's keep track of the new points it picks up and the old points that drop off. Two simple relations determine everything:

1. Between $t$ and $t+dt$ (for $dt \ge 0$), the expected number of new points is $\lambda\, dt$. (That's homogeneity.)

2. However, the expected number of points that are lost is $k\, dt$ because there are $k$ points randomly located within the interval. (This, too, is due to homogeneity.)

Almost surely at most one point is added or lost at any instant. (If not, there would be a positive lower bound to the chance of two or more points appearing in arbitrarily small intervals $[t, t+dt]$, but since the expected number of points in such intervals is only $\lambda dt$--which grows vanishingly small with small $dt$--this is impossible.) Accordingly, there are only two transitions that have nonzero chance of occurring: from $k$ points to $k+1$ points and from $k$ points to $k-1$ points. Their instantaneous rates are

$$\tau(k\to k+1) =\lim_{dt\to 0^{+}} \frac{\lambda\, dt}{dt} = \lambda$$

and, if $k \ge 1$,

$$\tau(k\to k-1) = \lim_{dt\to 0^{+}} \frac{k\, dt}{dt} = k.$$

This might look complicated, since it establishes a system of differential equations for the infinitely many probabilities $p_k$. However, because the process is homogeneous these probabilities are unchanging.

Look first at the case $k=0$. The expected rate at which $p_0$ changes (namely, zero) is the expected rate of transitions $1\to 0$ from states with $k=1$ minus the expected rate of transitions $0\to 1$ to states with $k=1$. Thus, by virtue of the simple relations (1) and (2),

$$(1)p_1 - (\lambda)p_0 = 0.$$

This enables us to find $p_1$ in terms of $p_0$:

$$p_1 = \lambda p_0.\tag{3}$$

Now consider the general situation $k\gt 0$. There are four transitions potentially changing $p_k$, all of which must balance out in expectation: $k\to k+1$ and $k\to k-1$ make it decrease, while $k+1\to k$ and $k-1\to k$ increase it. Thus, again using the simple relations (1) and (2) to compute the instantaneous rates,

$$[(\lambda)p_{k-1} - (k) p_k] + [(k+1)p_{k+1} - (\lambda)p_k] = 0.$$

We may assume inductively that the first term in brackets balances out (something we just showed for the case $k=0$), thereby automatically balancing out the second term in brackets and easily giving the solution

$$p_{k+1} = \frac{\lambda}{k+1}p_k.\tag{4}$$

Formulas $(3)$ and $(4)$ determine all the $p_k$ in terms of $p_0$: the solution is

$$p_k = p_0\frac{\lambda^k}{k!}.\tag{5}$$

(Proof: this formula satisfies the initial condition $(3)$ as well as the recursion $(4)$.)

The connection between the Exponential and Poisson parameters

There are two ways to find $p_0$. First, we may exploit what we learned previously about exponential waiting times: $p_0$ is the chance that an Exponential variable of parameter $\kappa$ exceeds $1$. This is

$$p_0 = e^{-\kappa}.\tag{6}$$

The second is the fact that the probabilities sum to unity, whence

$$1 = \sum_{k=0}^\infty p_0\frac{\lambda^k}{k!} = p_0\sum_{k=0}^\infty \frac{\lambda^k}{k!} = p_0 e^{\lambda}.$$

Therefore

$$p_0 = e^{-\lambda}\tag{7}$$

is the unique value that works. Consequently $(5)$ is now fully explicit:

$$p_k = e^{-\lambda} \frac{\lambda^k}{k!}.$$

This is the Poisson distribution.

Equating $(6)$ and $(7)$ reveals that

$$\kappa = \lambda.$$

This explicitly relates the parameter of the Exponential waiting times to the parameter of the Poisson distribution.

Simulation for greater understanding

The first figure did not show a sufficiently long time span to estimate the $p_k$ accurately. Consider, though, a portion another realization of $N(t)$ that is a hundred times longer:

(The vertical lines are no longer drawn, because they are so numerous.)

Here are the proportions of time spent for each $k$. Beneath them are the proportions for the Poisson$(2.5)$ distribution:

         0      1      2      3      4      5      6      7      8     9
y   0.0745 0.2068 0.2637 0.2215 0.1290 0.0660 0.0235 0.0128 0.0016 6e-04
fit 0.0821 0.2052 0.2565 0.2138 0.1336 0.0668 0.0278 0.0099 0.0031 9e-04

The agreement is apparent--although still imperfect, because this is still just a short initial segment of the realization.

---

Here is the R code used to produce the figures. Experiment with lambda and n to get a feel for this analysis.

lambda <- 2.5               # Poisson intensity
n <- 1000                   # Number of points to realize
x <- cumsum(rexp(n, lambda))# Accumulate the waiting times

# Compute the proportion of times for each `k` and compare to the Poisson distribution.
f <- ecdf(x)                # The ECDF does the work of computing N(t)
b <- max(x)
u <- seq(1, b, length.out=10*n)
y <- table(round(n*(f(u) - f(u-1)), 4))
y <- y / sum(y)
fit <- dpois(as.numeric(names(y)), lambda)
round(rbind(y, fit), 4)

# Plot N(t)
y.max <- max(as.numeric(names(y)))
curve(ifelse(x >= 1 & x <= b, n*(f(x) - f(x-1)), NA), 0,b, ylim=c(0, y.max*1.01),
      n=max(10001, 10*n), xlab="t", ylab="", col="#00000080",
      yaxp=c(0, y.max, y.max), bty="n", yaxt="n", yaxs="i")
rect(0, 0, b, y.max, col="#f4f4f4", border=NA)
abline(h=0:y.max, col="White")
if (n < 1000) {
  abline(v=x, lty=1, col="Gray")
  abline(v=x[x <= b-1]+1, lty=3, col="Blue")
}
curve(ifelse(x >= 1 & x <= b, n*(f(x) - f(x-1)), NA), add=TRUE,
      n=max(10001, 10*n), lwd=2, col="#802020")