原文地址：《如何做Gibbs采样(how to do gibbs-sampling)》
随机模拟
随机模拟（或者统计模拟）方法最早有数学家乌拉姆提出，又称做蒙特卡洛方法。蒙特卡洛是一个著名的赌场，赌博总是和统计有着密切的关系，所以这个命名风趣而贴切，被广为接受。
随机模拟的一个重要问题就是给定一个概率分布\(p(x)\)，如何生成它的样本。均匀分布\(Uniform(0,1)\)的样本可以通过线性同余发生器生成的伪随机数来模拟。常见的概率分布，无论是离散的还是连续的分布，都可以基于\(Uniform(0,1)\)的样本经过变换生成。
然而，当概率分布\(p(x)\)的形式很复杂或者\(p(x)\)是个高维分布的时候，样本的生成就会变得困难。此时就需要一些更加复杂的随机模拟方法来生成样本。接下来将会介绍下MCMC(Markov Chain Monte Carlo)和Gibbs Sampling算法。
马尔科夫链及其平稳分布
马尔科夫链的数学定义如下：
\[P(X_{t+1}=x| X_t, X_{t-1},...) = P(X_{t+1}| X_t)\]
定理 如果一个非周期马尔科夫链具有转移概率矩阵\(P\)，且它的任何两个状态是连通的，那么
$\lim_{n\rightarrow \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].$
- \(\pi(j)=\sum_{i=0}^{\infty}{\pi(i)P_{ij}}\).
- \(\pi\)是方程\(\pi P=\pi\)的唯一非负解，其中\(\pi =[\pi(1), \pi(2), ..., \pi(j), ...], \sum_{i=0}^{\infty}{\pi_i} = 1.\)
\(\pi\)称为马尔科夫链的平稳分布。
收敛到平稳分布
从初始概率分布\(\pi_0\)出发，我们在马尔科夫链上做状态转移，记\(X_i\)的概率分布为\(\pi_i\)，则有
\begin{array}{cc}X_0\sim \pi_0(x)& \\X_i\sim \pi_i(x),& \pi_i(x) = \pi_{i-1}(x)P=\pi_0(x)P^i\end{array}
由马尔科夫链收敛的定理，概率分布\(\pi_i(x)\)将收敛到平稳分布\(\pi(x)\)。假设到第\(n\)步的时候马尔科夫链收敛，则有
\begin{array}{l}X_0 \sim \pi_0(x)\\X_1\sim \pi_1(x)\\...\\X_n\sim \pi_n(x) = \pi(x)\\X_{n+1}\sim \pi(x)\\X_{n+2}\sim \pi(x)\\...\end{array}
所以\(X_n,X_{n+1},X_{n+2},...\sim \pi(x)\)都是属于同分布的随机变量，当然他们并不是独立的。如果从\(x_0\)开始沿着马尔科夫链按照概率转移矩阵跳转，那么我们将会得到一个转移序列\(\{x_0,x_1,...,x_n,x_{n+1},...\}\)，由于马尔科夫链具有收敛性，所以\(\{x_{n},x_{n+1},...\}\)都是平稳分布\(\pi(x)\)的样本。
MCMC
定理(细致平稳条件) 如果非周期马尔科夫链的转移矩阵P和分布\(\pi(x)\)满足
\[\pi(i)P_{ij}=\pi(j)P_{ji}, ~for~all~ i,j\]
则\(\pi(x)\)是马尔科夫链的平稳分布，上式被称为细致平稳条件。
假设我们已经有了一个转移矩阵为\(Q\)的马尔科夫链(\(q(i\rightarrow j)\) 表示从状态\(i\)转移到状态\(j\)的概率)，\(p(x)\)表示当前的概率分布\(\pi\)，显然通常情况下：
\[p(i)q(i\rightarrow j) \ne p(j)q(j\rightarrow i)\]
也就是细致平稳条件不成立，所以\(p(x)\)不太可能是马尔科夫链的平稳分布。为了能够对\(p(x)\)进行抽样，我们必须对马尔科夫链进行改造，使得细致平稳条件成立。比如，引入一个接受率\(\alpha(i,j)\)使得
\[p(i)q(i\rightarrow j)\alpha(i,j) = p(j)q(j\rightarrow i)\alpha(j,i)\]
最简单的按照对称性，我们可以取如下\(\alpha(i,j),\alpha(j,i)\)使得上式成立
\[\alpha(i,j)=p(j)q(j\rightarrow i), ~ \alpha(j,i)=p(i)q(i\rightarrow j)\]
我们可以理解这个改造过程为在原来的马尔科夫链上，从状态\(i\)以概率\(q(i\rightarrow j)\)转移到状态\(j\)时，以概率\(\alpha(i,j)\)接受这个转移。所以新的马尔科夫链\(Q^\prime\)的转移概率为\(q(i\rightarrow j)\alpha(i,j)\)。
上面的改造还有一个小小的问题：马尔科夫链的接受率\(\alpha(i,j)\)在转移过程中可能会偏小，如果这样子采样过程中马尔科夫链容易原地踏步，拒绝大量跳转，致使收敛速度太慢。为了提高采样的接受率，我们可以对接受率进行细微的修改：\[\alpha(i,j)=\min \left \\{ \frac{p(j)q(j\rightarrow i) }{p(i)q(i\rightarrow j )}, 1 \right \\}.\]
Metropolis-Hastings采样算法
- 初始化马尔科夫链初始状态\(X_0=x_0\)
- 对t=0,1,2...,循环以下过程进行采样
- 第t时刻马尔科夫链状态\(X_t=x_t\)，采样\(y\sim q(x|x_t)\)
- 从均匀分布中采样\(u\sim Uniform(0,1)\)
- 如果\(u\lt \alpha(x_t,y)\)则接受转移\(x_t \rightarrow y\)，即\(X_{t+1} = y\)
- 否则不接受转移，即\(X_{t+1}=x_t\)
Gibbs Sampling
上述的Metropolis-Hasting算法也可以被应用与高维的场景，但是由于接受率的存在，算法的效率仍然较低。Gibbs Sampling，找到了一个接受率为1的转移矩阵。考察二维的情形，假设概率分布为\(p(x,y)\)。给定\(x\)坐标相同的两个点\(A(x_1,y_1),B(x_1,y_2)\)不难发现
\[
\begin{array}{c}
p(x_1,y_1)p(y_2|x_1)=p(x_1)p(y_1|x_1)p(y_2|x_1)\\\\
p(x_1,y_2)p(y_1|x_1)=p(x_1)p(y_2|x_1)p(y_1|x_1)\\\\
p(x_1,y_1)p(y_2|x_1) = p(x_1,y_2)p(y_1|x_1)\\\\
\end{array}
\]
也就是\(p(A)p(y_2|x_1)=p(B)p(y_1|x_1)\)。也就是说给定\(x=x_1\)这条平行于\(y\)轴的直线，使用分布\(p(y|x_1\)作为这条直线上任意两点的转移概率，该分布满足细致平稳条件。同样，对于点\(A(x_1,y_1),C(x_2,y_1)\)我们也可以得到\(p(A)p(x_2|y1)=p(C)p(x_1|y_1)\)。
构造如下转移矩阵
\[
\begin{array}{l}
Q(A\rightarrow B) = p(y_B|x_1)\\\\
Q(A\rightarrow C) = p(x_C|y_1)\\\\
Q(A\rightarrow D) = 0
\end{array}
\]
显然上面的转移矩阵\(Q\)满足细致平稳条件。于是这个二维的马尔科夫链将会收敛到平稳分布\(p(x,y)\)。这个算法被称为Gibbs Sampling，是由德国物理学家Gibbs首先提出的。
二维的算法很容易扩展到\(n\)维空间。对于\(n\)维空间的概率分布\(p(x_1,x_2,...,x_n)\)的转移矩阵如下：
- 对于当前状态\((x_1,x_2,...,x_n)\)，马尔科夫链只能沿着坐标轴转移。比如，但状态沿着\(x_i\)转移时，其马尔科夫链转移概率为\(p(x_i|x_1,...,x_{i-1},x_{i+1},...,x_n)\)。
- 其它无法沿着坐标轴进行跳转的转移概率均为0。
n维Gibbs Sampling算法
- 随机选取一个初始状态\((x_1,x_2,...,x_n)\)
- 对t=0,1,2... 循环采样：
- \(x_1^{(t+1)}\sim p(x_1|x_2^{(t)},x_3^{(t)},...,x_n^{(t)})\)
- \(x_2^{(t+1)}\sim p(x_2|x_1^{(t+1)},x_3^{(t)},...,x_n^{(t)})\)
- ...
- \(x_j^{(t+1)}\sim p(x_j|x_1^{(t+1)},x_2^{(t+1)},...,x_{j-1}^{(t+1)},x_{j+1}^{(t)},...,x_n^{(t)})\)
- ...
- \(x_n^{(t+1)}\sim p(x_n|x_1^{(t+1)},x_2^{(t+1)},...,x_{n-1}^{(t+1)})\)