Metropolis–Hastings.py

# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
from datetime import datetime
import time
import os

# True Model
w = 5.
b = 3.
sd = 8.0
sample_size = 1000. #NOTE 데이터 갯수
iteration = 50000
burning_phase = iteration // 5
DIR = './experiment_{}/'.format(datetime.now().strftime("%m%d-%H%M%S"))
proposal_std = 1
bin_size = 20

os.makedirs(DIR)

# data
x = np.random.rand(sample_size)*20.
y = w*x + b + np.random.normal(0, sd, size=sample_size)

# show data
plt.scatter(x, y)
plt.savefig(DIR+'MH_data.png')
plt.close()

# likelihood
def loglikelihood(params, x, y):
    w, b, sd = params
    epsilon = y-w*x-b
    # np.sum(np.log(norm.pdf(point/sd)))
    return -sample_size*np.log(np.sqrt(2*np.pi)*sd)-np.sum((epsilon)**2)/(2.*(sd**2))

# prior (parameter 간의 독립을 가정//covariance matrix 사용 안함)
def logPrior(params):
    w = params[0]
    b = params[0]
    sd = params[0]

    # 파라미터 전부 gaussian mean = 0, std = 10 로 prior 줌
    return norm.pdf((w-0)/10) * norm.pdf((b-0)/10) * norm.pdf((sd-0)/10)

# run
alphas = [5]
trace = [[1., 1., 1.]]
accepted_count = 0

start_time = time.time()
for n in range(iteration):
    old_params = trace[-1]

    new_params = np.random.multivariate_normal(old_params, [[proposal_std,0,0],[0,proposal_std,0],[0,0,proposal_std]])

    alpha = np.exp(loglikelihood(new_params, x, y) + logPrior(new_params) - loglikelihood(old_params, x, y) - logPrior(old_params))

    r = min(1, alpha)
    u = np.random.uniform(low=0, high=1)

    if u < r:
        accepted_count += 1
        print "{}\taccepted({})\tunaccept({})\tu:{}\tr:{}".format(new_params, accepted_count, n-accepted_count, u, r)
        trace.append(new_params)
    else:
        print "{}\taccepted({})\tunaccept({})\tu:{}\tr:{}".format(old_params, accepted_count, n-accepted_count, u, r)
        trace.append(old_params)

trace = np.array(trace)
w_result = trace[burning_phase:][:, 0]
b_result = trace[burning_phase:][:, 1]
sd_result = trace[burning_phase:][:, 2]

with open(DIR+'meta.txt', 'w') as f:
    s ="Experiment end! make dir and save results, it takes {} seconds".format(time.time()-start_time)
    print s
    f.write("{0}\niteration: {1}\nburning_phase: {2}\n sample_size: {3}\n accepted_count: {4}\n gaussian prior w, b, sd as mean=10, std=1".format(s, iteration, burning_phase, sample_size, accepted_count))

plt.hist(w_result, bins=bin_size)
plt.savefig(DIR+'w_true_{}.png'.format(w))
plt.close()

plt.hist(b_result, bins=bin_size)
plt.savefig(DIR+'b_true_{}.png'.format(b))
plt.close()

plt.hist(sd_result, bins=bin_size)
plt.savefig(DIR+'sd_true_{}.png'.format(sd))
plt.close()