# Metropolis¶

In [51]:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.mlab as mlab
%matplotlib inline

###### Target distribution¶
In [52]:
def p(x, y):
g1 = mlab.bivariate_normal(x, y, 1.0, 1.0, -1, -1, -0.8)
g2 = mlab.bivariate_normal(x, y, 1.5, 0.8, 1, 2, 0.6)
return 0.6*g1+28.4*g2/(0.6+28.4)

###### Proposal distribution¶
In [53]:
def q(r):
return r + np.random.normal(size=2)

###### Metropolis¶
In [54]:
N = 100000
s = 10  #thinning
r = np.zeros(2)
p0 = p(r[0], r[1])
print p0
samples = []
for i in xrange(N):
rn = q(r)
pn = p(rn[0], rn[1])
if pn >= p0:
p0 = pn
r = rn
else:
u = np.random.rand()
if u < pn/p0:
p0 = pn
r = rn
if i % s == 0:
samples.append(r)

samples = np.array(samples)

0.00632453617818

In [55]:
plt.scatter(samples[:, 0], samples[:, 1], alpha=0.5, s=1)

'''Plot target'''
dx = 0.01
x = np.arange(np.min(samples), np.max(samples), dx)
y = np.arange(np.min(samples), np.max(samples), dx)
X, Y = np.meshgrid(x, y)
Z = p(X, Y)
CS = plt.contour(X, Y, Z, 10)
plt.clabel(CS, inline=1, fontsize=10)
plt.show()

In [56]:
print np.mean(samples[:,0]), np.mean(samples[:,1])

0.210361197903 0.826203546417

In [57]:
np.histogram(samples[:,0])

Out[57]:
(array([  24,  251, 1463, 2693, 2517, 1821,  929,  256,   41,    5]),
array([-4.90450004, -3.72145056, -2.53840108, -1.35535161, -0.17230213,
1.01074735,  2.19379683,  3.37684631,  4.55989579,  5.74294527,
6.92599475]))
In [58]:
plt.hist(samples[:,0], bins=80, histtype='step',normed=True, color='b', label='x')
plt.hist(samples[:,1], bins=80, histtype='step', normed=True, color='r', alpha=0.5, label='y')
plt.title("Marginals")
plt.xlabel("Value")
plt.ylabel("Probability")
plt.legend()
plt.show()

In [59]:
from sklearn import mixture

def fit_samples(samples):
gmix = mixture.GMM(n_components=2, covariance_type='full')
gmix.fit(samples)
print gmix.means_
colors = ['r' if i==0 else 'g' for i in gmix.predict(samples)]
ax = plt.gca()
ax.scatter(samples[:,0], samples[:,1], c=colors, alpha=0.8)
plt.show()

fit_samples(samples)


[[-0.96433308 -1.03744208]
[ 0.93701863  1.97904132]]

In [61]:
?mlab.bivariate_normal

In [ ]: