In [10]:
import sympy
import numpy as np
import matplotlib.pyplot as plt
import pynverse
%matplotlib inline
import pynverse
First, let’s compute $\Lambda(t)=\int_0^t \lambda(t) dt$
$\lambda(t)=-3t(t-5)+20$
Then we try to solve the inverse of $\Lambda(t)$
In [11]:
t=sympy.Symbol(“t”)
lm_t=-3*t*(t-5)+20
Lm_t=sympy.integrate(lm_t,(t,0,t))
#a to denote the arrival times from the homongeneous poisson with lam=1
a=sympy.Symbol(“a”)
sympy.N(sympy.solve(Lm_t-a,t)[2].subs({a:1}))
Out[11]:
$\displaystyle -2.12645858058961 + 2.0 \cdot 10^{-22} i$
In [12]:
#Use sympy.lambdify(args,expr) to convert
#a sympy funtion to a numeric function
Lm_t_func=sympy.lambdify((t),Lm_t)
In [13]:
def poisson_p():
arrival_time=np.array([])
t=0
T=float(Lm_t_func(6))
lmbda=1
#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 ):
arrival_time=np.append(arrival_time,t)
inter_a=-1/lmbda*np.log(np.random.rand())
t=t+inter_a
return arrival_time
homo_arrivals=poisson_p()
In [14]:
#implement inversion method
pynverse.inversefunc(Lm_t_func, homo_arrivals, [0,6] )
Out[14]:
array([0.04721955, 0.04870638, 0.0524114 , 0.11448708, 0.11952843,
0.17901235, 0.19808311, 0.2966206 , 0.31594621, 0.43454996,
0.45637323, 0.60895744, 0.65416302, 0.69505227, 0.81304122,
0.86698297, 0.8906286 , 0.94046899, 0.98647636, 1.13141834,
1.20807892, 1.21282713, 1.22670543, 1.26065658, 1.31216987,
1.32549402, 1.32568909, 1.3487755 , 1.44327526, 1.47187343,
1.53973211, 1.54251816, 1.56163448, 1.64914955, 1.73095213,
1.81631209, 1.81776031, 1.92560919, 1.93294536, 1.94849207,
1.99413433, 2.00982058, 2.08570291, 2.09394846, 2.11734499,
2.13493206, 2.17491529, 2.18437223, 2.20350998, 2.22396921,
2.311048 , 2.38803323, 2.3891933 , 2.43127914, 2.47032054,
2.50580318, 2.52336242, 2.55629706, 2.57341996, 2.57502774,
2.59431349, 2.60266776, 2.62164491, 2.63371295, 2.65509883,
2.76817327, 2.81326454, 2.81647597, 2.82196582, 2.84213484,
2.88904134, 2.92485642, 2.96156378, 2.99734254, 3.08261088,
3.09137355, 3.10697254, 3.13138636, 3.13533675, 3.15941951,
3.1612955 , 3.2131429 , 3.34587889, 3.36835294, 3.41415134,
3.48097327, 3.53636687, 3.60410213, 3.61715265, 3.64722625,
3.72073209, 3.76253422, 3.78394516, 3.91662238, 3.95231083,
3.95937412, 3.96043834, 3.96845012, 3.96955517, 3.98440627,
3.98804518, 4.04928466, 4.05518817, 4.10160164, 4.10342413,
4.13200953, 4.20779243, 4.22891935, 4.27509048, 4.27742102,
4.32143033, 4.38134013, 4.38620756, 4.40898299, 4.51879914,
4.55956404, 4.56795031, 4.57572577, 4.67996225, 4.90616888,
4.90930161, 5.02505229, 5.02853764, 5.27260634, 5.35794745,
5.36862851, 5.42946868, 5.53633561, 5.56804037, 5.70790834,
5.82616835])
In [15]:
In [16]:
In [17]:
In [18]:
In [19]: