In [19]:
import numpy as np
import matplotlib.pyplot as plt
import math
import scipy.special as spsp
import scipy.stats as spst
%matplotlib inline
In [20]:
#coding for this cell is not required
import plotly.graph_objs as go
from ipywidgets import widgets
p_value=widgets.FloatSlider(
value=0.9,
min=0.01,
max=0.99,
step=0.01,
describe=”p”,
disabled=False,
continuous_update=False,
orientation=’horizontal’,
readout=True,
readout_format=’.2f’,
)
p=p_value.value
Outcomes_P=np.arange(10)
PMF_P=(1-p)**Outcomes_P*p
Outcomes_T=np.arange(10)
PMF_T=np.array([0.4,0.25,0.15,0.1,0.05,0.05,0,0,0,0])
c=np.max(PMF_T/PMF_P)
trace1=go.Bar(x=Outcomes_T,y=PMF_T,
hoverinfo=”text”,text=np.round(PMF_T,3),name=”Target distribution”+’ ‘*10 )
trace2= go.Bar(x=Outcomes_P,y=PMF_P,
hoverinfo=”text”,text=np.round(PMF_P,3),name=”Proposal distribution”+’ ‘*10 )
g = go.FigureWidget(data=[trace1,trace2],
layout=go.Layout(
title=dict(
text=’rejection sampling with p=’+str(p_value.value)+”and c=”+str(np.round(c,2)),
),
hovermode=None,
margin={‘l’: 0, ‘r’: 0, ‘t’: 0, ‘b’: 0},width=1000, height=400 ),
)
g.update_layout(barmode=’group’,
title_x=0.5,
title_y=0.9,
yaxis=dict(range=[0,1]),
legend=dict(
x=0.7,
y=0.7,
traceorder=”normal”,
font=dict(
family=”sans-serif”,
size=12,
color=”black”
))
)
def response1(change):
p=p_value.value
PMF_P=(1-p)** Outcomes_P*p
c=np.max(PMF_T/PMF_P)
with g.batch_update():
g.data[1].y=PMF_P
g.layout.title.text=’rejection sampling with p=’+str(p_value.value)+”and c=”+str(np.round(c,2))
g.data[1].text=np.round(PMF_P,3)
g.data[0].text=np.round(PMF_T,3)
p_value.observe(response1,names=”value”)
Widget1=widgets.VBox([p_value,g] )
Widget1
In [21]:
In [22]:
In [23]:
##compute the best p
# We definitely do need to worry about the range of target(only to 5) and proposal distribution(goes to +oo)
# because by looking at the graph, above 5, pmf is always equal to 0
x=np.arange(6)
PMF_t=np.array([0.4,0.25,0.15,0.1,0.05,0.05])
def c_value(p):
PMF_p=(1-p)**x*p
return np.max(PMF_t/PMF_p)
p_array=np.linspace(0.01,0.99,10000)
c_array=[c_value(p) for p in p_array]
c_array=np.array(c_array)
# use np.argmin to find the index the min c_array, and use that index to find the corresponding p value
p_array[np.argmin(c_array)], np.min(c_array)
p=p_array[np.argmin(c_array)]
c=np.min(c_array)
In [29]:
#generate samples between 0 and 5
#if it is above 5, we discard
def proposal_0_5():
u=np.random.rand()
#geometric distribution
prop=int(np.log(u)/np.log(1-p)) #from direct conversion of inverse transform lecture
while prop>5:
u=np.random.rand()
prop=int(np.log(u)/np.log(1-p))
return prop
def sampling():
# give us a proposal between 0 and 5
prop=proposal_0_5()
#prop could be directly used as index since 0 is the first element
#if it is not from 0, we could do conversion
#(1-p)**prop*p): probability of geometric distribution
#accpetance rate: PMF_t[prop]/((1-p)**prop*p)/c
while np.random.rand()>=PMF_t[prop]/((1-p)**prop*p)/c:
prop=proposal_0_5() # until we get a proposal
return prop
In [30]:
Outcomes_T=np.arange(6)
PMF_T=PMF_t
samples=[sampling() for i in range(10000)]
values,counts=np.unique(samples,return_counts=True)
plt.scatter(values,counts/np.sum(counts),zorder=2)
plt.bar(Outcomes_T, PMF_T )
plt.show()

In [0]: