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)$。


mcmc-transition马氏链转移和接受概率

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

mcmc-algo-1

上述过程中 $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 算法。

mcmc-algo-2

对于分布 $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)$,也有如下等式

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

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

two-stage-gibbs二维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)$

gibbs-algo-2

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

作者 rickjin

《LDA-math-MCMC 和 Gibbs Sampling(2)》有28条评论
  1. 同样的谢谢 rickjin老师。有个问题,Gibbs采样收敛迭代次数是不是只能设置比较大数字,然后认为在这个次数后目标函数都是收敛的。有没有方法可以估计这个收敛次数?

    [回复]

    rickjin 回复:

    判断收敛是一个难题,MCMC 的一些书上有一些经验的方法, 不过我也没有仔细研究过,所以这个问题我暂时回答不了。

    [回复]

  2. 讲的深入浅出,且对于MCMC,Metropolis等的介绍加入了最初的应用背景,喜欢这种风格。有一个疑问:
    从Gibbs采样的过程来看,是不是只要是Markov Chain 上的各个状态是连通的,那么便可以通过这种\alpha = 1 的逐个状态之间的条件概率转移的方法来得到转移矩阵P,使得最终的平稳分布能够满足p(X)?

    [回复]

    rickjin 回复:

    是的, 你的理解没有问题

    [回复]

    richard 回复:

    random walk收敛除了强连通外,还有非周期性的要求..

    [回复]

  3. 谢谢rick的精妙讲解!

    有一个地方没有看明白,图 https://www.52nlp.cn/wp-content/uploads/2013/01/mcmc-transition.jpg 中,为什么拒绝转移之后,不停留在状态i而是转向$x_0$?

    [回复]

    rickjin 回复:

    你误解那张图的意思了, 黑线表示马氏链中可能的状态转移,红线表示t时刻的状态转移。
    那张图的意思是,在 t 时候, 马氏链从节点 i 沿着红色线以 q(i,j) 的概率转移向 j, 但是以 alpha(i,j) 的概率接受这个转移, 以 1-alpha(i,j) 的概率拒绝转移。

    看来那张图画得有问题,我得改改

    [回复]

    hunter 回复:

    非常感谢rick在百忙之中予回复!

    用图表示这个状态转移,确实比较麻烦.我理解,这里状态与时间是两类不同的值.如果以$\alpha(j|i)$概率接收跳转,t+1时刻的状态为j(如图所示)如果被拒绝,t+1 时刻的状态值仍为i(图中似乎仍没有这层意思,在以时间为横轴的情况下),不知这样理解对不对?

    [回复]

  4. Gibbs sampling是Geman and Geman (1984)提出来的,之所以叫做Gibbs是因为他们研究了Gibbs random field,而此时距离物理学家J. Gibbs逝世已经80年了。

    [回复]

    rickjin 回复:

    谢谢指正, 一会修改

    [回复]

  5. MH算法中的α(y,xt)中的p(y)q(xt|y)/p(xt)p(y|xt),是不是应该为p(y)q(xt|y)/p(xt)q(y|xt)

    [回复]

  6. 我觉得MH算法最后一步X_{t+1} = x_t似乎不对啊。我印象中如果reject,应该重新resample,直到成功而已...所以失败了应该是不更新X_{t+1}而不是让X_{t+1} = x_t吧

    [回复]

    luckyjet 回复:

    我糊涂了,我想的那个算法和这个不是同一个东西

    [回复]

    yeky 回复:

    这个系列非常有帮助。感谢rickjin。不过我也想不明白,MH的最后一步,如果reject, X_{t+1}应该选择除y之外的其他状态才对啊?为什么停留在x_{t}?

    [回复]

    yeky 回复:

    有点想明白了。不管有多少概率质量停留在x_t,整体而言p(x_t=x)是不会变的。但是p(x_{t+1}=x'|x_t=x)是根据不同的a(x,x')选择而改变的。也就是autocorrelation是不同的。

    [回复]

  7. 我是不是糊涂了。我怎么看不懂MH中的第一步了。“采样y~q(x|xt)”这个怎么实现呢?整个算法的目的是从给定分布中得到抽样,但是算法步骤中需要用到另外一个分布的抽样?望指点,非常感谢。

    [回复]

    yeky 回复:

    我觉得q(x|xt)应该是任意选取的吧。只要选择一个有解析式的就可以了。不知道理解的对不对。

    [回复]

    yeky 回复:

    但是Gibbs Sampling中的p(y|x_t)与p(x|y_{t+1})如何抽样呢?

    [回复]

  8. 博主您好,请恕我愚钝,本文刚上来就看不懂了,有两个疑问:
    1) 关于细致平稳条件,以LDA-math-MCMC 和 Gibbs Sampling(1)中提到的那个矩阵,平稳分布不是π=[0.286,0.489,0.225] 么,P矩阵是[ 0.6500 0.2800 0.0700; 0.1500 0.6700 0.1800; 0.1200 0.3600 0.5200],取i =1, j = 2, 那么π(1) * P_{1,2} != π(2)*P_{2,1}呐,难道这个π不是平稳分布吗?
    2)接下来就是关于这个定理的证明,看上去等式的第二项和第三项并无差别啊,
    ∑i=1∞π(i)Pij=∑i=1∞π(j)Pji=∑i=1∞π(j)Pji=π(j)⇒πP=π,请问是遗漏了什么吗?

    [回复]

    fangfei 回复:

    这个公式是π平稳分布的充分而非必要条件,这个公式成立,π一定是平稳分布,因为两边加上求和符号就可以了。但是π是平稳分布,不代表这个公式会成立

    [回复]

  9. 你好,写的很不错,不过可能有个小问题。
    对于Alg. 5,利用q矩阵和x_t来对x_t+1进行采样的时候,当随机的u不满足转移的条件(即跳转到第6行),MCMC算法不应该是把t+1时刻的采样值赋值为x_t,而应该是回到第3行,重新进行采样操作,直到随机的u满足转移条件,返回X_t+1 = y作为t+1时间的采样。
    同样对于Alg. 6,类似的问题也存在。我想作者的思路是没有错的,可能是在写法上不大对,当然也有可能是我自己的理解错误。总之,多多交流!

    [回复]

    Patrick 回复:

    我也这么觉得,求指点。

    [回复]

  10. 证明[细致平稳条件]时,那个连续几个等号的等式,后半部分,我想作者的意思应该是将π(j)提到∑前面,不然等号左右两边是一样的。

    [回复]

  11. 我们想要达到平稳分布,转移矩阵很重要,开篇楼主说我们想要采集p(x)的样本,只要让马链的平稳分布符合p(x)就可以了,但是我们确定下来的矩阵为什么能确保平稳分布是符合p(x)的呢?

    [回复]

  12. MCMC算法中:
    第t个时刻马氏链状态为Xt=xt,采样y〜q(x|xt), 然后再根据a(xt,y)接收或拒绝.

    其实Q'(i,j) = q(i,j)a(i,j) = q(i,j)p(j)q(j,i)可以计算出来,
    为什么不直接按照y〜q'(x|xt)采样. 就不用多一步计算接收拒绝了.

    [回复]

  13. 还有Q'(i,j) = q(i,j)a(i,j) = q(i,j)p(j)q(j,i)
    需证明 Σ_j Q'(i,j) =1 才是马氏转移矩阵

    ∵Q已经是马氏转移矩阵了.
    ∴Σ_j q(i,j) = 1

    ∵p(j)<1, q(j,i)<1
    ∴Σ_j q(i,j)p(j)q(j,i) <1

    我这里哪步错了么?

    [回复]

发表回复

您的电子邮箱地址不会被公开。 必填项已用 * 标注