Poisson process: Introduction

 We consider stochastic processes with continuous time. For example,

  • The number of births (or deaths) in a population by time \(t\).
  • The number of phone calls at an office by time \(t\).
  • The number of radioactive particles detected by a Geiger counter by time \(t\).
The Poisson process is a basic stochastic process that can be used to model above phenomena. Let's define it.

The Poisson process

Definition

Let \(N(t)\) be a time-dependent random variable of integer values, and \(\lambda > 0\) a constant. Here, \(t\) represents time. Suppose the probability mass function of \(N(t)\) is given by
\begin{equation} p_n(t) = \Pr(N(t) = n) = \frac{(\lambda t)^ne^{-\lambda t}}{n!}, ~~ (t \geq 0, n = 0, 1, 2,\cdots) \tag{Eq:PoissonProc}\end{equation}
A stochastic process \(\{N(t)\}\) characterized by the probability distribution (Eq:PoissonProc) is called a Poisson process.

Example: Let \(N(t)\) be the number of phone calls at an office by time \(t\). Then \(N(t)\) can be modeled as a Poisson process. (The same applies to other examples above.)

Why is it called the "Poisson" process? Because the probability distribution \(p_n(t)\) is the Poisson distribution with parameter (the mean of \(N(t)\)) \(\lambda t\). You should check that \begin{equation} \mu(t) = \mathbb{E}[N(t)] = \sum_{n=0}^{\infty}np_n(t) = \lambda t \tag{Eq:Mean}\end{equation} where \(\mathbb{E}[ \cdot ]\) represents the expectation value of the random variable in the argument. The variance is \begin{equation} \mathbb{V}[N(t)] = \mathbb{E}[(N(t) - \mu(t))^2] = \lambda t. \tag{Eq:Var}\end{equation}
What can we learn from this?
  • According to (Eq:Mean), the mean of \(N(t)\) is proportional to time \(t\). Thus, \(N(t)\) is a non-decreasing function of time, on average.
  • According to (Eq:Var), the variance of \(N(t)\) is also proportional to time \(t\). Thus, for a larger \(t\), the difference between different samples of \(N(t)\) becomes larger, on average.
Here is a sample path of a Poisson process (with \(\lambda = 0.1\)).



We can see step-wise increases of \(N(t)\) at random intervals.

Simulating Poisson processes

(See this page for details.)
How can we simulate a Poisson process with parameter \(\lambda\), such as in the above figure? Here, I explain a rough approximate method without details. In later posts, we will see why this method correctly (yet approximately) simulates the Poisson process and why it is a rough approximation. Since we cannot represent continuous variables on a computer, we must discretize them. In the present case, it is the time variable \(t\) that needs to be discretized. So, we consider discrete time steps:
\begin{equation}t = 0, \delta t, 2\delta t, 3\delta t, \cdots \end{equation} where \(\delta t\) is some small positive real number (floating-point number). As the initial condition, we assume that \begin{equation}N(0) = 0.\end{equation}
Now, suppose we have \(N(k\delta t) = n\) for some \(k = 0, 1, \cdots\). We determine the value of \(N((k+1)\delta t)\) in the following manner:
  1. Generate a random number \(r\) from the uniform distribution in [0, 1) (i.e., between 0 and 1, including 0, excluding 1).
  2. If \(r < \lambda \delta t\), then we set \(N((k+1)\delta t) = n + 1\), otherwise we set \(N((k+1)\delta t) = n\).
By this procedure, the value of \(N((k+1)\delta t)\) is incremented by 1 with probability \(\lambda t\), and stays the same (\(N((k+1)\delta t) = n\)) with probability \(1 - \lambda \delta t\). Let's see an example.
  1. Suppose we have \(\lambda = 0.1\) and \(\delta t = 0.1\). Thus, \(\lambda\delta t = 0.01\).
  2. We start from \(N(0) = 0\).
  3. Generate a random number from [0, 1). Say, we obtained \(r = 0.2358\). Since this is greater than \(\lambda\delta t = 0.01\), we set \(N(\delta t) = N(0) + 0 = 0\).
  4. Generate another random number. Say we obtained \(r = 0.0052103 < \lambda\delta t\). Then, we set \(N(2\delta t) = N(\delta t) + 1 = 1\).
  5. Generate another random number. Say we obtained \(r = 0.025139 > \lambda\delta t\). Then, we set \(N(3\delta t) = N(2\delta t) + 0 = 1\).
  6. We repeat this procedure over and over.

Since the quantity \(\lambda \delta t\) serves as a probability,  this method is applicable only when \(\lambda \delta t \leq 1\) (otherwise \(\lambda \delta t (> 1)\)cannot be interpreted as a probability). For example, if we have \(lambda = 20, \delta t = 0.1\), then \(\lambda \delta t = 2 > 1\) so that \(r < \lambda\delta t\) for any random number \(r \in [0, 1)\). This means \(N(k\delta t)\) is always incremented by 1 regardless of the random numbers, which is nonsense. To use a larger \(\lambda\) value, we replace \(r < \lambda \delta t\) with \[r < 1 - e^{-\lambda \delta t} (< 1)\] in Step 2. Note that when \(\lambda \delta t\) is sufficiently small, the Maclaurin expansion gives
\begin{equation} 1 - e^{-\lambda \delta t} = \lambda\delta t + o(\delta t^2). \end{equation}
That is, \(\lambda \delta t\) can be obtained as the first-order approximation of \(1 - e^{-\lambda \delta t}\). In a later post, we will see why \(1 - e^{-\lambda\delta t}\) is indeed the correct probability of increment. (Hint: \(1-e^{-\lambda \delta t}\) is the cumulative distribution function of the exponential distribution of \(\delta t\) with parameter \(\lambda\).)

Related videos:




Comments

Popular posts from this blog

Open sets and closed sets in \(\mathbb{R}^n\)

Euclidean spaces

Applications of multiple integrals