Simulating birth processes

 

The birth process $\{N(t)\}$ with birth rate $\lambda$ and initial population size $N(0)=n_0$ is characterized by the following differential-difference equations for the probability mass functions $p_n(t)=Pr(N(t)=n)$:

$$\begin{array}{}\text{(Eq:Birth0)}& \frac{d{p}_{n_0}(t)}{dt}& =& -\lambda n_0{p}_{n_0}(t),\text{(Eq:Birthn)}& \frac{d{p}_n(t)}{dt}& =& \lambda (n-1)p_{n-1}(t)-\lambda n{p}_n(t),(n\ge n_0+1)\end{array}$$

with the initial condition

$$p_n(0)={\delta }_{n,n_0}.$$

See also: Birth process

Here, we want to simulate this process numerically to obtain some concrete sample paths like this one:

To do so, we first discretize the time variable so that we consider time steps with some small interval $\delta t$: $N(0),N(\delta t),N(2\delta t),N(3\delta t),\cdots ,N(k\delta t),\cdots$. 

Next, we exploit the Markov property of the birth process: Every birth is the first birth since the last one. That is, given $N(t)=n_t$ at time $t$, we may regard this as the "initial" condition for the next step: $N(t+\delta t)=n_t$ or $n_t+1$. We can solve (Eq:Birth0) (with $n_0$ replaced with $n_t$) to find $p_{n_t}(t+\delta t)$ with the "initial" condition $p_n(t)={\delta }_{n,n_t}$. A bit of exercise gives

$$p_{n_t}(t+\delta t)=e^{-\lambda n_t\delta t},$$

which is the probability that $N(t+\delta t)=n_t$ (i.e., no birth during $\delta t$). Accordingly, we have the probability that $N(t+\delta t)=n_t+1$ (i.e., one birth during $\delta t$) as

$$p_{n_t+1}(t+\delta t)=1-e^{-\lambda n_t\delta t}.$$

Thus, the algorithm for simulating the birth process is the following:

1. Set $N(0)=n_0$; set k = 0.
2. Generate a uniformly distributed random number $r\in [0,1)$.
3. If $r<1-e^{-\lambda \delta tN(k\delta t)}$, then set $N((k+1)\delta t)=N(k\delta t)+1$; Otherwise, set $N((k+1)\delta t)=N(k\delta t)$.
4. Update $k:=k+1$, go to Step 2, and repeat.