LDA-math-MCMC 和 Gibbs Sampling(2)

3 LDA-math-MCMC 和 Gibbs Sampling(2)
3.2 Markov Chain Monte Carlo

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

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

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

定理：[细致平稳条件] 如果非周期马氏链的转移矩阵$P$和分布$\pi(x)$ 满足
\begin{equation}
\pi(i)P_{ij} = \pi(j)P_{ji} \quad\quad \text{for all} \quad i,j
\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}
= \sum_{i=1}^\infty \pi(j)P_{ji} = \pi(j) \\
& \Rightarrow \pi P = \pi
\end{align*}
由于$\pi$ 是方程 $\pi P = \pi$的解，所以$\pi$是平稳分布。

假设我们已经有一个转移矩阵为$Q$马氏链($q(i,j)$表示从状态 $i$转移到状态$j$的概率，也可以写为 $q(j|i)$或者$q(i\rightarrow j)$), 显然，通常情况下 也就是细致平稳条件不成立，所以 $p(x)$ 不太可能是这个马氏链的平稳分布。我们可否对马氏链做一个改造，使得细致平稳条件成立呢？譬如，我们引入一个 $\alpha(i,j)$, 我们希望
\begin{equation}
\label{choose-alpha}
p(i) q(i,j)\alpha(i,j) = p(j) q(j,i)\alpha(j,i)  \quad (*)
\end{equation}
取什么样的 $\alpha(i,j)$ 以上等式能成立呢？最简单的，按照对称性，我们可以取

于是(*)式就成立了。所以有
\begin{equation}
\label{detailed-balance}
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 (**)
\end{equation}
于是我们把原来具有转移矩阵$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)$的算法。

上述过程中 $p(x),q(x|y)$ 说的都是离散的情形，事实上即便这两个分布是连续的，以上算法仍然是有效，于是就得到更一般的连续概率分布 $p(x)$的采样算法，而 $q(x|y)$ 就是任意一个连续二元概率分布对应的条件分布。

以上的 MCMC 采样算法已经能很漂亮的工作了，不过它有一个小的问题：马氏链$Q$在转移的过程中的接受率 $\alpha(i,j)$ 可能偏小，这样采样过程中马氏链容易原地踏步，拒绝大量的跳转，这使得马氏链遍历所有的状态空间要花费太长的时间，收敛到平稳分布$p(x)$的速度太慢。有没有办法提升一些接受率呢?

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

上式两边扩大5倍，我们改写为

看，我们提高了接受率，而细致平稳条件并没有打破！这启发我们可以把细致平稳条件(**) 式中的$\alpha(i,j),\alpha(j,i)$ 同比例放大，使得两数中最大的一个放大到1，这样我们就提高了采样中的跳转接受率。所以我们可以取

于是，经过对上述MCMC 采样算法中接受率的微小改造，我们就得到了如下教科书中最常见的 Metropolis-Hastings 算法。

对于分布 $p(x)$,我们构造转移矩阵 $Q'$ 使其满足细致平稳条件

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

那么以上的 Metropolis-Hastings 算法一样有效。

3.2 Gibbs Sampling

对于高维的情形，由于接受率 $\alpha$的存在(通常 $\alpha < 1$), 以上 Metropolis-Hastings 算法的效率不够高。能否找到一个转移矩阵Q使得接受率 $\alpha=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) \\
p(x_1,y_2)p(y_1|x_1) = p(x_1)p(y_2|x_1)p(y_1|x_1)
\end{align*}
所以得到
\begin{equation}
\label{gibbs-detailed-balance}
p(x_1,y_1)p(y_2|x_1) = p(x_1,y_2)p(y_1|x_1)  \quad (***)
\end{equation}
即

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

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

于是我们可以如下构造平面上任意两点之间的转移概率矩阵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*}

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

于是这个二维空间上的马氏链将收敛到平稳分布 $p(x,y)$。而这个算法就称为 Gibbs Sampling 算法,是 Stuart Geman 和Donald Geman 这两兄弟于1984年提出来的，之所以叫做Gibbs Sampling 是因为他们研究了Gibbs random field, 这个算法在现代贝叶斯分析中占据重要位置。

二维Gibbs Sampling 算法中的马氏链转移

以上采样过程中，如图所示，马氏链的转移只是轮换的沿着坐标轴 $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 period。额外说明一下，我们看到教科书上的 Gibbs Sampling 算法大都是坐标轴轮换采样的，但是这其实是不强制要求的。最一般的情形可以是，在$t$时刻，可以在$x$轴和$y$轴之间随机的选一个坐标轴，然后按条件概率做转移，马氏链也是一样收敛的。轮换两个坐标轴只是一种方便的形式。

以上的过程我们很容易推广到高维的情形，对于(***) 式，如果$x_1$ 变为多维情形$\mathbf{x_1}$，可以看出推导过程不变，所以细致平稳条件同样是成立的
\begin{equation}
\label{gibbs-detailed-balance-n-dimen}
p(\mathbf{x_1},y_1)p(y_2|\mathbf{x_1}) = p(\mathbf{x_1},y_2)p(y_1|\mathbf{x_1})
\end{equation}
此时转移矩阵 Q 由条件分布 $p(y|\mathbf{x_1})$ 定义。上式只是说明了一根坐标轴的情形，和二维情形类似，很容易验证对所有坐标轴都有类似的结论。所以$n$维空间中对于概率分布 $p(x_1,x_2,\cdots, x_n)$ 可以如下定义转移矩阵

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

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

以上算法收敛后，得到的就是概率分布$p(x_1,x_2,\cdots, x_n)$的样本，当然这些样本并不独立，但是我们此处要求的是采样得到的样本符合给定的概率分布，并不要求独立。同样的，在以上算法中，坐标轴轮换采样不是必须的，可以在坐标轴轮换中引入随机性，这时候转移矩阵 $Q$ 中任何两个点的转移概率中就会包含坐标轴选择的概率，而在通常的 Gibbs Sampling 算法中，坐标轴轮换是一个确定性的过程，也就是在给定时刻$t$，在一根固定的坐标轴上转移的概率是1。