MCMC和Gibbs Sampling

本文主要摘抄整理自Rickjin的《LDA数学八卦》

统计模拟中一个很重要的问题就是给定一个概率分布$p(x)$,我们如何在计算机中生成它的样本。一般而言，均匀分布$Uniform(0,1)$的样本是相对容易生成的：我们可以通过随机数生成器生成伪随机数序列，作为平均分布的样本。

一些常见的概率分布，可以基于$Uniform(0,1)$的样本生成，比如正态分布可以通过著名的Box-Muller变化得到

Box-Muller变换 如果随机变量$U_1,U_2$独立且$U_1,U_2~Uniform[0,1]$,令

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

其它几个著名的连续分布包括指数分布、Gamma分布、t分布、F分布、Beta分布、Dirichlet分布等等，也都可以通过类似的数学变换得到。

但是，我们并不总是那么幸运，有的时候p(x)的形式十分复杂，或者是个高维分布的时候，样本的生成就十分困难了，此时就需要用到我们下面介绍的MCMC(Markov Chain Monte Carlo)和Gibbs采样算法了。而要了解这两个算法，我们需要对马氏链的平稳分布有所认识。

马氏链及其平稳分布

马尔科夫链是随机变量$X_1,X_2,X_3,...$的一个数列，在这个数列中$X_{n+ 1}$的分布仅跟$X_n$有关,即

我们来看一个马氏链的具体例子。

社会学家经常把人按照经济状况分成下、中、上三层（我们用1、2、3表示），他们发现，决定一个人的收入阶层的最重要的因素就是其父母的收入阶层：如果一个人的收入属于下层类别，那么它的孩子属于下层收入的概率是0.65,属于中层收入的概率是0.28，属于上层收入的概率是0.07.其转移概率矩阵如下：

如上图，其状态转移矩阵为

假设第一代人处在下、中、上的概率是$\pi_0=[\pi_0(1),\pi_0(2),\pi_0(3)]$,那么他们的子女分布比例将是$\pi _2 = \pi _1 P = \pi _0 P^2,......,$第n代子孙的收入分布比例将是$\pi _n=\pi _{n-1}P=\pi _0P^n$。

假设初始概率分布为$\pi_0 = [0.21,0.68,0.11]$,则我们计算前n带人的分布状况如下

我们发现从第7代人开始，这个分布就稳定不变了，这是偶然的吗？如果你还一个初始的概率分布，就会发现，经过几代人的迭代之后，分布又收敛了。更为奇特的是，虽然两次给定的是不同的初始概率分布，但最终都会收敛到相同的概率分布$\pi = [0.286,0.489,0.225]$,也就是说收敛的结果和初始概率分布$\pi_0$无关。这说明收敛是由概率转移矩阵P决定的。

我们发现，当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. $\begin{equation}\lim_{n \to \infty}P^n= \left[ \begin{array}{ccccc}\pi(1)&\pi(2)&...&\pi(j)&...\\\pi(1)&\pi(2)&...&\pi(j)&...\\...&...&...&...&...\\pi(1)&\pi(2)&...&\pi(j)&...\\...&...&...&...&...\end{array}\right]\end{equation}$;
2. $\pi(j)=\sum_{i=0}^\infty \pi(i)P_{ij}$;
3. $\pi$是方程$\pi P=\pi$的唯一非负解

其中，

$\pi$称为马氏链的平稳分布。

马氏链收敛定理非常重要，所有的MCMC(Markov Chain Monte Carlo)方法都是以这个定理作为理论基础的。几个补充说明：

1. 虽然我们在示例中用到的马氏链的状态空间是有限的，但是，收敛定理其实并不要求状态空间有限，这意味着，对于连续型变量，收敛定理依然有效。
2. 两个状态i,j是连通并非指i可以直接一步转移到j($P_{ij}>0$),而是指i可以通过有限的n步转移到达j。也就是说，马氏链的任何两个状态连通是指存在一个n，使得矩阵$P^n$中的任何一个元素的数值都大于0.
3. 我们不打算证明整个定理，但可以小证一下第二个结论，柿子专检软的捏`(∩_∩)′

上式两边取极限即得

由马氏链的收敛定理，概率分布$\pi_i(x)$将收敛到平稳分布$\pi(x)$,假设到第n步的时候马氏链收敛

所以$X_n,X_{n+1},X_{n+2},...\sim \pi(x)$都是同分布(但是并不独立)的随机变量，如果我们从一个具体的初始状态$x_0$开始，沿着马氏链按照概率转移矩阵做跳转，那么我们得到一个转移序列$x_0,x_1,x_2,...,x_n,x_{n+1},...$,其中，$x_n,x_{n+1},...$都将是平稳分布$\pi(x)$的样本。

Markov Chain Monte Carlo

马氏链能收敛到平稳分布，而我们的目的要生成概率分布p(x)的样本，一个很直观的想法就是: 如果我们能够构造一个平稳分布$\pi(x)=p(x)$的马氏链，那么我们从任何一个初始状态$x_0$出发，得到一个序列$x_0,x_1,x_2,...,x_n,x_{n+1}...，$如果马氏链在第n步已经收敛了，于是我们就得到了$\pi(x)$的样本$x_n,x_{n+1},......$

马氏链的收敛性质主要由转移矩阵P决定，所以基于马氏链做采样的关键问题是如何构造转移矩阵P，使得平稳分布$\pi (x)=p(x)$。我们使用马氏链的细致平稳条件来构造转移矩阵P。

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

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

这个定理的物理意义是显而易见的: 对于任何两个状态i,j, 从i转移出去到j丢失的概率，恰好会被从j转移回i的概率补充回来，所以状态i上的概率$\pi(i)$是稳定的。数学上的证明也十分简单:

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

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

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

取什么样的$\alpha(i,j)$能满足上述条件呢？最简单的，按照对称性，我们可以取

于是我们有

于是我们把原来具有转移矩阵Q的一个很普通的马氏链，改造为了具有转移矩阵Q’的马氏链，而Q’满足细致平稳条件，因此马氏链Q’的平稳分布就是p(x)(即马氏链的初始状态分布！)。

在改造Q的过程中引入的$\alpha(i,j)$称为接受率，它的物理意义是从状态i以q(i,j)的概率跳转到状态j的时候，我们以$\alpha(i,j)$的概率接受这个转移，以$1-\alpha(i,j)$的概率留在原状态。

利用马尔科夫链采样概率分布p(x)的MCMC算法如下所示：

• 初始化马尔科夫链初始状态$X_0=x_0$;
• 对$t=0,1,2,...,$循环以下过程就行采样
  1. 第t个时刻马氏链状态为$X_t=x_t$, 采样$y \sim q(x\|x_t)$
  2. 从均匀分布采样$u \sim Uniform[0,1]$
  3. 如果$u < \alpha(x_t,y) = p(y)q(x_t\|y)$,则接受转移$x_t \to y$,即$X_{t+1} = y$;否则不接受转移，$X_{t+1}=x_t$

Gibbs Sampling

以上的MCMC算法已经能够work了，但是对于高维的情形，由于接受率$\alpha$的存在(通常\alpha <1),以上算法的效率不够高，能够找到一个转移矩阵Q使得接受率$\alpha =1$呢？

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

也就是说，细致平稳条件成立！事实上，在$x=x_1$这条平行于y轴的直线上，任何两个点之间的转移满足细致平稳条件。同样的，如果我们在$y=y_1$这条直线上任意取两点$A(x_1,y_1),C(x_2,y_1)$,这两个点之间的转移也满足细致平稳条件。

于是我们可以按照如下规则构造平面上任意两点之间的转移概率矩阵Q

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

于是这个二维空间上的马氏链将收敛到平稳分布$p(x,y)$。这个算法就称为Gibbs算法，是由物理学家Gibbs首先提出来的。

二维Gibbs Sampling算法

1. 随机初始化$X_0=x_0,Y_0=y_0$;
2. 对$t=0,1,2,...$循环采样
   • $y_{t+1} \sim p(y\|x_t)$;
   • $x_{t+1} \sim p(x\|y_{t+1})$。

如上图所示，马氏链的转移只是轮换的沿着坐标轴x轴和y轴做转移，于是得到样本$(x_0,y_0),(x_0,y_1),(x_1,y_2),(x_2,y_2),...$马氏链收敛后，最终得到的样本就是$p(x,y)$的样本。我们上边的算法虽然是坐标轴轮换采样的，但是这其实是不强制要求的;一般的情形是，在t时刻，在x轴和y轴之间随机的选一个坐标轴，然后按条件概率做转移，马氏链也一样收敛。

由以上过程我们可以看出，如果$x_1$变为多为情形$\vec x_1$，细致平稳条件同样成立

此时转移矩阵Q由条件分布$p(y\|x_1)$定义，所以n为空间中对于概率分布$p(x_1,x_2,...,x_n)$可以如下定义转移矩阵

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

于是我们可以把Gibbs Sampling算法从二维推广到n维，当以上算法收敛后，得到的就是概率分布$p(x_1,x_2,...,x_n)$的样本，当然这些样本并不独立，但是我们要求的是采样得到的样本符合给定的概率分布，并不要求独立。

后记

1. MCMC算法和Gibbs Sampling都能得到指定的概率分布p(x)的分布样本，但是样本之间并不独立;
2. 由于MCMC算法存在接受率$\alpha$,使得在高维的时候效率低下，Gibbs Sampling采用坐标轴轮换的采样方法,巧妙的避开了接受率$\alpha$(其实相当于接受率为1)。