#!/usr/bin/env python
# coding: utf-8
# # MCMC: Slice sampling
# Florent Leclercq,
# Imperial Centre for Inference and Cosmology, Imperial College London,
# florent.leclercq@polytechnique.org
# In[1]:
import numpy as np
from scipy.integrate import quad
from scipy.stats import norm, cauchy
from math import pi
from matplotlib import pyplot as plt
from matplotlib.ticker import NullFormatter
np.random.seed(123456)
get_ipython().run_line_magic('matplotlib', 'inline')
plt.rcParams.update({'lines.linewidth': 2})
# ## Slice sampling tutorial
# In[2]:
gaussian=norm(0,1)
a=-5.
b=5.
x_arr=np.linspace(a,b,100)
# ### 1- Drawing $y$
# Starting from an initial $x$, we draw $y$ uniformly in $[0,f(x)]$.
# In[3]:
x=1.3
f_of_x=gaussian.pdf(x)
y=np.random.uniform(0.,f_of_x)
# In[4]:
plt.figure(figsize=(10,6))
plt.xlim([a,b])
plt.plot(x_arr,gaussian.pdf(x_arr))
plt.plot([x,x],[0.,f_of_x],color='black')
plt.plot([a,b],[f_of_x,f_of_x],color='black')
plt.plot([a,b],[y,y],color='black',linestyle='--')
plt.fill_between([a,b],0.,f_of_x,facecolor='grey',alpha=0.3, linewidth=0.)
plt.ylim(bottom=0.)
plt.show()
# ### 2- Drawing a new $x$
# We draw $x$ uniformly in the "slice" where $f(x)\geq y$. In the case of the Gaussian distribution $G(0,1)$, this is drawing $x$ uniformly in $[-x_0,x_0]$ where $x_0=\sqrt{-2 \ln(y \sqrt{2\pi})}$. This proposed sample is accepted or rejected according to the usual Metropolis-Hastings rule.
# In[5]:
x0=np.sqrt(-2*np.log(y*np.sqrt(2*pi)))
x_new=np.random.uniform(-x0,x0)
# In[6]:
x_slice_arr=np.linspace(-x0,x0,100)
plt.figure(figsize=(10,6))
plt.xlim([-5.,5.])
plt.plot(x_arr,gaussian.pdf(x_arr))
plt.plot([-5.,5.],[y,y],color='black')
plt.plot([x0,x0],[0.,y],color='black')
plt.plot([-x0,-x0],[0.,y],color='black')
plt.plot([x_new,x_new],[0.,y],color='black',linestyle='--')
plt.fill_between(x_slice_arr,0.,y,facecolor='grey',alpha=0.3, linewidth=0.)
plt.ylim(bottom=0.)
plt.show()
# ### 3- Iterate!
# In[7]:
def slice_sampler_gaussian(Nsteps,x_start):
Naccepted=0
samples=np.zeros(Nsteps+1)
samples[0]=x_start
x=x_start
for i in range(Nsteps):
f_of_x=gaussian.pdf(x)
y=np.random.uniform(0.,f_of_x)
x0=np.sqrt(-2*np.log(y*np.sqrt(2*pi)))
x_p=np.random.uniform(-x0,x0)
a = min(1, gaussian.pdf(x_p)/gaussian.pdf(x))
u = np.random.uniform()
if u < a:
Naccepted+=1
x=x_p
samples[i+1]=x
return Naccepted,samples
# In[8]:
x_start=2.
Nsteps=99
Naccepted,samples=slice_sampler_gaussian(Nsteps,x_start)
# In[9]:
fraction_accepted=float(Naccepted)/Nsteps
fraction_accepted
# In[10]:
plt.figure(figsize=(10,6))
plt.xlim([-5.,5.])
plt.plot(x_arr,gaussian.pdf(x_arr))
markerline, stemlines, baseline = plt.stem(samples,gaussian.pdf(samples),linefmt='-k',markerfmt='k.')
baseline.set_visible(False)
plt.title("Slice sampling")
plt.ylim(bottom=0.)
plt.show()
# ## Slice sampling in practice and multimodal distributions
# In[11]:
def target_pdf(x):
return cauchy(scale=0.5,loc=0.8).pdf(x)+0.5*norm(2.8,0.3).pdf(x)
target_pdf=np.vectorize(target_pdf)
# In[12]:
a=-2.
b=5.
x_arr=np.linspace(a,b,200)
f_arr=target_pdf(x_arr)
# In[13]:
plt.figure(figsize=(10,6))
plt.xlim([a,b])
plt.plot(x_arr,f_arr,color='C2')
plt.title("Target pdf")
plt.ylim(bottom=0.)
plt.show()
# In[14]:
def slice_sampler(target_pdf,Nsteps,x_start,x_width):
Naccepted=0
samples=np.zeros(Nsteps+1)
samples[0]=x_start
x=x_start
for i in range(Nsteps):
y=np.random.uniform(0, target_pdf(x))
lb=x
rb=x
# we build the approximate slice by expanding around the current x
while yy:
# x_p was in the slice, we keep it as a proposed sample
# slice sampling satisfies detailed balance
a = min(1, target_pdf(x_p)/target_pdf(x))
u = np.random.uniform()
if u < a:
Naccepted+=1
x=x_p
samples[i+1]=x
else:
# x was not in the slice, we adjust the boundaries of the approximate slice
if np.abs(x-lb)