In [17]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import scipy.special as spsp
In [18]:
def start(lm):
start_x=int(lm)
#compute the starting point outside to improve efficiency
pmf=np.exp(-lm)
cdf=pmf
for x in range(1, start_x+1):
pmf=pmf*lm/x
cdf=cdf+pmf
return cdf,pmf
def poisson_sampling(lm,cdf_start, pmf_start):
#starting value of search
x=int(lm)
# cdf corresponding to the starting value?
cdf=cdf_start
pmf=pmf_start
# compare u with cdf
u=np.random.rand()
# if u>= cdf: upward
if u>=cdf:
x=x+1
pmf=pmf*lm/x
cdf=cdf+pmf
while u>=cdf:
x=x+1
pmf=pmf*lm/x
cdf=cdf+pmf
return x
#else downward
else:
cdf=cdf-pmf
while u