In [1]:
import numpy as np
import matplotlib.pyplot as plt
import sympy
%matplotlib inline
In [2]:
#Compute the expected number of arrivals between 0 and T
#Your function should return a float
#You can use float(result) to the sympy result into a float number
def Ex_arrivals(T):
t=sympy.Symbol(“t”)
lm_t=42*t*sympy.exp(-t**2)
Ex=sympy.integrate(lm_t,(t,0,T))
return float(Ex)
Ex_arrivals(T = 6)
Out[2]:
20.999999999999996
In [3]:
#simulate the actual number of arrivals between 0 and T
#hint: you can need to call Ex_arrivals(T) function you
#constructed in the last step
def N_arrivals(T):
return np.random.poisson(Ex_arrivals(T))
N_arrivals(T = 6)
Out[3]:
16
In [4]:
# compute the inverse of cdf of the arrival time for sampling the arrival time
# Your sympy result can be more than one possible inverse
# return only the valid sympy result in our setting (i.e should be positive > 0)
def sympy_invcdf(T):
t=sympy.Symbol(“t”)
lm_t=42*t*sympy.exp(-t**2)
Ex=sympy.integrate(lm_t,(t,0,T))
#pdf of the target distribution
pdf=lm_t/Ex
#computing the cdf of the target (use inverse transform method)
cdf=sympy.integrate(pdf,(t,0,t))
u=sympy.Symbol(“u”)
#checking the validity of the results
#print(sympy.N(sympy.solve(cdf-u,t)[0].subs({u:0.8})))
#print(sympy.N(sympy.solve(cdf-u,t)[1].subs({u:0.8})))
return sympy.solve(cdf-u,t)[1]
inv_cdf=sympy_invcdf(T = 6)
inv_cdf
Out[4]:
$\displaystyle \sqrt{\log{\left(\frac{1}{- u e^{36} + u + e^{36}} \right)} + 36}$
In [5]:
# Combine all the previous functions to simulate the process up to time T
# use sympy.lambdify function to convert the sympy function coming from sympy_invcdf()
# to a numpy function that can directly give arrival time
np.random.seed(1996)
def Simulate(T):
#decide the number of arrival samples needed
n=N_arrivals(T = 6)
#converting sympy expression to a sampling function
inv_cdf=sympy_invcdf(T = 6)
u=sympy.Symbol(“u”)
sample=sympy.lambdify((u),inv_cdf,”numpy”)
#generating n arrival times
atimes=sample(np.random.rand(n))
#sort the arrival times from smallest to largest
return np.sort(atimes)
result_simulation = Simulate(T = 6)
result_simulation
Out[5]:
array([0.07162409, 0.250669 , 0.26133756, 0.31571496, 0.40889667,
0.50773498, 0.50932239, 0.53008839, 0.60658651, 0.62853664,
0.66058952, 0.76430313, 0.88371247, 0.97831058, 1.06344273,
1.10980383, 1.18246613, 1.3000786 , 1.30693752, 1.53319418,
1.81310792, 1.85130138, 1.96527166])
In [7]:
In [0]: