尤度比によるメトロポリス・ヘイスティングス法の実験

mh.py

# coding: utf-8
import numpy as np
import math

def compute_log_likelihood(x, mean, stddev):
	return math.log(1 / (stddev * math.sqrt(2.0 * math.pi))) - (x - mean) ** 2 / (2.0 * stddev ** 2)

def main():
	prior_mean = 2.5
	prior_stddev = 0.15
	random_walk = 0.1
	burn_in = 500
	samples = []
	num_acceptance = 0
	num_rejection = 0
	x = -10.0	# 初期点
	for i in xrange(5000):
		ll_old = compute_log_likelihood(x, prior_mean, prior_stddev)
		new_x = x + np.random.normal(0, random_walk)
		ll_new = compute_log_likelihood(new_x, prior_mean, prior_stddev)
		acceptance_rate = math.exp(ll_new - ll_old)
		bernoulli = np.random.uniform(0, 1)
		if bernoulli <= acceptance_rate:
			num_acceptance += 1
			x = new_x
		else:
			num_rejection += 1
		if i % 4 == 0 and i > burn_in:
			samples.append(x)
	samples = np.asanyarray(samples)
	mean = np.mean(samples)
	stddev = np.std(samples)
	print mean, stddev, num_acceptance / float(num_acceptance + num_rejection)
	print samples

if __name__ == "__main__":
	main()