#!/usr/bin/env python
# coding: utf-8

# In[1]:


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


# In[2]:


# 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


# In[3]:


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()


# In[4]:


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()


# In[5]:


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()


# In[ ]: