import matplotlib.pyplot as plt
import numpy as np
import scipy.stats
np.random.seed(42)
Generate Poisson process data and some rate and also at half that rate.
def poisson_process_wait_times(rate=1.0, num_events=100000):
"""Generate waiting times between events in a Poisson process."""
scale = 1.0 / rate
waiting_times = np.random.exponential(scale=scale, size=num_events)
return waiting_times
base_rate = 2.0
half_rate = base_rate / 2.0
wait_times = poisson_process_wait_times(rate=base_rate)
wait_times_half = poisson_process_wait_times(rate=half_rate)
Sanity check: Visualize the first 100 event times from each process.
fig, ax = plt.subplots(figsize=(12, 1.5))
ax.scatter(wait_times[:100].cumsum(), [1]*100, s=1, color='blue', label='Rate {0}'.format(base_rate))
ax.scatter(wait_times_half[:100].cumsum(), [0]*100, s=1, color='red', label='Rate {0}'.format(half_rate))
ax.set_xlim(0, 110)
ax.tick_params(axis='y', left='off', right='off', labelleft='off')
ax.set_xlabel('Time')
ax.set_title('First 100 event times from each process')
ax.legend(loc='upper right')
<matplotlib.legend.Legend at 0xaec576c>
Plot the pdf's of the waiting times between (a) every other Poisson event, and (b) every Poisson event at half the rate.
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))
wait_times_every_other = wait_times[0::2] + wait_times[1::2]
hist_kwargs = {'normed': True, 'bins': 100, 'range': (0, 6), 'alpha': 0.5}
ax1.hist(wait_times_every_other, label='Wait times (every other)', **hist_kwargs)
ax2.hist(wait_times_half, label='Wait times (half rate)', **hist_kwargs)
x = np.linspace(0, 6, num=1000)
plot_kwargs = {'linewidth': 2, 'alpha': 0.8}
for ax in (ax1, ax2):
label = r'$Gamma(\alpha=2, \beta={0})$'.format(base_rate)
gamma_pdf = scipy.stats.gamma(2, scale=(1.0 / base_rate)).pdf
ax.plot(x, gamma_pdf(x), label=label, color='red', **plot_kwargs)
label = r'$Exponential(\lambda={0})$'.format(half_rate)
expon_pdf = scipy.stats.expon(scale=(1.0 / half_rate)).pdf
ax.plot(x, expon_pdf(x), label=label, color='green', **plot_kwargs)
ax.legend()