蒙特卡洛马尔科夫链（MCMC）

在以贝叶斯方法为基础的机器学习技术中，通常需要计算后验概率，然后通过最大后验概率（MAP）等方法进行参数推断和决策。然而，在很多时候，后验分布的形式可能非常复杂，这个时候寻找其中的最大后验估计或者对后验概率进行积分等计算往往非常困难，此时可以通过采样的方法来求解。

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

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

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

---

从蒙特卡洛积分开始

在此前的文章《蒙特卡洛（Monte Carlo）法求定积分》一文中我们已经详细介绍了基于随机算法进行定积分求解的技术。这里主要用到其中的平均值法。为了便于后续介绍，这里先做简单回顾。

在计算定积分$\displaystyle\int_a^bg(x)dx$时，我们会把$g(x)$拆解成两项的乘积，即$g(x)=f(x)p(x)$，其中$f(x)$是某种形式的函数，而$p(x)$是关于随机变量$X$的概率分布（也就是PDF或PMF）。如此一来，上述定积分就可以表示为求$f(x)$的期望，即

$$\displaystyle\int_a^bg(x)dx=\int_a^bf(x)p(x)dx=E[f(x)]$$
当然，$g(x)$的分布函数可能具有很复杂的形式，仍然无法直接求解，这时就可以用采样的方法去近似。这时积分可以表示为
$$\displaystyle\int_a^bg(x)dx=E[f(x)]=\frac1n\sum_{i=1}^nf(x_i)$$
在贝叶斯分析中，Monte Carlo 积分运算常常被用来在对后验概率进行积分时做近似估计。比如要计算
$$I(y)=\int f(y\ |\ x)p(x)dx$$
便可以使用下面这个近似形式
$$\hat I(y)=\frac1n\sum_{i=1}^nf(y\ |\ x_i)$$

---

采样方法

不难发现，在利用蒙特卡洛法进行积分求解时，非常重要的一个环节就是从特定的分布中采样。这里的“采样”的意思也即是生成满足某种分布的观测值。就此，本博客之前已经介绍过“逆采样”和“拒绝采样”等方法。如果你对此不甚了解，可以参考：

• 《蒙特卡洛采样之拒绝采样（Reject Sampling） 》

我们采样的目的很多时候都是为了近似积分运算，前面的采样方法（逆采样和拒绝采样）都是先对分布进行采样，然后再用采样的结果近似计算积分。下面要介绍的另外一种方法“重要性采样”（Importance sampling）则两步并做一步，直接的近似计算积分。

我们现在的目标是计算下面的积分

$$E[f(x)]=\int f(x)p(x)dx$$
按照蒙特卡洛求定积分的方法，我们将从满足$p(x)$的概率分布中独立地采样出一系列随机变量$x_i$，然后便有
$$E[f(x)]\approx\frac{1}{N}\sum_{i=1}^Nf(x_i)$$
但是现在的困难是对满足$p(x)$的概率分布进行采样非常困难，毕竟实际中很多$p(x)$的形式都相当复杂。这时我们该怎么做呢？于是想到做等量变换，于是有
$$\int f(x)p(x)dx=\int f(x)\frac{p(x)}{q(x)}q(x)dx$$ 如果把其中的 $\displaystyle f(x)\frac{p(x)}{q(x)}$看成是一个新的函数$h(x)$，则有
$$\int f(x)\frac{p(x)}{q(x)}q(x)dx=\int h(x)q(x)dx\approx\frac{1}{N}\sum_{i=1}^Nf(x_i)\frac{p(x_i)}{q(x_i)}$$ 其中$\displaystyle\frac{p(x_i)}{q(x_i)}$被称为是weights of $x_i$或importance weights。所以这种采样的方法就被称为是重要性采样（Importance sampling）。

如下图所示，当使用Importance Sampling时，我们并不会像在拒绝采样那样拒绝掉某些采样点。此时，所有的采样点都将为我们所用，但是它们的权重是不同的。因为权重为$\displaystyle\frac{p(x_i)}{q(x_i)}$，所以在下图中，当$p(x_i)>q(x_i)$时，采样点$x_i$的权重就大（红色线在绿色在上方时），反之就小。重要性采样就是通过这种方式来从一个“参考分布”$q(x)$中获得“目标分布”$p(x)$的。

---

蒙特卡洛马尔科夫链之基本概念

In statistics, Markov chain Monte Carlo (MCMC) methods are a class of algorithms for sampling from a probability distribution based on constructing a Markov chain that has the desired distribution of its equilibrium distribution.

MCMC构造Markov chain，使其稳态分布等于我们要采样的分布 ，这样我们就可以通过Markov chain来进行采样。这种等价如何来理解是我们深入探讨具体操作方法之前需要先攻克的一个问题。在此之前，希望你对马尔科夫链有一个比较清晰的认识，为此你可以参考

• 《矩阵的极限与马尔科夫链》

我们用下面的式子来表示每一步（时刻推进）中从状态$s_i$到状态$s_j$的转移概率：$P(i,j)=P(i\to j)=Pr(X_{t+1}=s_j\ | \ X_t=S_i)$。这里的一步是指都时间从时刻$t$过渡到下一时刻$t+1$。

Markov chain在时刻$t$处于状态$j$的概率可以表示为$\pi_j(t)=Pr(X_t=s_j)$。用向量$\mathbf{\pi}({t})$来表示在时刻$t$各个状态的概率向量。在初始时刻，需要选取一个特定的$\mathbf{\pi}(0)$，通常情况下我们可以使向量中一个元素为1，其他元素均为0，即从某个特定的状态开始。随着时间的推移，状态会通过转移矩阵，向各个状态扩散。

Markov chain在$t+1$时刻处于状态$s_i$的概率可以通过$t$时刻的状态概率和转移概率来求得，并可通过Chapman-Kolmogorov Equation来得到（可以参考《Chapman-Kolmogorov Equation》）

$$\pi_i(t+1)=Pr(X_{t+1}=s_i)=\sum_k Pr(X_{t+1}=s_i\ |\ X_t=s_k)\cdot Pr(X_t=s_k)\=\sum_kP(k\rightarrow i)\pi_k(t)=\sum_k P(i\ |\ k)\pi_k(t)$$
用转移矩阵写成矩阵的形式如下
$$\mathbf{\pi}(t+1)=\mathbf{\pi}(t)\mathbf{P}$$
其中转移矩中的元素$P(i,j)$表示$P(i\rightarrow j)$。因此 $\mathbf{\pi}(t)=\mathbf{\pi}(0)\mathbf{P}^t$。定义$P_{ij}^{(n)}$为$Pr(X_{t+n}=s_j\ |\ X_t=s_i)$，易见$P_{ij}^{(n)}$表示矩阵$\mathbf{P}^n$中第$ij$个元素。

一条Markov chain有一个平稳分布$\pi^*$ ，是指给定任意一个初始分布$\pi(0)$ ，经过有限步之后，最终都会收敛到平稳分布$\pi^*$ 。平稳分布具有性质$\pi^*=\pi^* \mathbf{P}$ 。可以结合之前文章中谈到的矩阵极限的概念来理解Markov chain的平稳分布。

一条Markov chain拥有平稳分布的一个充分条件是对于任意两个状态$i$和$j$， 其满足detailed balance ：$P(j\rightarrow i)\pi_j=P(i\rightarrow j)\pi_i$ 。可以看到，此时

$$(\pi\mathbf{P})_j=\sum_i \pi_i P(i\rightarrow j)=\sum_i\pi_j P(j\rightarrow i)=\pi_j\sum_iP(j\rightarrow i)=\pi_j$$ 所以$\pi=\pi \mathbf{P}$，明显满足平稳分布的条件。

如果一条Markov chain满足detailed balance，我们就说它是reversible的。在Markov chain中，对于随机变量的取值是连续的情况，上面的这些定义和性质都是类似的。比如这时转移概率为 $P(x\ |\ y)$，那么$\displaystyle\int P(x\ |\ y)dy=1$。Chapman-Kolmogorov Equation 这时可以表示为

$$\pi_t(y)=\int \pi_{t-1}(y)P(x\ |\ y)dy$$
平稳分布的性质表示为（Obviously, at equilibrium, that stationary distribution satisfies）
$$\pi^*(y)=\int \pi^*(y)P(x\ |\ y)dy$$
The detailed balance condition holds when $\pi(x)P(y\ |\ x)=\pi(y)P(x\ |\ y)$
$$\int \pi(x)P(y\ |\ x)dx=\int \pi(y)P(x\ |\ y)dx=\pi(y)\int P(x\ |\ y)dx=\pi(y)$$

最后我们来总结一下MCMC的基本思想。在拒绝采样和重要性采样中，当前生成的样本点与之前生成的样本点之间是没有关系的，它的采样都是独立进行的。然而，MCMC是基于马尔科夫链进行的采样。这就告诉我们，当前的样本点生成与上一时刻的样本点是有关的。如下图所示，假设我们当前时刻生成的样本点是$x$，下一次时刻采样到$x‘$的（转移）概率就是$P(x\ |\ x‘)$，或者写成$P(x\rightarrow x‘)$。我们希望的是这个过程如此继续下去，最终的概率分布收敛到$\pi(x)$，也就是我们要采样的分布。

易见，转移概率就是$P(x\ |\ x‘)$与采样分布$\pi(x)$之间必然存在某种联系，因为我们希望依照$P(x\ |\ x‘)$来进行转移采样最终能收敛到$\pi(x)$。那这个关系到底是怎么样的呢？首先，我们知道给定所有可能的$x$，然后求从$x$转移到$x‘$的概率，所得之结果就是$x‘$出现的概率，如果用公式表示即
$$\pi(x‘)=\int \pi(x)P(x‘\ |\ x)dx$$ 这其实也就是Chapman-Kolmogorov Equation 所揭示的关系。如果马尔科夫链是平稳的，那么基于该马尔科夫链的采样过程最终将收敛到平稳状态
$$\pi^*(x‘)=\int \pi^*(x‘)P(x\ |\ x‘)dx‘$$ 而平稳状态的充分条件又是满足detailed balance：$\pi(x)P(x‘\ |\ x)=\pi(x)P(x\ |\ x‘)$。可见detailed balance的意思是说从$x$转移到$x‘$的概率应该等于从$x‘$转移到$x$的概率，在这样的情况下，采样过程将最终收敛到目标分布。也就是说，我们设计的MCMC算法，只要验证其满足detailed balance的条件，那么用这样的方法做基于马尔科夫链的采样，最终就能收敛到一个稳定状态，即目标分布。

实际应用中有两个MCMC采样算法非常常用，即Metropolis–Hastings算法与吉布斯采样，可以证明吉布斯采样是Metropolis–Hastings算法的一种特殊情况，而且二者都满足detailed balance的条件。

我们将在本系列的最后一篇文章里介绍Metropolis–Hastings算法与吉布斯采样有关的内容。彼时亦将通过一些具体的例子和R语言程序来体会一下这两种采样方法的威力。

---

参考文献

[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] https://theoreticalecology.wordpress.com/2010/09/17/metropolis-hastings-mcmc-in-r/
[6] Christopher Bishop, Pattern Recognition and Machine Learning, Springer, 2007