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\)).
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:
- Generate a random number \(r\) from the uniform distribution in [0, 1) (i.e., between 0 and 1, including 0, excluding 1).
- 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.
- Suppose we have \(\lambda = 0.1\) and \(\delta t = 0.1\). Thus, \(\lambda\delta t = 0.01\).
- We start from \(N(0) = 0\).
- 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\).
- 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\).
- 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\).
- 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\).)
Comments
Post a Comment