In [2]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import sympy
In [3]:
#visualize the relationship between t and lambda
t=np.linspace(0,6,1000)
lm_t=-0.1*t*(t-5)+1
plt.plot(t,lm_t)
plt.show()

In [4]:
#find out the highest lambda between 0 and T
lmda_max=np.max(lm_t)
lmda_max
plt.plot(t,lm_t/lmda_max)
plt.show()

In [5]:
In [6]:
In [7]:
#generate the arrivals using thinning method
def poisson_np():
arrival_time=np.array([])
t=0
T=6
lmbda=lmda_max
#inter arrival time is a random sample from an exponential distribution
inter_a=-1/lmbda*np.log(np.random.rand())
t=t+inter_a
while (t<=T ):
#check if the proposal arrival gets accepted or not
if np.random.rand()<(-0.1*t*(t-5)+1)/lmbda:
arrival_time=np.append(arrival_time,t)
inter_a=-1/lmbda*np.log(np.random.rand())
t=t+inter_a
return arrival_time
poisson_np()
Out[7]:
array([0.59925791, 0.72935337, 1.05599228, 2.60365973, 3.43398877,
4.40890569, 4.4417671 ])
In [8]: