从随机过程到马尔科夫链蒙特卡洛方法

1. Introduction

第一次接触到 Markov Chain Monte Carlo (MCMC) 是在 theano 的 deep learning tutorial 里面讲解到的 RBM 用到了 Gibbs sampling，当时因为要赶着做项目，虽然一头雾水，但是也没没有时间仔细看。趁目前比较清闲，把 machine learning 里面的 sampling methods 理一理，发现内容还真不少，有些知识本人也是一知半解，所以这篇博客不可能面面俱到详细讲解所有的 sampling methods，而是着重讲一下这个号称二十世纪 top 10 之一的算法—— Markov chain Monte Carlo。在介绍 MCMC 之前，我们首先了解一下 MCMC 的 Motivation 和在它之前用到的方法。本人也是初学者，错误在所难免，欢迎一起交流。

这篇博客从零开始，应该都可以看懂，主要内容包括：

随机采样

拒绝采样

重要性采样

Metropolis-Hastings Algorithm

Gibbs Sampling

2. Sampling

我们知道，计算机本身是无法产生真正的随机数的，但是可以根据一定的算法产生伪随机数（pseudo-random numbers)。最古老最简单的莫过于 Linear congruential generator：

# -*- coding=utf8 -*-

# Code from Chapter 14 of Machine Learning: An Algorithmic Perspective
# A simple Gibbs sampler

from pylab import *
from numpy import *

def pXgivenY(y,m1,m2,s1,s2):
return random.normal(m1 + (y-m2)/s2,s1)

def pYgivenX(x,m1,m2,s1,s2):
return random.normal(m2 + (x-m1)/s1,s2)

def gibbs(N=5000):
k=20
x0 = zeros(N,dtype=float)
m1 = 10
m2 = 20
s1 = 2
s2 = 3
for i in range(N):
y = random.rand(1)
# 每次采样需要迭代 k 次
for j in range(k):
x = pXgivenY(y,m1,m2,s1,s2)
y = pYgivenX(x,m1,m2,s1,s2)
x0[i] = x

return x0

def f(x):
return exp(-(x-10)**2/10)

# 画图
N=10000
s=gibbs(N)
x1 = arange(0,17,1)
hist(s,bins=x1,fc='k')
x1 = arange(0,17, 0.1)
px1 = zeros(len(x1))
for i in range(len(x1)):
px1[i] = f(x1[i])
plot(x1, px1*N*10/sum(px1), color='k',linewidth=3)

show()

ps:

modified in 2013.11.1: 偶然发现统计之都有一篇类似的博客，gibbs采样解释得更加详细更加恰当，^_^，请点击这里

参考文献：

[1] PRML, chapter 11

[2] An introduction to MCMC for machine learning, Andrieu, Christophe

[3] 随机模拟的基本思想和常用采样方法（sampling）

[4] youtube 上的讲解 MCMC 的视频