In [0]:
import numpy as np
import scipy.stats as spst
import matplotlib.pyplot as plt
%matplotlib inline
How to generate the samples using rejection sampling method?
Collect 1000 samples from the following distribution:
segment 1: P=0.4, Poisson(1)
segment 2: Poisson(10)
In [0]:
p=0.4
lm1=1
lm2=10
theta=0.08
x=np.arange(100)
PMF_t=spst.poisson.pmf(x,lm1)*p+spst.poisson.pmf(x,lm2)*(1-p)
PMF_p=(1-theta)**x*theta
In [0]:
plt.bar(x,PMF_t,label=”target distribution”)
plt.scatter(x,PMF_p,label=”proposal distribution”)
plt.legend()
plt.show()
In [0]:
#compute c value
In [0]:
def pmf_p(x):
return (1-theta)**x*theta
def pmf_t(x):
return spst.poisson.pmf(x,lm1)*p+spst.poisson.pmf(x,lm2)*(1-p)
#rejection sampling to get a sample
def rejection():
return proposal
samples=[rejection() for i in range(1000)]
In [0]:
samples=np.array(samples)
values,counts=np.unique(samples,return_counts=True)
x=np.arange(np.max(samples))
pmf_1=spst.poisson.pmf(x,1)
pmf_2=spst.poisson.pmf(x,10)
p=0.4
pmf=pmf_1*p+pmf_2*(1-p)
plt.bar(x,pmf,label=”theoretical”)
plt.scatter(values,counts/np.sum(counts),zorder=2, label=”samples”)
plt.legend()
plt.show()
In [0]: