Instantly share code, notes, and snippets.

Embed
What would you like to do?
Metropolis-Hasting Algorithm
# coding:utf-8
# one-dimensional Metropolis-Hasting Algorithm
from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.mlab as mlab
def q(x):
return mlab.normpdf(x, 0, 2)
N = 100000
s = 10
x = 0
p = q(x) # p(0)
samples = np.zeros(N // s)
for i in range(N):
xn = x + np.random.normal() # get j
pn = q(xn) # get p(j)
if pn >= p: # it means p(j) >= p(i),
p = pn # also means alpha >= 1
x = xn
else: # it means p(j) < p(i), also means alpha < 1
u = np.random.rand() # u ~ Uniform(0,1)
if u < pn/p: # u < alpha, accept change
p = pn
x = xn
if i % s == 0:
samples[i // s] = x
plt.hist(samples, bins='auto', normed=True)
dx = 0.01
x = np.arange(np.min(samples), np.max(samples), dx)
y = q(x)
plt.plot(x, y)
plt.show()
# coding:utf-8
# two-dimensional Metropolis-Hasting Algorithm
from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.mlab as mlab
def q(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)
N = 100000
s = 10
r = np.zeros(2)
p = q(r[0], r[1]) # p(0)
samples = np.zeros(shape=(N // s, 2))
for i in range(N):
rn = r + np.random.normal(size=2) # get j
pn = q(rn[0], rn[1]) # get p(j)
if pn >= p: # it means p(j) >= p(i), also means alpha >= 1
p = pn
r = rn
else: # it means p(j) < p(i), also means alpha < 1
u = np.random.rand() # u ~ Uniform(0,1)
if u < pn/p: # u < alpha, accept change
p = pn
r = rn
if i % s == 0:
samples[i // s, :] = r
plt.scatter(samples[:, 0], samples[:, 1], alpha=0.5, s=1)
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 = q(X, Y)
CS = plt.contour(X, Y, Z, 10)
plt.clabel(CS, inline=1, fontsize=10)
plt.show()
@linw1995

This comment has been minimized.

Show comment
Hide comment
Owner

linw1995 commented May 21, 2017

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment