马尔可夫链蒙特卡洛（ MCMC ）采样

【摘要】传统的蒙特卡洛方法采用随机抽样的方式获得样本，其中大量随机抽取的样本要么被拒绝（拒绝采样）、要么被加权（重要性采样），样本效率不高。因此科学家在思考是否存在一种接受率为 $100%$ 的采样方法。马尔可夫链蒙特卡洛方法真是满足此要求的一种高效采样方法，它充分利用马尔可夫链的可逆性和平稳分布收敛特性，通过一段时间的老化后，所得到的样本能够实现 $100%$ 的接受率。

【原文】 MCMC and Gibbs Sampling

1 问题的提出

随机模拟(或者统计模拟)方法有一个很酷的别名是蒙特卡罗模拟(Monte Carlo Simulation)。这个方法始于20世纪40年代，和原子弹制造的曼哈顿计划密切相关，当时乌拉姆、冯.诺依曼、费米、费曼、Nicholas Metropolis 等，在美国洛斯阿拉莫斯国家实验室研究裂变物质的中子连锁反应的时候，开始使用统计模拟的方法，并在最早的计算机上进行编程实现。

simulation

图 1: 随机模拟与计算机

现代的统计模拟方法最早由数学家乌拉姆提出，被 Metropolis 命名为蒙特卡罗方法，蒙特卡罗是著名的赌场，赌博总是和统计密切关联的，所以这个命名风趣而贴切，很快被大家广泛接受。说起蒙特卡罗方法的源头，可以追溯到18世纪，布丰当年用于计算π的著名的投针实验就是蒙特卡罗模拟实验。统计采样的方法其实数学家们很早就知道，但是在计算机出现以前，随机数生成的成本很高，所以该方法也没有实用价值。随着计算机技术在二十世纪后半叶的迅猛发展，随机模拟技术很快进入实用阶段。对那些用确定算法不可行或不可能解决的问题，蒙特卡罗方法常常为人们带来希望。

monte-carlo-simulation

图 2：蒙特卡罗方法

统计模拟中有一个重要的问题就是给定一个概率分布 $p(x)$，如何在计算机中生成它的样本。一般而言均匀分布 $\operatorname{Uniform}(0,1)$ 的样本是相对容易生成的。例如通过线性同余发生器就可以生成伪随机数。我们用某些确定性算法生成 $\left[0,1\right]$ 之间的伪随机数序列后，这些序列的各种统计指标通常和均匀分布 $\operatorname{Uniform}(0,1)$ 的理论计算结果非常接近。因此，伪随机序列就有比较好的统计性质，可以被当成真实的随机数使用。

sampling

图 3：生成一个概率分布的样本

在统计学领域中常见的经典概率分布（无论是连续的还是离散的分布），基本上都可以基于 $\operatorname{Uniform}(0,1)$ 的样本生成。

例如：正态分布可以通过著名的 $\operatorname{Box-Muller}$ 变换得到：

[Box-Muller 变换]

如果随机变量 $U_1, U_2$ 独立，并且两者皆服从均匀分布 $U_1, U_2 \sim \operatorname{Uniform}[0,1]$，如果随机变量 $Z_0,Z_1$ 满足：

$$
\begin{align*}
Z_0 & = \sqrt{-2\ln U_1} cos(2\pi U_2) \
Z_1 & = \sqrt{-2\ln U_1} sin(2\pi U_2) \tag{1}
\end{align*}
$$

则 $Z_0,Z_1$ 相互独立且服从标准正态分布。

其它几个著名的连续分布，包括指数分布、Gamma 分布、t 分布、F 分布、Beta 分布、Dirichlet 分布等，也都可以通过类似的数学变换得到；离散分布通过均匀分布更加容易生成。

更多通过均匀分布的变换生成统计分布的知识，可以参考统计计算的书籍，其中 Sheldon M. Ross 的《统计模拟》是写得非常通俗易懂的一本。

但是，事实上我们并没有那么幸运，在实际需要解决的大多数问题中，$p(x)$ 的形式很复杂，有些 $p(x)$ 还是高维分布（这在大数据时代更为表现更为突出），利用均匀分布的方法生成样本几乎不太可能。譬如如下常见情况：

在贝叶斯统计中，生成后验分布 $\displaystyle p(x) = \frac{\tilde{p}(x)}{\int \tilde{p}(x) dx}$ 的样本, 此分布中的分子 $\tilde{p}(x)$ 可计算，但分母中的积分无法显式计算，此时无法用均匀分布来抽取样本。
在生成式模型中，生成联合概率分布 $p(x,y)$ 的样本，该分布本身计算很困难，但其条件分布 $p(x|y),p(y|x)$ 的计算可能相对容易; 如果 $p(\mathbf{x})$ 是高维的，这种情形就更加明显。

上述这种困难情况出现时，就需要使用一些新的随机模拟方法来生成样本。其中一种简易的方法，就是构造一个具有解析累积密度函数（ CDF ）的简单形式分布，在该 CDF 上做均匀的随机采样，然后对所抽取的样本做拒绝采样或重要性加权。这就是前面文章中介绍的：接受-拒绝采样、重要性采样。

拒绝采样和重要性采样都采取了样本抽取、接受拒绝（或重要性加权）两个步骤，其中在样本抽取环节大都采取了在均匀分布上随机游走的方法。这两种方法存在的问题是：

拒绝采样方法的接受概率几乎永远不等于 $1$，也就是抽取了大量无效样本，但其优点是最后得到的样本效率高。换句话说，就是：所抽取初始样本的接受概率虽然小于 $100%$，但最终得到的样本频率正比于概率密度（样本效率高）；
重要性采样方法所有样本照单全收，虽然接受概率等于 $1$，但所有样本都要参与计算，其缺点是最后得到的样本效率不高。换句话说，就是：所抽取初始样本的接受概率虽然等于 $100%$，但最终得到的样本频率与概率密度无关（样本效率低）。

因此，科学家们进一步研究了一些新型采样方法，希望既能够像重要性采样那样照单全收（接受概率 $100%$，又能像拒绝采样那样收获高效样本（即样本频率正比于概率密度）。这就是本节中将要介绍的 MCMC (Markov Chain Monte Carlo) 系列算法，其中比较典型的是基础 Metropolis Hastings 算法和面向高维分布的 Gibbs Sampling 算法，在现代贝叶斯分析中被广泛使用。

要了解 MCMC 系列算法，首先要对马尔可夫链的特性有基本的认识。

2 马氏链及其平稳分布

2.1 马氏链的概念

马氏链的数学定义很简单

$$
P(X_{t+1}=x|X_t, X_{t-1}, \cdots) =P(X_{t+1}=x|X_t) \tag{2}
$$

从物理上理解，就是从马氏链上的前一个状态转移到当前状态的概率只依赖于前一个状态（所谓状态指随机量 $X$ 的一个确定性取值 $x$，$X$ 的取值范围可以是无穷）。这意味着：

前一个时刻的确定状态 $i$ 转移到当前时刻的任意新状态，均存在一个转移概率。
前一个时刻的确定状态到当前时刻的所有可能的新状态，可以构成一个转移概率分布。
前一时刻所有可能的确定状态到当前时刻的所有可能的新状态，可以构成一个转移概率矩阵。
对于离散状态，转移概率矩阵有可能是有限维的；对于连续状态，转移概率可能是无限维的，只能用函数表示。

2.2 直观看马氏链的特点

我们来看一个马氏链的具体例子，以便直观了解马氏链的特点。

社会学家经常把人按其经济状况分成三类：下层、中层、上层（对应于三种离散状态），我们用 $1,2,3$ 分别代表这三个阶层。社会学家们发现决定一个人的收入阶层的最重要的因素是其父母的收入阶层。如果一个人的收入属于下层，那么其孩子属于下层收入的概率是 $0.65$, 属于中层收入的概率是 $0.28$, 属于上层收入的概率是 $0.07$。事实上，从父代到子代，收入阶层的变化的转移概率如下：

table-1

markov-transition

使用矩阵的表示方式，转移概率矩阵记为

$$
P = \begin{bmatrix} 0.65 & 0.28 & 0.07 \ 0.15 & 0.67 & 0.18 \ 0.12 & 0.36 & 0.52 \ \end{bmatrix} \tag{3}
$$

这里大家可以注意：

矩阵中每一个 $P_{ij}$ 对应于一种状态转移概率
每一行对应了前一时刻某个确定状态对应的转移概率分布
整个矩阵描述前一时刻所有可能的确定状态对应的所有可能的转移概率。

此时，如果前一时刻的所有可能状态构成了概率分布（或概率密度函数） $\pi_i$ ，则通过转移概率矩阵 $\pi_i P$，我们可以得到当前时刻所有可能状态的概率分布（或概率密度函数） $\pi_{i+1}$。例如，当前某 $0$ 代人处在下层、中层、上层的比例服从状态概率分布 $\pi_0=[\pi_0(1), \pi_0(2), \pi_0(3)]$ 时，则其子女状态概率分布将是 $\pi_1=\pi_0 P$, 其孙子代的状态概率分布比例将是 $\pi_2 = \pi_1 P=\pi_0P^2$, ……, 第 $n$ 代子孙的状态概率分布比例将是 $\pi_n = \pi_{n-1}P = \pi_0P^n$ 。

话不多数，用数据说话。假设初始概率分布为 $\pi_0 = [0.21,0.68,0.11]$，则利用 式（3） 的转移概率矩阵，可以计算前 $n$ 代人的分布状况如下：

table-2

我们发现从第 $7$ 代人开始，这个分布就稳定不变了，这个是偶然的吗？我们换一个初始的状态概率分布 $\pi_0 = [0.75,0.15,0.1]$ ，继续计算前 $n$ 代人的分布状况如下：

table-3

我们发现，到第 $9$ 代人的时候, 分布又收敛了。最神奇的是，两次给定不同的初始概率分布，最终都收敛到概率分布 $\pi=[0.286, 0.489, 0.225]$，也说明收敛行为和初始概率分布 $\pi_0$ 无关，而主要由概率转移矩阵 $P$ 决定。如果我们继续计算 $P^n$，会发现：

$$
P^{20} = P^{21} = \cdots = P^{100} = \cdots =
\begin{bmatrix}
0.286 & 0.489 & 0.225 \
0.286 & 0.489 & 0.225 \
0.286 & 0.489 & 0.225\
\end{bmatrix} \tag{4}
$$

也就是说，当 $n$ 足够大的时候，$P^n$ 矩阵的每一行都稳定地收敛到了 $\pi=[0.286, 0.489, 0.225]$ 。这种收敛现象是绝大多数马氏链的共同行为，下面给出马氏链定理。

2.3 马氏链的性质：收敛定理

如果一个非周期马氏链具有转移概率矩阵 $P$, 且它的任何两个状态之间是连通的（可能的状态可以取无穷多），那么当 $n\rightarrow \infty$ 时，转移概率 $P_{ij}^n$ 存在，且该概率仅与目标状态 $j$ 有关，与前一状态 $i$ 无关，记为 $\displaystyle \lim_{n\rightarrow\infty}P_{ij}^n = \pi(j)$，即有：

1 ) 最终各时刻的状态概率分布都会趋同，即有：
$$
\displaystyle \lim_{n\rightarrow\infty}P_{ij}^n = \pi(j) \tag{5}\
$$

$$
\displaystyle \lim_{n \rightarrow \infty} P^n =\begin{bmatrix} \pi(1) & \pi(2) & \cdots & \pi(j) & \cdots \ \pi(1) & \pi(2) & \cdots & \pi(j) & \cdots \ \cdots & \cdots & \cdots & \cdots & \cdots \ \pi(1) & \pi(2) & \cdots & \pi(j) & \cdots \ \cdots & \cdots & \cdots & \cdots & \cdots \ \end{bmatrix} \tag{6}
$$

2 ) 最终某个状态 $X_j$ 的概率只与其他状态的平稳概率及其转移概率有关：
$$
\displaystyle \pi(j) = \sum_{i=0}^{\infty}\pi(i)P_{ij} \tag{7}
$$

3 ) 最终的平稳分布 $\pi$ 仅与转移概率矩阵 $P$ 有关，与初始状态无关，且 $\pi$ 是方程 $\pi P = \pi$ 的唯一非负解。其中，$\pi = [\pi(1), \pi(2), \cdots, \pi(j),\cdots ], \quad \sum_{i=0}^{\infty} \pi_i = 1$

可以很容易看出： $\pi$ 代表最终稳定状态下的状态概率分布，因此也被称为马氏链的平稳分布。

马氏链的收敛定理非常重要，所有的 MCMC(Markov Chain Monte Carlo) 方法都是以此定理作为理论基础的。定理的证明相对复杂，一般随机过程课本中也不给证明，所以我们就直接用这个定理的结论就好了。

我们对这个定理的内容做一些解释说明：

马氏链支持的状态不要求有限，可以是有无穷多个的；
定理中的 “非周期” 概念不做解释，我们遇到的绝大多数马氏链都是非周期的；
两个状态 $X_i,X_j$ 连通并非指状态 $X_i$ 可以直接一步转移到 $X_j$ ( 即 $P_{ij} > 0$ )，而是指状态 $X_i$ 可以通过有限的 $n$ 步转移到达状态 $X_j$ ( 即 $P_{ij}^n > 0$ )。更进一步的，马氏链的任何两个状态之间是连通的，则暗指存在一个 $n$，使得矩阵 $P^n$ 中的每一个元素 $P^n_{ij}$ 的数值都大于零。

我们用 $X_i$ 表示在马氏链上跳转第 $i$ 步后的状态，如果 $\displaystyle \lim_{n\rightarrow\infty}P_{ij}^n = \pi(j)$ 存在，很容易证明以上定理的第二个结论。首先，根据马尔可夫条件：

$$
\begin{align*}
P(X_{n+1}=j) & = \sum_{i=0}^\infty P(X_n=i) P(X_{n+1}=j|X_n=i) \
& = \sum_{i=0}^\infty P(X_n=i) P_{ij}
\end{align*} \tag{8}
$$

上式两边取极限就得到 $\displaystyle \pi(j) = \sum_{i=0}^{\infty}\pi(i)P_{ij}$

2.4 状态转移过程

从初始的状态概率分布 $\pi_0$ 出发，我们在马氏链上做状态转移。如果将 $i$ 时刻的随机变量记为 $X_i$，其状态概率分布记为 $\pi_i(x)$ ，则有

$$
\begin{align*}
X_0 & \sim \pi_0(x)\
X_1 & \sim \pi_1(x)\
\ldots\
X_i & \sim \pi_i(x), \
\ldots\
\pi_i(x) &= \pi_{i-1}(x)P = \pi_0(x)P^n
\end{align*} \tag{9}
$$

由马氏链的收敛定理，状态概率分布 $\pi_i(x)$ 最终将收敛于平稳分布 $\pi(x)$。假设到第 $n$ 步的时候马氏链收敛了，则有

$$
\begin{align*}
X_0 & \sim \pi_0(x) \ X_1 & \sim \pi_1(x) \ & \cdots \ X_n & \sim \pi_n(x)=\pi(x) \ X_{n+1} & \sim \pi(x) \ X_{n+2}& \sim \pi(x) \ & \cdots
\end{align*} \tag{10}
$$

所以 $X_n,X_{n+1},X_{n+2},\cdots \sim \pi(x)$ 都是同分布的随机变量，当然由于状态转移矩阵的存在，他们之间可能并不独立。

此时，如果我们从任意一个具体的初始状态 $x_0$ 开始，沿着马氏链按照概率转移矩阵做跳转，将得到一个转移序列 $x_0, x_1, x_2, \cdots x_n, x_{n+1}\cdots$，此时由于马氏链的收敛特性，我们会发现到达某个平衡点 $n$ 后，$x_n, x_{n+1},\cdots$ 都将是服从于平稳分布 $\pi(x)$ 的样本，也就是说其接受概率为 $100%$，而这些恰恰是我们之前的诉求：既能像重要性采样一样所有样本照单全收（即接受概率 $100%$），又能像拒绝采样那样得到高效率样本（样本频率与状态的概率值呈正比）。

3 Metropolis Hastings 算法

对于给定的概率分布 $p(x)$, 我们当然希望能有便捷的方式生成它的样本。由于马氏链能收敛到平稳分布，于是产生了一个很漂亮的想法：如果能构造一个转移矩阵为 $P$ 的马氏链，使得该马氏链的平稳分布恰好是 $p(x)$，那么无论我们从任何一个初始状态 $x_0$ 出发沿着马氏链转移，都会得到一个转移序列 $x_0, x_1, x_2, \cdots x_n, x_{n+1}\cdots$，如果马氏链在第 $n$ 步实现了收敛，则该步以后的状态 $x_n, x_{n+1}\cdots$ 都是 $\pi(x)$ 的样本。

此绝妙想法在 1953 年被 Metropolis 想到了，为了研究粒子系统的平稳性质， Metropolis 考虑了物理学中常见的波尔兹曼分布的采样问题，首次提出了基于马氏链的蒙特卡罗方法，即 Metropolis 算法，并在最早的计算机上编程实现。Metropolis 算法是首个普适的采样方法，并启发了一系列 MCMC 方法，所以人们把它视为随机模拟技术腾飞的起点。 Metropolis 的这篇论文被收录在《统计学中的重大突破》中， Metropolis 算法也被遴选为二十世纪的十个最重要的算法之一。

我们接下来介绍的 MCMC 算法是 Metropolis 算法的一个改进变种，即常用的 Metropolis-Hastings 算法。

由上一节的例子和定理，我们期望解决如下问题：

如何能够构造一条能够收敛到目标分布的马尔可夫链？
该链如何能够快速进入平稳状态（即收敛）？

同样根据收敛定理，我们发现似乎转移矩阵 $P$ 是唯一能够决定上述两个问题结果的因素, 所以 MCMC 的关键问题是如何构造转移矩阵 $P$, 以使得该链最终的平稳分布恰好是 $p(x)$，并且收敛更快。如何能做到这一点呢？

3.1 如何能够收敛到目标分布？

科学家们发现，对于给定的分布 $\pi(x)$ ，如果转移概率矩阵 $P$ 能够满足如下 细致平衡条件 ，则其平稳分布就必然是目标分布 $\pi(x)$ ：

如果对于所有可能的状态 $i,j$，非周期马氏链的转移矩阵 $P$ 和分布 $\pi(x)$ 满足

$$
\begin{equation}
\pi(i)P_{ij} = \pi(j)P_{ji} \tag{11}
\end{equation}
$$

则 $\pi(x)$ 是马氏链的平稳分布，而上式被称为 细致平衡条件 (detailed balance condition)。

这个定理是显而易见的，因为细致平衡条件的物理含义就是：对于任何两个状态 $i,j$，从 $i$ 转移到 $j$ 而丢失的概率质量，恰好会被从 $j$ 转移回 $i$ 的概率质量补充回来，所以状态 $i$ 上的概率质量 $\pi(i)$ 是稳定的，从而 $\pi(x)$ 是马氏链的平稳分布。

数学上的证明也很简单，由细致平衡条件可得

$$
\begin{align*}
& \sum_{i=1}^\infty \pi(i)P_{ij} = \sum_{i=1}^\infty \pi(j)P_{ji} = \pi(j) \sum_{i=1}^\infty P_{ji} = \pi(j) \ & \Rightarrow \pi P = \pi \tag{12}
\end{align*}
$$

由于 $\pi$ 是方程 $\pi P = \pi P^{(n)}= \pi$ 的解，所以 $\pi$ 是平稳分布。

小结：

（1）当目标分布和状态转移矩阵之间满足细致平衡条件时，该状态转移矩阵对应的马尔可夫链一定会收敛到目标分布 $\pi$。

（2）细致平衡条件仅仅是平稳分布将收敛至 $\pi$ 的充分条件，而非充要条件。

3.2 如何构造状态转移矩阵？

假设我们已经有了一条马氏链，其对应的转移概率矩阵为 $Q$，$q(i,j)$ 表示从状态 $i$ 转移到状态 $j$ 的概率，也可以写为 $q(j|i)$ 或者 $q(i\rightarrow j)$ 。在一般情况下 $p(i) q(i,j) \neq p(j) q(j,i)$ ，也就是通常不满足细致平衡条件，因此无法保证 $p(x)$ 就是该马氏链的平稳分布。我们可否对该马氏链做一个改造，从而使细致平衡条件成立呢？

譬如，我们可以考虑引入一个新的 $\alpha(i,j)$，使得以下细致平衡条件成立：

$$
\begin{equation*}
p(i) q(i,j)\alpha(i,j) = p(j) q(j,i)\alpha(j,i) \quad \tag{13}
\end{equation*}
$$

取什么样的 $\alpha(i,j)$ 能使以上等式成立呢？按照对称性，最简单的想法可以取：

$$
\alpha(i,j)= p(j) q(j,i)， \quad \alpha(j,i) = p(i) q(i,j) \tag{14}
$$

这样 式 13 式就成立了。即有：

$$
p(i) \underbrace{q(i,j)\alpha(i,j)}{Q’(i,j)} = p(j) \underbrace{q(j,i)\alpha(j,i)}{Q’(j,i)} \quad \tag{15}
$$

也就是说，通过引入一个 $\alpha(i,j)$，我们把一个普通马氏链，改造成了一个具有转移矩阵 $Q^\prime$ 的马氏链，并且 $Q^{\prime}$ 满足细致平衡条件，从而使得转移矩阵 $Q^{\prime}$ 的马氏链对应的平稳分布就是 $p(x)$ 。

具体来说，我们只需在状态转移过程中，将原马尔可夫链每一个状态转移步骤 $\pi_i Q_{ij}$ 改造为 $\pi_i [Q_{ij} \pi_j Q_{ji}]$ 即可保证 $\pi$ 是新马尔可夫链的平稳分布了，而新马尔可夫链的状态转移矩阵是 $Q^\prime = [Q_{ij} \pi_j Q_{ji}]$。即改造后的马尔可夫链满足： $\pi Q^\prime = \pi$

上述过程中，改造 $Q$ 时引入的 $\alpha(i,j)$ 可以被称为接受率，起到了类似于 “接受-拒绝采样” 的效果。其物理意义可以理解为：在原马氏链上，从状态 $i$ 以 $q(i,j)$ 的转移概率跳转到状态 $j$ 时，我们并不是全盘接受转移，而是以 $\alpha(i,j)$ 的概率部分（或全部）接受此转移，进而得到一个基于新转移概率 $Q^{\prime} = [q(i,j)\alpha(i,j)]$ 的马氏链，且其对应的平稳分布正是我们期望的 $p(x)$。

mcmc-transition

马氏链的状态转移和接受概率示意图

3.3 MCMC 算法设计

假设我们已经有一个转移矩阵 $Q$，对应元素为 $q(i,j)$，把以上过程整理一下，就可以得到如下用于概率分布 $p(x)$ 的采样算法。

mcmc-algo-1

为了便于理解，在上述过程中的 $p(x)$ 和 $q(x|y)$ 都是离散的，不过即便这两个分布是连续的，上述算法仍然是有效，进而可以得到更一般的、对于连续概率分布 $p(x)$ 的采样算法，其中 $q(x|y)$ 是任意一个连续二元概率分布对应的条件分布。

以上的 MCMC 采样算法已经能很漂亮的工作了，但存在一个小问题： $Q$ 的马氏链在转移过程中的接受率 $\alpha(i,j)$ 可能偏小，容易造成马氏链原地踏步、拒绝大量的转移，进而导致马氏链遍历完整状态空间所需时间太长，收敛到平稳分布 $p(x)$ 的速度太慢。

3.4 加速转移 – Metropolis-Hastings 算法

有没有办法提升一些接受率呢?

假设 $\alpha(i,j)=0.1, \alpha(j,i)=0.2$, 此时满足细致平衡条件，于是

$$
p(i)q(i,j)\times 0.1 = p(j)q(j,i) \times 0.2 \tag{16}
$$

如果此时将上式两边扩大 $5$ 倍，改写为

$$
p(i)q(i,j) \times 0.5 = p(j)q(j,i) \times 1 \tag{17}
$$

可以看到，接受率被提高了，但细致平衡条件并没有打破！这启发我们可以把细致平衡条件 式 15 中的 $\alpha(i,j),\alpha(j,i)$ 同比例放大，直至其中一个放大到 $1$ ( 满足概率条件 )，这样就可以提高采样过程中的转移接受率。所以可以取：

$$
\alpha(i,j) = \min\left{\frac{p(j)q(j,i)}{p(i)q(i,j)},1\right} \tag{18}
$$

于是，经过对 MCMC 算法中接受率的微小改造，就得到了如下 Metropolis-Hastings 算法 ：

mcmc-algo-2

对于分布 $p(x)$，我们构造转移矩阵 $Q^{\prime}$ 使其满足细致平衡条件

$$
p(x) Q’(x\rightarrow y) = p(y) Q’(y\rightarrow x) \tag{19}
$$

此处 $x$ 并不要求是一维的，对于高维空间的 $p(\mathbf{x})$，如果满足细致平衡条件：

$$
p(\mathbf{x}) Q’(\mathbf{x}\rightarrow \mathbf{y}) = p(\mathbf{y}) Q’(\mathbf{y}\rightarrow \mathbf{x}) \tag{20}
$$

上述 Metropolis-Hastings 算法 一样有效。

4 思考

上一篇文章中，我们指出直接从$p(\boldsymbol{x})$中采样是很困难的，就算是拒绝采样，如果$\boldsymbol{x}$的维度过高，那么从近似分布$q(\boldsymbol{x})$采样后其接受率通常也会很低，导致拒绝采样的效率极低而无法使用。那么很自然的一个问题是：

在MCMC方法中，同样要从转移概率$q(\boldsymbol{y}\leftarrow \boldsymbol{x})$中采样，同样有接受概率$\mathcal{A}(\boldsymbol{y}\leftarrow \boldsymbol{x})$，那为什么MCMC就是实用的，而普通的拒绝采样就是不实用的？

其实，对于直接拒绝采样来说，它要从近似分布$q(\boldsymbol{x})$中直接采样的是整个高维序列，如果$q(\boldsymbol{x})$与精确的$p(\boldsymbol{x})$差得比较大，那么接受概率往往是指数衰减的，因此接受概率极低；而对于MCMC方法来说，我们对$q(\boldsymbol{y}\leftarrow \boldsymbol{x})$的形式没有过多的限制，那么我们可以适当地设计$q(\boldsymbol{y}\leftarrow \boldsymbol{x})$，使得概率分布只集中在跟$\boldsymbol{x}$相似的那些$\boldsymbol{y}$中，换句话说，只有当$\boldsymbol{y}$与$\boldsymbol{x}$很相似时，$q(\boldsymbol{y}\leftarrow \boldsymbol{x})$才不为0，否则就为0。这样一来，从$q(\boldsymbol{y}\leftarrow \boldsymbol{x})$采样的结果跟输入$\boldsymbol{x}$只有小的差别，变化不大的情况下，接受率往往更高些，所以使得拒绝采样成为可能。

所以，说白了，MCMC方法就是将“直接生成$\boldsymbol{x}$”，转换为了“从$\boldsymbol{x}_0$出发，反复对它微调、润色，直到生成符合条件的$\boldsymbol{x}$”这么一个渐进式过程，从而使得它行之有效，这跟上一篇文章提到的“一步步来”的思想是一致的。

接下来，我们将介绍两个MH采样的例子：Gibbs采样和模拟退火，它们都体现了MCMC的这种微调、渐进的思想。

5 吉布斯（Gibbs）采样

虽然 Metropolis-Hastings 算法 适用于高维情形，但由于接受率 $\alpha$ 的存在 ( 通常 $\alpha < 1$ )，高维情况下的算法效率不够高。能否找到一个转移矩阵 $Q$ 使得接受率 $\alpha=1$ 呢？

Gibbs 采样背后的想法是：依靠分布中所有其他变量的值进行采样。例如，如果我们有 $3$ 个随机变量，我们可以使用

$x_1^{s+1} \sim p\left(x_1 \mid x_2^s, x_3^s\right)$
$x_2^{s+1} \sim p\left(x_2 \mid x_1^{s+1}, x_3^s\right)$
$x_3^{s+1} \sim p\left(x_3 \mid x_1^{s+1}, x_2^{s+1}\right)$

这很容易推广到 $D$ 维情况。请注意，如果 $x_i$ 是一个值已知的变量，我们将不对其采样，但它可以用作另一个条件分布的输入。

表达式 $p\left(x_i \mid \boldsymbol{x}_{-i}\right)$ 称为第 $i$ 个变量的 完全条件。一般来说，$x_i$ 可能只依赖于其他一些变量。如果将 $p(\boldsymbol{x})$ 表示为图模型，我们可以通过查看 $i$ 的马尔可夫毯来推断依赖关系（参见无向概率图模型中的条件依赖性）：

$$
x_i^{s+1} \sim p\left(x_i \mid \boldsymbol{x}{-i}^s\right)=p\left(x_i \mid \boldsymbol{x}{\mathrm{mb}(i)}^s\right)
$$

我们可以并行采样一些节点，而不影响正确性。特别是，假设可以创建无向图的着色，使得没有任何相邻节点具有相同的颜色，则可以并行采样相同颜色的所有节点，并依次循环遍历所有颜色。

事实证明，吉布斯采样是 Metropolis-Hastings 的特殊情况，我们使用一系列形式的提议：

$$
q_i\left(\boldsymbol{x}^{\prime} \mid \boldsymbol{x}\right)=p\left(x_i^{\prime} \mid \boldsymbol{x}{-i}\right) \mathbb{I}\left(\boldsymbol{x}{-i}^{\prime}=\boldsymbol{x}_{-i}\right)
$$

也就是说，我们进入了一个新状态，在该状态下，$x_i$ 从其完全条件下采样，但 $x_{-i}$ 保持不变。

现在，我们证明每个此类提议的接受率为 $100%$，因此总体算法的接受率也为 $100%$。我们有

6 模拟退火

另一个例子是模拟退火（Simulated Annealing）算法。前一篇文章中，我们提到了最大化一个目标，其实也可以看成是从这个目标中进行随机采样，其结果便对应着模拟退火算法。

首先，设我们要最大化的函数是$f(\boldsymbol{x})$，并且假设存在某个常数$T > 0$，使得对于任意$0 < \tau \leq T$，都有

$$
\begin{equation}\sum_{\boldsymbol{x}} e^{f(\boldsymbol{x})/\tau} < \infty\end{equation}
$$

这里的$\tau$的物理意义就是温度。这个假设看起来是额外的约束，但事实上模拟退火也就只适用于满足这个条件的场景。由于满足该条件，那么我们就可以构造分布

$$
\begin{equation}p_{\tau}(\boldsymbol{x}) = \frac{e^{f(\boldsymbol{x})/\tau}}{\sum\limits_{\boldsymbol{x}} e^{f(\boldsymbol{x})/\tau}}\end{equation}
$$

假设最大值点是唯一的，那么当$\tau\to 0$时，$p_{\tau}(\boldsymbol{x})$就会退化为one hot分布，也就是只有最大值点的概率为1，如果从$p_{\tau}(\boldsymbol{x})$中采样，那么采样结果就是最大值点了。即使$\tau > 0$，那么依然是最大值点的概率最大，所以从$p_{\tau}(\boldsymbol{x})$中采样，也会有趋于最大值点的趋势。所以，我们可以先固定温度$\tau$，构造一个从$p_{\tau}(\boldsymbol{x})$进行MH采样的随机过程，然后再慢慢让$\tau$下降，那么最终的收敛结果就是最大值点附近了，这就是所谓的“模拟退火”了。

为了进行MH采样，我们需要构建转移矩阵，模拟退火的选择比较简单，就是根据当前状态$\boldsymbol{x}$设计固定数目的$\boldsymbol{y}$作为新的$\boldsymbol{x}$的候选值（这步通常称为“变异”，往往只是对$\boldsymbol{x}$进行简单修改来得到$\boldsymbol{y}$），然后$q(\boldsymbol{y}\leftarrow \boldsymbol{x})$就是从候选中均允随机地选一个。既然是均匀的，那么$q(\boldsymbol{y}\leftarrow \boldsymbol{x})$就是一个常数，并且$q(\boldsymbol{y}\leftarrow \boldsymbol{x})=q(\boldsymbol{x}\leftarrow \boldsymbol{y})$，因此就有

$$
\begin{equation}\mathcal{A}(\boldsymbol{y}\leftarrow\boldsymbol{x}) = \min\left(1, \frac{q(\boldsymbol{x}\leftarrow\boldsymbol{y})p_{\tau}(\boldsymbol{y})}{q(\boldsymbol{y}\leftarrow\boldsymbol{x})p_{\tau}(\boldsymbol{x})}\right) = \min\left(1, e^{[f(\boldsymbol{y}) - f(\boldsymbol{x})]/\tau}\right)\end{equation}
$$

所以，模拟退火就是一种搜索策略，如果$f(\boldsymbol{y}) \geq f(\boldsymbol{x})$，那么就更新$\boldsymbol{x}=\boldsymbol{y}$，就算不是，那么还有一定的概率更新$\boldsymbol{x}=\boldsymbol{y}$。在整个搜索过程中，我们慢慢进行退火（$\tau$慢慢趋于0），在适当的退火策略下，模拟退火几乎总是能找到最大值点的。当然，如何根据具体问题制定“适当的”退火策略，那又是一个值得思考的问题了。

模拟退火

初始状态为$\boldsymbol{x}_0$，初始温度为$\tau_0$，温度按照事先制定的策略下降，$t$时刻状态为$\boldsymbol{x}_t$，温度为$\tau_t$。

通过如下流程采样出$\boldsymbol{x}_{t+1}$：
1、采样$\boldsymbol{y}\sim q(\boldsymbol{y}\leftarrow\boldsymbol{x}_t)$；
2、采样$\varepsilon\sim U[0,1]$；
3、计算$\mathcal{A}(\boldsymbol{y}\leftarrow\boldsymbol{x}_t) = \min\left(1, e^{[f(\boldsymbol{y}) - f(\boldsymbol{x}t)]/\tau_t}\right)$；
4、如果$\varepsilon\leq \mathcal{A}(\boldsymbol{y}\leftarrow\boldsymbol{x}t)$，那么$\boldsymbol{x}{t+1} = \boldsymbol{y}$，否则$\boldsymbol{x}{t+1}=\boldsymbol{x}_t$。

如果简单修改一下，将接受策略改为“如果$f(\boldsymbol{y}) \geq f(\boldsymbol{x}_t)$则接受，否则就拒绝”，那么就变成了更朴素的“爬山法（Hill Climbing）”了。显然，爬山法目标更明确，前期可能收敛更快，但是一旦陷入了局部最大值后，通常就跳不出来了，最终的收敛结果没有模拟退火那么好。当然，也可以通过选择不同的初始值来重复实验，提升爬山法的效果。具体选择什么算法，那就看具体问题而定了。

7 案例

4.1 案例：二维分布的采样

我们先看看二维的情形，假设有一个概率分布 $p(x,y)$，考察 $x$ 坐标相同的两个点 $A(x_1,y_1)$ 和 $B(x_1,y_2)$，我们发现：

$$
\begin{align*}
p(x_1,y_1)p(y_2|x_1) = p(x_1)p(y_1|x_1)p(y_2|x_1) \tag{21}\
p(x_1,y_2)p(y_1|x_1) = p(x_1)p(y_2|x_1)p(y_1|x_1) \tag{22}
\end{align*}
$$

所以得到

$$
\begin{equation*}
p(x_1,y_1)p(y_2|x_1) = p(x_1,y_2)p(y_1|x_1) \quad \tag{23}
\end{equation*}
$$

即

$$
p(A)p(y_2|x_1) = p(B)p(y_1|x_1) \tag{24}
$$

基于以上等式，我们发现，在 $x=x_1$ 这条平行于 $y$ 轴的直线上，如果使用条件分布 $p(y|x_1)$ 做为两点之间的转移概率，那么任何两个点之间的转移均满足细致平衡条件。同样的，如果在 $y=y_1$ 这条直线上任意取两个点 $A(x_1,y_1)$ 和 $C(x_2,y_1)$，也有如下等式：

$$
p(A)p(x_2|y_1) = p(C)p(x_1|y_1) \tag{25}
$$

gibbs-transition

平面上马氏链转移矩阵的构造

于是可以构造平面上任意两点之间的转移概率矩阵 $Q$ 如下：

$$
\begin{align*}
Q(A\rightarrow B) & = p(y_B|x_1) & \text{如果} \quad x_A=x_B=x_1 & \
Q(A\rightarrow C) & = p(x_C|y_1) & \text{如果} \quad y_A=y_C=y_1 & \
Q(A\rightarrow D) & = 0 & \text{其它} &
\end{align*} \tag{26}
$$

有了如上转移矩阵 $Q$，我们很容易验证对平面上任意两点 $X,Y$，均满足细致平衡条件

$$
p(X)Q(X\rightarrow Y) = p(Y) Q(Y\rightarrow X) \tag{27}
$$

于是这个二维空间上的马氏链将收敛到平稳分布 $p(x,y)$。这个算法就是大名鼎鼎的 Gibbs 采样算法，是 Stuart Geman 和 Donald Geman 两兄弟于 1984 年提出来的，之所以叫做 Gibbs 采样 是因为他们将其用于 Gibbs 随机场 研究。 Gibbs 采样算法 在现代贝叶斯分析中占据着非常重要的位置。

gibbs-algo-1

two-stage-gibbs

Gibbs 采样算法中的马氏链的转移示例

以上采样过程中（如上图所示），马氏链的转移只是轮换的沿着坐标轴 $x$ 轴和 $y$ 轴做转移，得到样本 $(x_0,y_0),(x_0,y_1),(x_1,y_1), (x_1,y_2),(x_2,y_2), \cdots$，马氏链收敛后得到的样本就是平稳分布 $p(x,y)$ 的样本，而收敛之前的阶段称为 burn-in 阶段（不仅是吉布斯采样，所有 MCMC 采样方法均这么称呼）。

注意：教科书上的 Gibbs 采样算法大都是坐标轴轮换采样的，但其实这并不是强制要求。最一般的情形可以是，在某个时刻 $t$，可以在 $x$ 轴和 $y$ 轴之间随机的选一个坐标轴，然后按条件概率做转移，马氏链同样是可以收敛的，不过轮换坐标轴的做法是一种比较方便的形式。

4.2 案例：多维分布的采样

以上的过程我们很容易推广到高维的情形，对于 式 23，如果 $x_1$ 变为多维情形 $\mathbf{x_1}$，可以看出推导过程不变，细致平衡条件同样成立:

$$
p(\mathbf{x_1},y_1)p(y_2|\mathbf{x_1}) = p(\mathbf{x_1},y_2)p(y_1|\mathbf{x_1}) \tag{28}
$$

此时转移矩阵 $Q$ 由条件分布 $p(y|\mathbf{x_1})$ 定义。上式只是说明了一根坐标轴的情形，和二维情形类似，很容易验证对所有坐标轴都有类似的结论。所以 $n$ 维空间中对于概率分布 $p(x_1,x_2,\cdots, x_n)$ 可以按照如下方式定义转移矩阵：

如果当前状态为 $(x_1,x_2,\cdots,x_n)$，马氏链转移的过程中，只能沿着坐标轴做转移。沿着 $x_i$ 这根坐标轴做转移的时候，转移概率由条件概率 $p(x_i|x_1, \cdots, x_{i-1}, x_{i+1}, \cdots, x_n)$ 定义；
其它无法沿着单根坐标轴进行的转移，转移概率都设置为 $0$。

于是可以把 Gibbs 采样算法从二维的 $p(x,y)$ 推广到 $n$ 维的 $p(x_1,x_2,\cdots, x_n)$ ：

gibbs-algo-2