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 [ ]: