In [5]:
import numpy as np
In [6]:
def lmbda(t):
return np.sin(t/60) + 10
def f(t):
return 60 – 60 * np.cos(t/60) + 10 * t
def rejection_sampling(N):
res = []
for i in range(N):
x = np.random.rand()*480
while np.random.rand() >= lmbda(x) / 11:
x = np.random.rand() * 480
res.append(x)
return np.array(res)
def order_statistics(T):
# Can integrate lambda(t) via sympy or by hand.
Lmbda_T = f(T)
N = np.random.poisson(Lmbda_T)
arrivals = rejection_sampling(N)
return np.sort(arrivals)
arrivals = order_statistics(480)
inter_arrivals = np.diff(np.insert(arrivals, 0, 0))
print(np.mean(inter_arrivals))
print(len(inter_arrivals))
0.0988737112798129
4854