如何做Gibbs采样

随机模拟

随机模拟（或者统计模拟）方法最早有数学家乌拉姆提出，又称做蒙特卡洛方法。蒙特卡洛是一个著名的赌场，赌博总是和统计有着密切的关系，所以这个命名风趣而贴切，被广为接受。
随机模拟的一个重要问题就是给定一个概率分布$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)})$