原理

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

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

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

而我们常见的概率分布，无论是连续的还是离散的分布，都可以基于Uniform(0,1)的样本生成，例如正态分布可以通过著名的Box-Muller变换得到

Box-Muller变换：如果随机变量\(U_1, U_2\)独立且\(U_1, U_2 ~ Uniform(0,1)\)，则

\[Z_0 = \sqrt{-2 ln U_1} cos(2\pi U_2)\] \[Z_1 = \sqrt{-2 ln U_1} sin(2\pi U_2)\]

则，\(Z_0, Z_1\)独立且服从标准正态分布。

其它几个著名的连续分布，包括指数分布、Gamma分布、t分布、F分布、Beta分布、Dirichlet分布等等，也都可以通过类似的数学变换得到；离散的分布通过均匀分布更加容易生成（更多可参看Sheldon M. Ross 《统计模拟》）。

不过我们并不是总是这么幸运的，当\(p(x)\)的形式很复杂，或者\(p(x)\)是个高维的分布的时候，样本的生成就可能很困难了。譬如有如下的情况

  1. \(P(x) = \frac{\widetilde{p}(x)}{\int \widetilde{p}(x)dx}\)，而\(\widetilde{p}(x)\)我们是可以计算的，但是底下的积分式无法显式计算。

  2. \(p(x,y)\)是一个二维的分布函数，这个函数本身计算很困难，但是条件分布\(p(x \mid y) p(y \mid x)\)的计算相对简单；如果p(x)是高维的，这种情形就更加明显。

此时就需要使用一些更加复杂的随机模拟的方法来生成样本。而本文介绍的MCMC和之后介绍的Gibbs Sampling算法就是最常用的两种。这两个方法在现代贝叶斯分析中被广泛使用。要了解这两个算法，我们首先要对马氏链的平稳分布的性质有基本的认识。

马氏链及其平稳分布

马氏链的数学定义很简单

\[P(X_{t+1}=x \mid X_t, X_{t-1}, ...) = P(X_{t+1}=x \mid X_t)\]

也就是状态转移的概率只依赖于前一个状态。

我们先来看马氏链的一个具体例子。社会学家经常把人按其经济状况分成3类：下层（lower-class）、中层（middle-class）、上层（upper-class），我们用1，2，3分别代表这三个阶层。社会学家们发现决定一个人的收入阶层的最重要的因素就是其父母的收入阶层。如果一个人的收入属于下层类别，那么他的孩子属于下层收入的概率是0.65，属于中层收入的概率是0.28，属于上层收入的概率是0.07。事实上，从父代到子代，收入阶层的变化的转移概率如下：

      子代  
  State 1 2 3
  1 0.65 0.28 0.07
父代 2 0.15 0.67 0.18
  3 0.12 0.36 0.52

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

\[P = \left[ \begin{matrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \end{matrix} \right] \tag{3}\]

假设当前这一代人处在下层、中层、上层的人的比例是概率分布向量\(\pi_0=[\pi_o(1), \pi_0(2), \pi_0(3)]\)，那么他们的子女的分布比例将是\(\pi_1 = \pi_0 P\)，他们的孙子代的分布比例将是\(\pi_2 = \pi_1 P = \pi_0 P^2, ...\)，第n代子孙的收入分布比例将是\(\pi_n = \pi_{n-1} P = \pi_0 P^n\)。

假设初始概率分布为\(\pi_0 = [0.21, 0.68, 0.11]\)，则我们可以计算前n代人的分布状况如下

第n代人 下层 中层 上层
0 0.210 0.680 0.110
1 0.252 0.554 0.194
2 0.270 0.512 0.218
3 0.278 0.497 0.225
4 0.282 0.490 0.226
5 0.285 0.489 0.225
6 0.286 0.489 0.225
7 0.286 0.489 0.225
8 0.286 0.489 0.225
9 0.286 0.489 0.225
10 0.286 0.489 0.225
… … … …

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

第n代人 下层 中层 上层
0 0.75 0.15 0.1
1 0.522 0.347 0.132
2 0.407 0.426 0.167
3 0.349 0.459 0.192
4 0.318 0.475 0.207
5 0.303 0.482 0.215
6 0.295 0.485 0.220
7 0.291 0.487 0.222
8 0.289 0.488 0.225
9 0.286 0.489 0.225
10 0.286 0.489 0.225
… … … …

我们发现，到第9代人时候，分布又收敛了。最为奇特的是，两次给定不同的初始概率分布，最终都收敛到概率分布\(\pi = [0.286, 0.489, 0.225]\)，也就是说收敛行为和初始概率分布\(\pi_0\)无关。这说明这个收敛行为主要是由概率转移矩阵P决定的。我们计算一下\(P^n\)

\[P^20 = P^21 = ... = P^100 = ... = \left[ \begin{matrix} 0.286 & 0.489 & 0.225 \\ 0.286 & 0.489 & 0.225 \\ 0.286 & 0.489 & 0.225 \end{matrix} \right] \tag{3}\]

我们发现，当n足够大的时候，这个\(P^n\)矩阵的每一行都是稳定地收敛到\(\pi = [0.286, 0.489, 0.225]\)这个概率分布。自然地，这个收敛现象并非是我们这个马氏链独有的，而是绝大多数马氏链的共同行为，关于马氏链的收敛我们有如下定理：

如果一个非周期马氏链具有转移概率矩阵P，且它的任何两个状态是连通的，那么\(\lim_{n \to \infty} P_{ij}^n\)存在且与i无关，记\(\lim_{n \to \infty} P_{ij}^n = \pi(j)\)，我们有

  1. \[\lim_{n \to \infty} P^n = \left[ \begin{matrix} \pi(1) & \pi(2) & ... & \pi(j) & ... \\ \pi(1) & \pi(2) & ... & \pi(j) & ... \\ ... & ... & ... &... & ... \\ \pi(1) & \pi(2) & ... & \pi(j) & ... \\ ... & ... & ... &... & ... \end{matrix} \right] \tag{3}\]
  2. \[\pi(j) = sum_{i=0}^\infty \pi(i) P_{ij}\]
  3. \(\pi\)是方程\(\pi P = \pi\)的唯一非负解，其中
\[\pi = [\pi(1), \pi(2), ..., \pi(j), ...], \sum_{i=0}^\infty \pi_i =1\]

\(\pi\)成为马氏链的平稳分布。

这个马氏链的收敛定理非常重要，所有的MCMC（Markov Chain Monte Carlo）方法都是以这个定理作为理论基础的。

Markov Chain Monte Carlo (MCMC)

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

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

我们接下来介绍的MCMC算法是Metropolis算法的一个改进变种，即常用的Metropolis-Hastings(M-H)算法。由上面的例子和定理我们看到了，马氏链的收敛性质主要由转移矩阵P决定，所以基于马氏链做采样的关键问题是如何构造转移矩阵P，使得平稳分布恰好是我们要的分布p(x)。如何能做到这一点呢？我们主要使用如下定理：

定理（细致平稳条件） 如果非周期马氏链的转移矩阵P和分布\(\pi(x)\)满足

\[\pi(i) P_{ij} = \pi(j) P_{ji}, for all i,j\]

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

其实这个定理是显而易见的，因为细致平稳条件的物理含义就是对于任何两个状态i,j，从i转移出去到j而丢失的概率质量，恰好会被从j转移回i的概率质量补充回来，所以状态i上的概率质量\(\pi(i)\)是稳定的，从而\(\pi(x)\)是马氏链的平稳分布。数学上的证明也很简单，由细致平稳条件可得

\[\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\]

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

假设我们已经有一个转移矩阵为Q的马氏链，\(q(i,j)\)表示从状态i转移到状态j的概率，也可以写为\(q(j\mid i)\)或者\(q(i \to j)\)，显然，通常情况下

\[p(i)q(i,j) \neq p(j)q(j,i)\]

也就是细致平稳条件不成立，所以p(x)不太可能是这个马氏链的平稳分布。我们可否对马氏链做一个改造，使得细致平稳条件成立呢？譬如，我们引入一个\(\alpha(i,j)\)，我们希望

\[p(i)q(i,j)\alpha(i,j) \neq p(j)q(j,i)\alpha(j,i)\]

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

\[\alpha(i,j)=p(j)q(j,i), \alpha(j,i)=p(i)q(i,j)\]

于是等式就成立了。所以有

\[Q'(i,j) = q(i,j)\alpha(i,j), Q'(j,i)=q(j,i)\alpha(j,i)\]

于是我们把原来具有转移矩阵Q的一个很普通的马氏链，改造为了具有转移矩阵\(Q'\)的马氏链，而\(Q'\)恰好满足细致平稳条件，由此马氏链\(Q'\)平稳分布就是p(x).

在改造Q的过程中引入的\(\alpha(i,j)\)称为接受率，物理意义可以理解为在原来的马氏链上，从状态i以q(i,j)的概率跳转到状态j的时候，我们以\(\alpha(i,j)\)的概率接受这个转移，于是得到新的马氏链\(Q'\)的转移概率为\(q(i,j)\alpha(i,j)\)。

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

MCMC采样算法：

  1. 初始化马氏链初始状态\(X_0 = x_0\)
  2. 对于t=0,1,2…,循环一下过程进行采样：
    • 第t个时刻马氏链状态为\(X_t = x_t\)，采样\(y ~ q(x\mid x_t)\)。
    • 从均匀分布采样\(u ~ Uniform(0,1)\)
    • 如果\(u < \alpha(x_t, y) = p(y)q(x_t \mid y)\)则接受转移\(x_t \to y\)，即\(X_{t+1} = y\)
    • 否则不接受转移，即\(X_{t+1} = x_t\)

参考文献