Simulating a Poisson Distribution
Simulating a Poisson distribution (starting with a uniform RV)
This is a simple example of simulating a Poisson distribution.
import numpy as np
from matplotlib import pyplot as plt
Recall that a Poisson distribution is defined as having exponential interarrival times. This is useful in simulating a Poisson distribution, as we can easily get an exponential RV from uniforms with inversion.
def gen_exp(size=1):
uniforms = np.random.uniform(0,1,size=size)
return -np.log(uniforms)
plt.hist(gen_exp(1000), bins=30)
Now that we have our exponential RV we can use it to generate interarrival times for our Poisson RV.
def poisson(time):
arrivals = 0
total = 0
while total < time:
total += gen_exp()
arrivals += 1
return arrivals - 1
The below histograms illustrate that our Poisson RV does indeed behave as expected, given various expected numbers of arrivals \lambda
# expected Poisson values, for use with overlaying onto histograms
def pmf(lam, k):
return (lam**k * np.exp(-lam))/np.math.factorial(k)
sample = [poisson(50) for _ in range(10_000)]
exp_space = np.linspace(25, 80, 100)
expected = pmf(exp_space, 50)
plt.plot(exp_space, expected)
plt.title(r'$\lambda = {}$'.format(50))
plt.hist(sample, density=True, bins=30)
sample = [poisson(10) for _ in range(10_000)]
exp_space = np.linspace(0, 25, 100)
expected = pmf(exp_space, 10)
plt.plot(exp_space, expected)
plt.title(r'$\lambda = {}$'.format(10))
plt.hist(sample, density=True, bins=20)
sample = [poisson(5) for _ in range(10_000)]
exp_space = np.linspace(0, 16, 100)
expected = pmf(exp_space, 5)
plt.plot(exp_space, expected)
plt.title(r'$\lambda = {}$'.format(5))
plt.hist(sample, density=True, bins=15)
sample = [poisson(1) for _ in range(10_000)]
exp_space = np.linspace(0, 7, 100)
expected = pmf(exp_space, 1)
plt.plot(exp_space, expected)
plt.title(r'$\lambda = {}$'.format(1))
plt.hist(sample, density=True, bins=8)
