3. The Poisson Counting Process

Contents

3. The Poisson Counting Process#

A Poisson counting process \(N(t)\) is a continuous time stochastic process that assumes values on the nonnegative integers, and which can be defined in the following ways. Let \(p_n(t) = \text{Prob}\, \{N(t) = n\},\ t \geq 0,\ n = 0,\ 1,\ 2,\ \ldots\). We define \(N(t)\) by a collection of differential equations for \(p_n(t)\):

(13)#\[\begin{split}\begin{aligned} \frac{d}{dt}\ p_0(t) &= -\lambda p_0(t), p_0(0) = 1 \\ \frac{d}{dt}\ p_1(t) &= -\lambda p_1(t) + \lambda p_0(t) \\ \frac{d}{dt}\ p_j(t) &= -\lambda p_j(t) + \lambda p_{j-1}(t),\ j \geq 1 \\ \end{aligned}\end{split}\]

where \(\lambda > 0\), and \(p_j(0) = 0\) for \(j \geq 1\). These equations can be arranged conveniently in the form of the infinite matrix equation.

\[\begin{split} \frac{d}{dt}\ \begin{pmatrix} p_0(t) \\ p_1(t) \\ p_2(t) \\ \vdots \end{pmatrix}\ =\ \begin{pmatrix} -\lambda & 0 & 0 & \ldots \\ \lambda & -\lambda & 0 & \ldots \\ 0 & \lambda & -\lambda & \ldots \end{pmatrix}\ \begin{pmatrix} p_0(t) \\ p_1(t) \\ p_2(t) \\ \vdots \end{pmatrix} \end{split}\]

The positive parameter \(\lambda\) is called the rate of the process. The solution of this recursive system of differential equations is readily found to be

\[\begin{split} \begin{aligned} p_0(t) &= e^{-\lambda t} \\ p_1(t) &= \lambda t e^{-\lambda t} \\ \ \vdots \\ p_j(t) &= \frac{(\lambda t)^j}{j!}\ e^{-\lambda t},\, j \geq 1 \\ \end{aligned} \end{split}\]

Evidently, we could define the Poisson counting process by simply starting with the definitions of probabilities

\[ p_n(t) = \text{Prob } \{N(t) = n\} = \frac{(\lambda t)^n}{n!}\ e^{-\lambda t},\ n = 0,\ 1,\ \ldots \]

The mean of the Poisson counting process is readily calculated as

\[\begin{split} \begin{aligned} EN(t) &= \sum_{n=0}^{\infty}\ n p_n(t) \\ &= \sum_{n=0}^{\infty}\ n\ \frac{(\lambda t)^n}{n!}\ e^{-\lambda t} = \lambda t\ \sum_{n=0}^{\infty}\ \frac{(\lambda t)^{n-1}}{(n-1)!}\ e^{-\lambda t} \\ \end{aligned} \end{split}\]

or

\[ EN(t) = \lambda t. \]

This formula motivates the interpretation of \(\lambda\) as the rate or arrival rate of the process. Similar calculations show that

\[ EN(t)^2 = \lambda^2\, t^2 + \lambda t, \]

which implies that

\[ E\, \big(N(t) - EN(t)\big)^2 = \lambda t. \]

An example of a realization of a Poisson counting process is given in Fig. 1. The random points in time \(t_1,\ t_2,\ \ldots\) at which \(N(t)\) jumps are called arrival times.

../_images/fig-3-1_poisson_sample_path.png

Fig. 1 Figure 1. A realization (sample path) of a Poisson counting process \(N(t)\) with rate \(\lambda\). The process starts at \(N(0) = 0\) and increases by one unit at each random arrival time \(t_1, t_2, \ldots\). It is right-continuous: at an arrival time \(t_i\) it takes the post-jump value \(N(t_i)\) (filled dot), while the pre-jump left limit \(N(t_i^-)\) is shown as an open dot. The interarrival times \(t_i - t_{i-1}\) are independent and exponentially distributed with mean \(1/\lambda\).#

Notice that the event that the first arrival time \(t_1 \leq T_1\) is the event that \(N(T_1) \geq 1\), which is the complement of the event that \(N(T_1) = 0\). Therefore, we have that

\[ \text{Prob } \{t_1 \leq T_1\} = 1 - e^{-\lambda T_1}, \]

which states that the first arrival time is exponentially distributed.

From the above, it follows that a Poisson process \(N(t)\) is characterized by the properties

(i) \(N(0) = 0\)

(ii) For \(t_1 > t_2\),

\[\begin{split} \begin{aligned} \text{Prob } &\big\{N(t_1) - N(t_2) = k \big\} =\ \frac{(\lambda (t_1 - t_2))^k}{k!}\ e^{-\lambda (t_1 - t_2)}, \\ k &= 0,\ 1,\ \ldots \\ \end{aligned} \end{split}\]

(iii) If \(t_1 > t_2 > t_3 > t_4\), then \((N(t_1) - N(t_2))\) is statistically independent of \(N(t_3) - N(t_4)\), and \(E(N(t_1) - N(t_2))\, (N(t_3) - N(t_4)) = \lambda^2 (t_1 - t_2)\, (t_3 - t_4)\), so that the expectation of the products of the increments \((N(t_1) - N(t_2))\) and \((N(t_3) - N(t_4))\) is the product of their expectations.

Property (iii) is the “independent increments” property of a Poisson process.

For \(t_1 > t_2\), we have

\[\begin{split} \begin{aligned} E\, \big\{N(t_1) - N(t_2)\big\} &= \lambda (t_1 - t_2) \\ E\, \big\{N(t_1) - N(t_2)\big\}^2 &= \lambda^2 (t_1 - t_2)^2 + \lambda (t_1 - t_2). \\ \end{aligned} \end{split}\]

Assuming again that \(t_1 > t_2 > t_3 > t_4\), and using \(N(t_1) - N(t_3) = (N(t_1) - N(t_2)) + (N(t_2) - N(t_3))\), and \(N(t_2) - N(t_4) = (N(t_2) - N(t_3)) + (N(t_3) - N(t_4))\), we obtain

\[ E\, \big\{ (N(t_1) - N(t_3))\, (N(t_2) - N(t_4)) \big\} = \lambda^2 (t_1 - t_3)\ (t_2 - t_4) + \lambda (t_2 - t_3). \]

It follows that

\[ \text{cov}\ \big( N(t_1) - N(t_3),\ N(t_2) - N(t_4) \big) = \lambda (t_2 - t_3), \]

so that the covariance between two increments is proportional to their overlap in time \((t_2 - t_3)\).

Using the preceding formulas, we can calculate the autocorrelation function \(R(t_1,\, t_2)\) and the autocovariance function \(C(t_1,\, t_2)\) for the Poisson process. We find

\[\begin{split} R(t_1,\, t_2) = \begin{cases} \lambda t_2 &+ \lambda^2 t_1\, t_2,\ t_1 \geq t_2 \\ \lambda t_1 &+ \lambda^2 t_1\, t_2,\ t_1 \leq t_2 \end{cases} \end{split}\]

or

(14)#\[ R(t_1,\, t_2) = \lambda\, \min\, (t_1,\, t_2) + \lambda^2 t_1\, t_2. \]

Note that \(\min\, (t_1,\, t_2)\) is the length of overlap of \(N(t_2)\) and \(N(t_1)\). It follows from (14) that

(15)#\[ C\, (t_1,\, t_2) = \lambda\, \min\, (t_1,\, t_2). \]

From (15) it follows that \(N(t)\) is mean square continuous. (Why?) From (14), we also have that

\[\begin{split} \frac{\partial R(t_1,\, t_2)}{\partial t_1}\ =\ \begin{cases} 0 &+ \lambda^2 t_2,\ t_1 \geq t_2 \\ \lambda &+ \lambda^2 t_2,\ t_2 \geq t_1 \end{cases} \end{split}\]

or

(16)#\[ \frac{\partial R(t_1,\, t_2)}{\partial t_1}\ =\ \lambda u(t_2 - t_1) + \lambda^2 t_2 \]

where \(u(t)\) is the Heaviside unit step function defined by \(u(t) = 1\) for \(t \geq 0,\ u(t) = 0\) for \(t < 0\). It follows from (16) that

(17)#\[ \frac{\partial^2 R(t_1,\, t_2)}{\partial t_1 \partial t_2}\ =\ \lambda \delta(t_2 - t_1) + \lambda^2 \]

where \(\delta(\ \ )\) is the Dirac delta generalized function defined by

\[ \int_{-\infty}^{\infty} g(\tau)\, \delta(\tau)\, d\tau = g(0) \]

for all test functions \(g(t)\) that are continuous at \(t = 0\), and that go to zero sufficiently quickly as \(|t| \to \infty\).

Equation (17) states that \(\partial^2 R(t_1,\, t_2)/\partial t_1 \partial t_2\) does not exist as an ordinary function, so that \(R(t_1,\, t_2)\) is not twice differentiable. It follows that \(N(t)\) is not mean square differentiable, so that \(dN(t)/dt\) does not exist as an ordinary stochastic process.

Although \(dN(t)/dt\) does not exist in the mean square sense, it does exist as a “generalized stochastic process.” The sample paths of the Poisson process can be represented as

\[ N(t) = \sum_{i=1}^{\infty}\, u(t - t_i) \]

where \(\{ t_1,\, t_2,\, \ldots\}\) are the Poisson arrival times and \(u(t)\) is again the Heaviside unit step function. Taking the generalized derivative of the sample path \(N(t)\) and using the definition \(\frac{d}{dt}\, u(t - t_i) = \delta(t - t_i)\), we obtain

\[ z(t) \equiv\ \frac{dN(t)}{dt} = \sum_{i=1}^{\infty}\ \delta(t - t_i), \]

so that \(dN(t)/dt\) is a sum of \(\delta\)-impulses at the random arrival times \(t_i\). From our preceding results, it follows that

\[ Ez(t) = \lambda \]
\[ Ez(t_1) z(t_2) = \lambda^2 + \lambda \delta(t_1 - t_2) \]

The zero mean generalized stochastic process \(z(t) - \lambda\), which has autocorrelation function \(\lambda \delta(t_1 - t_2)\) is an example of a white noise. Any stochastic process with an autocorrelation function proportional to \(\delta(t_1 - t_2)\) is known as a white noise. White noise is the elementary building block from which the processes of this book are assembled: it drives the stochastic differential equations of 5. Stochastic Processes Driven by a Poisson Counting Process, reappears in the limit construction of the Wiener process in 6. The Wiener Process, and is the fundamental input \(w(t)\) of the Wold representation in 8. Spectral Densities.

The process \(z(t) = dN(t)/dt\) does not exist as an ordinary stochastic process, but can be regarded as a particular kind of limit point of a process that does exist as an ordinary stochastic process. In particular, for \(\epsilon > 0\), define

\[ y(t)\ =\ \frac{N(t + \epsilon) - N(t)}{\epsilon} \]

where \(N(t)\) is a Poisson counter with rate \(\lambda\). It follows that

\[ y(t) = k/\epsilon \]

where \(k = N(t + \epsilon) - N(t) =\) number of arrivals of the \(N(t)\) process in the interval \((t,\, t + \epsilon)\). Therefore,

\[ \text{Prob } \big\{y(t) = k/\epsilon\big\}\ =\ \frac{e^{-\lambda \epsilon} (\lambda \epsilon)^k}{k!} \]

Using this Poisson probability distribution, and the above results on moments of a Poisson process, we can deduce that

\[ Ey(t) = \lambda \]
\[\begin{split} R(t_1,\, t_2) = \begin{cases} \lambda^2 & |t_1 - t_2| > \epsilon \\ \lambda^2 &+ \frac{\lambda}{\epsilon}\ -\ \frac{\lambda |t_1 - t_2|}{\epsilon^2}\, ,\ \text{ for }\ |t_1 - t_2| \leq \epsilon \end{cases} \end{split}\]

The function \(\lambda/\epsilon - \lambda |t_1 - t_2|/\epsilon^2\) for \(|t_1 - t_2| \leq \epsilon\) is plotted in Fig. 2, and inscribes a triangle of area \(\lambda\) above the horizontal axis. It is known that the limit as \(\epsilon \to 0\), of the ordinary function \(\max\ (0,\, \frac{\lambda}{\epsilon}\ -\ \frac{\lambda |t_1 - t_2|}{\epsilon^2})\) defines the Dirac delta generalized function with mass \(\lambda\), \(\lambda \delta(t_1 - t_2)\).

../_images/fig-3-2_delta_triangle.png

Fig. 2 Figure 2. The triangular function \(\lambda/\epsilon - \lambda|t_1 - t_2|/\epsilon^2\) on \(|t_1 - t_2| \leq \epsilon\), shown for three window widths \(\epsilon\). Each triangle has height \(\lambda/\epsilon\), base \([-\epsilon, \epsilon]\), and area exactly \(\lambda\). As \(\epsilon \to 0\) the triangles become taller and narrower while preserving their area, so in the limit they define the Dirac delta generalized function with mass \(\lambda\), namely \(\lambda\,\delta(t_1 - t_2)\).#

It is a “spike” of “mass” \(\lambda\) at \(t_1 - t_2 = 0\), and is equal to zero for \(t_1 - t_2 \neq 0\).

As an example of one use of the generalized stochastic process \(z(t) = dN(t)/dt\), let \(L(\tau),\ \tau \in [0,\, \infty)\) be a continuous and square integrable function, i.e., \(\int_0^{\infty} L(\tau)^2\, d\tau < +\infty\). Then consider a stochastic process defined by the distributed lag

\[ Y(t) = \int_0^{\infty} L(\tau) z(t - \tau)\, d\tau. \]

Substituting \(z(t) = \sum_{i=1}^{\infty}\, \delta(t - t_i)\) and using the definition of the delta generalized function, we find

(18)#\[\begin{split} \begin{aligned} Y(t) &= \int_0^{\infty} L(\tau) \sum_{i=1}^{\infty} \delta(t - t_i - \tau)\, d\tau \\ &= \sum_{i=1}^{\infty} \int_0^{\infty} L(\tau) \delta(t - t_i - \tau)\, d\tau \\ Y(t) &= \sum_{i=1}^{\infty} L(t - t_i),\ t \geq 0 \\ \end{aligned} \end{split}\]

The process \(Y(t)\) defined by (18) is an ordinary stochastic process, consisting of a sum of the function \(L\) shifted by the random arrival times \(t_i\). Such a process \(Y(t)\) is called shot noise. In 8. Spectral Densities we shall show how to calculate its first and second order moments. Shot noise is the prototype of the lesson made precise in 4. The Concept of “Physical Realizability”: the generalized white noise \(dN/dt\) lives only under the integral sign, while the integrated process \(Y(t)\) is perfectly ordinary.

Exercises#

The following exercises use simulation to illustrate the properties of the Poisson counting process derived above. They use only numpy and matplotlib.

import numpy as np
import matplotlib.pyplot as plt

The key fact we exploit is that a Poisson process of rate \(\lambda\) can be simulated by drawing independent exponential interarrival times with mean \(1/\lambda\) and accumulating them: the cumulative sums are the arrival times \(t_1, t_2, \ldots\), and \(N(t)\) counts how many arrival times do not exceed \(t\).

def poisson_arrivals(lam, T, rng):
    """Arrival times of a rate-lam Poisson process on the interval [0, T]."""
    times, t = [], 0.0
    while True:
        t += rng.exponential(1.0 / lam)
        if t > T:
            break
        times.append(t)
    return np.array(times)

Exercise 1

Take a Poisson process with rate \(\lambda = 2\).

(a) Simulate and plot three sample paths of \(N(t)\) on \([0, 10]\), together with the mean function \(E\,N(t) = \lambda t\).

(b) By simulating many independent paths, verify the moment formulas \(E\,N(T) = \lambda T\) and \(\operatorname{Var} N(T) = \lambda T\) derived in the text. (Recall that \(\operatorname{Var} N(T) = E\,N(T)\) is a signature property of the Poisson distribution.)

Exercise 2

The text shows that the first arrival time is exponentially distributed, \(\operatorname{Prob}\{t_1 \le T_1\} = 1 - e^{-\lambda T_1}\). More generally, the interarrival times \(t_{i} - t_{i-1}\) of a Poisson process are independent and exponentially distributed with mean \(1/\lambda\).

Simulate one long sample path (rate \(\lambda = 2\) on \([0, 5000]\)), extract its interarrival times, and compare their histogram with the exponential density \(\lambda e^{-\lambda u}\). Report the sample mean interarrival time and compare it to \(1/\lambda\).

Contents