import numpy as np
import seaborn as sns
sns.set_style('ticks')
import matplotlib.pyplot as plt
from pylab import rcParams
rcParams['figure.figsize'] = 12, 3
# Code ported from slides by Iain Murray, MLSS 2009
# http://homepages.inf.ed.ac.uk/imurray2/teaching/09mlss/slides.pdf
def metropolis(x_init, log_ptilde, iters, sigma):
samples = np.zeros(iters)
n_accept = 0
# initialize current state and unnormalized log prob
x = x_init
log_p_x = log_ptilde(x)
for i in xrange(iters):
# make proposal
prop = x + sigma * np.random.randn()
log_p_prop = log_ptilde(prop)
# compute acceptance probability
p_accept = np.minimum(1., np.exp(log_p_prop - log_p_x))
if np.random.rand() < p_accept:
# accept
x = prop
log_p_x = log_p_prop
n_accept += 1
samples[i] = x
return samples, 1. * n_accept / iters
N = 1000
sigma = 0.1
samples, accept_ratio = metropolis(0, lambda x: -0.5*pow(x, 2), N, sigma)
plt.plot(np.arange(1, N+1), samples)
plt.title('sigma = {:0.1f}, acceptance ratio = {:0.1f}%'.format(sigma, 100*accept_ratio))
sns.despine()
plt.show()
sns.distplot(samples)
sns.despine()
plt.show()
N = 1000
sigma = 1
samples, accept_ratio = metropolis(0, lambda x: -0.5*pow(x, 2), N, sigma)
plt.plot(np.arange(1, N+1), samples)
plt.title('sigma = {:0.1f}, acceptance ratio = {:0.1f}%'.format(sigma, 100*accept_ratio))
sns.despine()
plt.show()
sns.distplot(samples)
sns.despine()
plt.show()
N = 1000
sigma = 100
samples, accept_ratio = metropolis(0, lambda x: -0.5*pow(x, 2), N, sigma)
plt.plot(np.arange(1, N+1), samples)
plt.title('sigma = {:0.1f}, acceptance ratio = {:0.1f}%'.format(sigma, 100*accept_ratio))
sns.despine()
plt.show()
sns.distplot(samples)
sns.despine()
plt.show()