转载

MCMC中的Metropolis–Hastings算法与吉布斯采样

Metropolis–Hastings算法是一种具体的MCMC方法，而吉布斯采样（Gibbs Sampling）是Metropolis–Hastings算法的一种特殊形式。二者在机器学习中具有重要作用，Bishop在他的机器学习经典之作PRML中也专门用了一章的篇幅来介绍随机采样方法。本文将结合R语言实例来探讨这两种算法的相关话题。本文是这个系列文章的最后一篇，主要介绍MCMC中的Metropolis–Hastings算法和吉布斯采样（Gibbs Sampling）方面的内容。

作为本系列文章的组成部分，也作为你阅读本文所必须的预备知识，希望各位读者确认已经对如下文章所谈之话题了然于心：

蒙特卡洛（Monte Carlo）法求定积分
蒙特卡洛采样之拒绝采样（Reject Sampling）
矩阵的极限与马尔科夫链（上）
矩阵的极限与马尔科夫链（下）
Chapman-Kolmogorov Equation
蒙特卡洛马尔科夫链（MCMC）

欢迎关注白马负金羁的CSDN博客 http://blog.csdn.net/baimafujinji ，为保证公式、图表得以正确显示，强烈建议你从该地址上查看原版博文。本博客主要关注方向包括：数字图像处理、算法设计与分析、数据结构、机器学习、数据挖掘、统计分析方法、自然语言处理。

Metropolis–Hastings算法

In statistics, the Metropolis–Hastings algorithm is a MCMC method for obtaining a sequence of random samples from a probability distribution for which direct sampling is difficult.

听到Metropolis–Hastings这个名字，有没有点似曾相识，或者在其他地方也听过。如果你做过数学建模，或者学过优化算法方面的内容，那么在研究模拟退火算法时就应该听过Metropolis–Hastings算法这个名字。本质上来说，模拟退火中的Metropolis–Hastings和MCMC中的Metropolis–Hastings确实是一回事。The original algorithm used in simulated annealing and MCMC’s is due to Metropolis. Later generalized by Hastings. Hastings showed that it is not necessary to use a symmetric proposal distribution, and proposed that the proposed new state can be generated from any q ( y | )

$q(y/ |/ x)$ . Of course, the speed with which we reach the equilibrium distribution will depend on the choice of the proposal function.

Metropolis–Hastings算法的执行步骤如下，注意其中的q就是前面文字描述中的P，即后验概率：

Metropolis–Hastings算法的执行步骤是先随便指定一个初始的样本点 0 $x^0$ ， u $u$ 是来自均匀分布的一个随机数， ∗ $x^*$ 是来自q分布的样本点，样本点是否从前一个位置跳转到 ∗ $x^*$ ，由 α ( ∗ ) $/alpha(x^*)$ 和 u $u$ 的大小关系决定。如果接受，则跳转，即新的采样点为 ∗ $x^*$ ，否则就拒绝，即新的采用点仍为上一轮的采样点。这个算法看起来非常简单，你一定会疑问按照如此的步骤来执行，最终采样分布能收敛到目标分布吗？根据之前的讲解，可知要想验证这一点，就需要检查detailed balance是否满足。下面是具体的证明过程，我们不再赘言。
MCMC中的Metropolis–Hastings算法与吉布斯采样

既然已经从理论上证明Metropolis–Hastings算法的确实可行，下面我们举一个简单的例子来看看实际中Metropolis–Hastings算法的效果。稍微有点不一样的地方是，我们这里并未出现 ( ) P ( ′ | )

$/pi(x)P(x'/ |/ x)$ 这样的后验概率形式，作为一个简单的示例，我们只演示从柯西分布中进行采样的做法。

柯西分布也叫作柯西-洛伦兹分布，它是以奥古斯丁·路易·柯西与亨德里克·洛伦兹名字命名的连续概率分布，其概率密度函数为

( ; 0 , ) = 1 [ 1 + ( ) 2 ] = [ ( − 0 ) 2 + 2 ]

其中 0 $x_0$ 是定义分布峰值位置的位置参数， $/gamma$ 是最大值一半处的一半宽度的尺度参数。 0 = 0 $x_0 = 0$ 且 = 1 $/gamma = 1$ 的特例称为标准柯西分布，其概率密度函数为

( ; 0 , 1 ) = 1 ( 1 + 2 )

下面是在R中进行基于Metropolis–Hastings算法对柯西分布进行采样的示例代码。

rm(list=ls()) ## 清除全部对象
set.seed(201912)
reps=40000

# target density: the cauchy distribution
cauchy<-function(x, x0=0, gamma=1){

  out<-1/(pi*gamma*(1+((x-x0)/gamma)^2))
  return(out)

}

chain<-c(0)

for(i in 1:reps){

    proposal<-chain[i]+runif(1, min=-1, max=1)
    accept<-runif(1)<cauchy(proposal)/cauchy(chain[i])
    chain[i+1]<-ifelse(accept==T, proposal, chain[i])
}

# plot over time
plot(chain, type="l")


plot(density(chain[1000:reps]), xlim=c(-5,5), ylim=c(0,0.4), col="red")
den<-cauchy(seq(from=-5, to=5, by=0.1), x=0, gamma=1)
lines(den~seq(from=-5, to=5, by=0.1), lty=2, col="blue")

下图是每次采样点的数值分布情况。

为了更清晰明确的看到采样结果符合预期分布，我们也绘制了分布的密度图，如下图所示，红色实线是采样分布的密度图，而蓝色虚线则是实际柯西分布的密度图，可见二者吻合地相当好。
MCMC中的Metropolis–Hastings算法与吉布斯采样

当然，作为一个简单的开始，上面的例子中我们并没有涉及到 q ( y | ) $q(y/ |/ x)$ ，接下来的这个例子则会演示涉及后验概率的基于Metropolis–Hastings算法的MCMC。我们将用Metropolis–Hastings算法从瑞利分布（Rayleigh Distribution）中采样，瑞利分布的密度为

( ) = e − 2 2 2 , ≥ 0 , > 0

然后取自由度为 $X_t$ 的卡方分布为proposal function。实现代码如下：

rm(list=ls()) ## 清除全部对象

f <- function(x, sigma) {
    if (any(x < 0)) return (0)
    stopifnot(sigma > 0)
    return((x / sigma^2) * exp(-x^2 / (2*sigma^2)))
}

m <- 40000
sigma <- 4
x <- numeric(m)
x[1] <- rchisq(1, df=1)
k <- 0
u <- runif(m)

for (i in 2:m) {
    xt <- x[i-1]
    y <- rchisq(1, df = xt)
    num <- f(y, sigma) * dchisq(xt, df = y)
    den <- f(xt, sigma) * dchisq(y, df = xt)
    if (u[i] <= num/den) x[i] <- y else {
         x[i] <- xt
         k <- k+1     #y is rejected
     }
}

然后我们要验证一下，我的生产的采样数据是否真的符合瑞利分布。注意在R中使用瑞利分布的相关函数，需要加装VGAM包。下面的代码可以让我们直观地感受到采样结果的分布情况。

> curve(drayleigh(x, scale = 4, log = FALSE), from = -1, to = 15, xlim=c(-1,15), ylim=c(0,0.2), 
+ lty = 2, col="blue",xlab = "", ylab="")
> par(new=TRUE)
> plot(density(x[1000:m]) , xlim=c(-1,15), ylim=c(0,0.2),col="red")

从下图中不难看出，我们的采样点分布确实符合预期。当然，这仅仅是一个演示用的小例子。显然，它并不高效。Œ因为我们采用自由度为

$X_t$ 的卡方分布来作为proposal function，大约有40%的采样点都被拒绝掉了。如果换做其他更加合适proposal function，可以大大提升采样的效率。
MCMC中的Metropolis–Hastings算法与吉布斯采样

吉布斯采样（Gibbs Sampling）

In statistics, Gibbs sampling is a MCMC algorithm for obtaining a sequence of observations which are approximated from a specified multivariate probability distribution , when direct sampling is difficult. 吉布斯采样 can be seen as a special case of the Metropolis-Hastings algorithm. 而这个特殊之处就在于，我们在使用吉布斯采样时，通常是对多变量分布进行采样的。比如说，现在我们有一个分布 p ( z 1 , z 2 , z 3 )

$p(z_1,z_2,z_3)$ ，可见它是定义在三个变量上的分布。这种分布很多，比如一维的正态分布也是关于其均值和方差这两个变量的分布。在算法的第 $t$ 步，我们选出了三个值 z ( ) 1 , z ( ) 2 , z ( ) 3 $z_1^{(t)},z_2^{(t)},z_3^{(t)}$ 。这时既然 z ( ) 2 $z_2^{(t)}$ 和 z ( ) 3 $z_3^{(t)}$ 是已知的，那么我们可以由此获得一个新的 z ( + 1 ) 1 $z_1^{(t+1)}$ ，并用这个新的值替代原始的 z ( ) 1 $z_1^{(t)}$ ，而新的 z ( + 1 ) 1 $z_1^{(t+1)}$ 可以由下面这个条件分布来获得

p ( z 1 | z ( ) 2 , z ( ) 3 )

接下来同样的道理再用一个新的 z ( + 1 ) 2 $z_2^{(t+1)}$ 来替代原始的 z ( ) 2 $z_2^{(t)}$ ，而这个新的 z ( + 1 ) 2 $z_2^{(t+1)}$ 可以由下面这个条件分布来获得

p ( z 2 | z ( + 1 ) 1 , z ( ) 3 )

可见刚刚获得的新值 z ( + 1 ) 1 $z_1^{(t+1)}$ 已经被直接用上了。同理，最后用 z ( + 1 ) 3 $z_3^{(t+1)}$ 来更新 z ( ) 3 $z_3^{(t)}$ ，而新的值由下面这个条件分布来获得

p ( z 3 | z ( + 1 ) 1 , z ( + 1 ) 2 )

按照这个顺序如此往复即可。下面给出更为完整的吉布斯采样的算法描述。
MCMC中的Metropolis–Hastings算法与吉布斯采样

当然如果要从理论上证明吉布斯采样确实可以得到预期的分布，我们就应该考察它是否满足detailed balance。但是前面也讲过，吉布斯采样是Metropolis–Hastings算法的一个特例。所以其实我们甚至无需费力去考察其是否满足detailed balance，如果能够证明吉布斯采样就是Metropolis–Hastings算法，而Metropolis–Hastings算法是满足detailed balance，其实也就得到了我们想要的。下面是证明的过程，其中 x − $/mathbf{x}_{-k}$ 表示一个向量 ( 1 , ⋯ , − 1 , + 1 , ⋯ , D ) $(x_1, /cdots, x_{k-1},x_{k+1},/cdots,x_D)$ ，也就是从中剔除了 $x_k$ 元素。
MCMC中的Metropolis–Hastings算法与吉布斯采样

吉布斯采样的例子网上有很多资料可以参考，这里就不再赘述了。最后说明一点就是，以吉布斯采样为基础实现的MCMC在自然语言处理中的LDA里有重要应用。我们花费了很长的篇幅，用了一个系列的文章最后才把MCMC介绍完毕，如果你正在或者后续准备研究LDA的话，这些内容将是非常必要的基础，也是相当不错的开端。

参考文献

[1] 本文中的英文定义来自于维基百科
[2] http://zhfuzh.blog.163.com/blog/static/1455393872012822854853/
[3] 悉尼科技大学徐亦达博士的机器学习公开课授课材料
[4] 莱斯大学Justin Esarey助理教授( http://jee3.web.rice.edu/index.htm )的公开课资料( https://www.youtube.com/watch?v=j4nEAqUUnVw )
[5] http://staff.ustc.edu.cn/~zwp/teach/Stat-Comp/Lec8.R
[6] Christopher Bishop, Pattern Recognition and Machine Learning, Springer, 2007

原文 http://blog.csdn.net/baimafujinji/article/details/53946367

正文到此结束