This Jupyter notebook is the Python equivalent of the R code in section 13.5 R, pp. 534 - 536, Introduction to Probability, 1st Edition, Blitzstein & Hwang.
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline
In Chapter 5, we discussed how to simulate a specified number of arrivals from a one-dimensional Poisson process by using the fact that the interarrival times are i.i.d. Exponentials. In this chapter, Story 13.2.3 tells us how to simulate a Poisson process within a specified interval $(0, \left.L\right]$. We first generate the number of arrivals $N(L)$, which is distributed $Pois(\lambda L)$. Conditional on $N(L) = n$, the arrival times are distributed as the order statistics of $n$ i.i.d. $Unif(0, L)$ r.v.s. The following code simulates arrivals from a Poisson process with rate 10 in the interval $(0, \left. 5\right]$:
np.random.seed(2971215073)
from scipy.stats import poisson, uniform
L = 5
lambd = 10
n = poisson.rvs(lambd*L, size=1)[0]
t = np.sort(uniform.rvs(loc=0, scale=(L-0), size=n))
To visualize the Poisson process we have generated, we can plot the cumulative number of arrivals $N(t)$ as a function of $t$:
plt.figure(figsize=(12, 4))
nt = np.arange(1, n+1)
plt.plot(t, nt, linestyle='-', drawstyle='steps')
plt.xlabel(r'$t$')
plt.xlim([0,5])
plt.ylabel(r'$N(t)$')
plt.title(r'Number of arrivals in Poisson process, rate=10, in interval (0,5]')
plt.show()
This produces a staircase plot as in Figure 13.8.
The following code starts with an array of arrival times t
and the corresponding number of arrivals n
, generated according to the previous section. For each arrival, we flip a coin with probability p
of Heads; these coin tosses are stored in the array y
. Finally, the arrivals for which the coin landed Heads are labeled as type-1; the rest are labeled as type-2. The resulting arrays of arrival times $t_1$ and $t_2$ are realizations of independent Poisson processes, by Theorem 13.2.13. To perform this with $p = 0.3$, we can execute
np.random.seed(1)
from scipy.stats import binom
p = 0.3
y = binom.rvs(1, p, size=n)
t1 = t[np.where(y==1)]
t2 = t[np.where(y==0)]
As before, we can plot the number of arrivals in each Poisson process, $N_{1}(t)$ and $N_{2}(t)$, as a function of $t$.
_, (ax1, ax2) = plt.subplots(1, 2, sharey=True, figsize=(12, 4))
# graph for t1
nt1 = np.arange(1, t1.size+1)
ax1.plot(t1, nt1, linestyle='-', drawstyle='steps', color='#ef8a62', label='rate=3')
ax1.set_xlabel(r'$t$')
ax1.set_xlim([0,5])
ax1.set_ylabel(r'$N_{1}(t)$')
ax1.legend()
# graph for t2
nt2 = np.arange(1, t2.size+1)
ax2.plot(t2, nt2, linestyle='-', drawstyle='steps', color='#67a9cf', label='rate=7')
ax2.set_xlabel(r'$t$')
ax2.set_xlim([0,5])
ax2.set_ylabel(r'$N_{2}(t)$')
ax2.legend()
plt.suptitle('Number of arrivals in two independent Poisson processes')
plt.show()
Simulating a 2D Poisson process is nearly as straightforward as simulating a 1D Poisson process. In the square $(0, \left. L\right] \times (0, \left. L \right]$, the number of arrivals is distributed $Pois(\lambda L^2)$. Conditional on the number of arrivals, the locations of the arrivals are Uniform over the square, so by Example 7.1.22, their horizontal and vertical coordinates can be generated independently:
np.random.seed(1)
L = 5
lambd = 10
n = poisson.rvs(lambd*(L**2), size=1)[0]
x = uniform.rvs(loc=0, scale=(L-0), size=n)
y = uniform.rvs(loc=0, scale=(L-0), size=n)
We can then easily plot the locations of the arrivals as in Figure 13.7:
plt.figure(figsize=(6, 6))
plt.plot(x, y, 'x')
plt.xlabel(r'$x$')
plt.ylabel(r'$y$')
plt.title(r'Arrivals in simulated $Pois(\lambda L^2)$, $\lambda={}$, $(0,L] \times (0,L]$'.format(lambd))
plt.show()
© Blitzstein, Joseph K.; Hwang, Jessica. Introduction to Probability (Chapman & Hall/CRC Texts in Statistical Science).