白话马尔科夫链蒙特卡罗方法(MCMC)

前言

你清茶园不是人待的地方！
里面的个个都是人才，说话又好听——就是我太菜了啥也听不懂，这次期中还考的贼**烂，太让人郁闷了。
最近课上讲这个马尔科夫链蒙特卡罗方法，我也学得一塌糊涂。这时我猛然想起了自己的博客园密码（雾），来更个博客吧。

[Warning] 本人数学水平差劲，下文用词不严谨、缺少部分证明，请酌情阅读。若出锅，欢迎指正。

啥是马尔科夫链？

马尔科夫链(Markov Chain)，简单来说就是一个用来随机游走的有向图，每条边(u, v)的边权\(p_{uv}\)代表“当前在u，下一步走到v”的概率，显然需要

\[p_{uv}\ge 0, \sum_{v}p_{uv}=1. \]

下文中我们假设这个有向图是强连通的，即任取两个点u和v，都存在从u到v的、边权都大于0的路径（当然从v到u的路径也要存在）。

马尔科夫链（也就是这个随机游走过程）的美妙性质在于它收敛。怎么个收敛法呢？

设这个图有n个点，令\(n\)维行向量\(\mathbf p(t)\)表示随机游走了t步之后的概率分布（在时间t分别位于每个点的概率），\(\mathbf p(t)_i\) 就是第t步到点i的概率。初始状态\(\mathbf p(0)\)是随便钦定的。

再定义一个量\(\mathbf a(t)\)，名叫“长期平均概率分布(long-term average probability distribution)”，

\[\mathbf a(t) = \frac1t \sum_{i=0}^{t-1} \mathbf p(i). \]

顾名思义，就是把前\(t\)个时间点的概率分布取个平均。

定理（马尔科夫链基本定理）

• 存在唯一概率分布向量\(\mathbf \pi\)使得\(\mathbf \pi P = \mathbf \pi\)。（这个P就是由\(p_{xy}\)构成的矩阵。）
• \(\lim_{t\rightarrow \infty} \mathbf a(t)\) 存在，而且就等于\(\mathbf \pi\)。

(证明待补充……放假再敲？）

这个定理就很令人开心了——不管钦定初始状态\(\mathbf p(0)\)的你是欧皇还是非酋，只要游走足够多步，\(\mathbf a(t)\)肯定会收敛到唯一的答案。

诶？为啥要定义这个a(t)呢？p(t)自己它不收敛么？
还真不收敛。考虑一下\(w_{12} = 1, w_{21}=1, p(0)=(1, 0)\)。在这个图上游走就是在两个点之间反复横跳，p(t)是不会收敛的！但是a(t)就能收敛到\((\frac12, \frac12)\)。

啥是蒙特卡罗方法？

蒙特卡罗(Monte Carlo)是谁？不是谁，这是个赌场的名字 (=_=|||)。取这个名字大概只是因为……它是随机的，赌场也是随机的？不得不说这个洋气名字实在太劝退了。

其实蒙特卡罗方法就是……抽样估计。小学学的撒豆子求面积啊，蒲丰投针计算圆周率啊，都可以视作蒙特卡罗算法。

啥是马尔科夫链蒙特卡罗方法(MCMC)？

I have a Markov Chain, I have a Monte Carlo, ah, Markov Chain Monte Carlo!

我们已经知道，马尔科夫链是个好东西，它保证能收敛到某个概率分布。现在我们已知一个概率分布，想要构造出相应的马尔科夫链。这有什么用呢？有些时候，概率分布长得比较复杂，直接根据它生成随机变量非常困难，例如想在一个高维空间中的凸多边形中随机取一个点——这时候我们就可以先构造出这个概率分布对应的马尔科夫链，然后在上面随机游走来取点。

不同的马尔科夫链可能收敛到相同的概率分布，但是收敛的速率有快有慢。在应用中，我们肯定希望马尔科夫链收敛得快一点，少游走几步就能获得和想要的概率分布差不多的结果。下面是两个比较好用的构造方法。

Metropolis-Hasting 算法

有的时候又叫Metropolis算法，反正是个一个名字挺长的算法。以下简称MH算法。

和它不好记的名字相反，MH算法描述起来非常简单、非常符合直觉！

首先你要有一个无向连通图。先不管这个图是怎么构建的，很有可能是出题人送给你的（假如你在做相关的练习题的话）。

设概率分布是\(\mathbf p\)，点\(i\)的概率是\(p_i\)。假如你当前在点\(i\)：

• 先从与点\(i\)相邻的所有点中，等概率随机取一个点\(j\);
• 如果\(p_j \ge p_i\)，则立刻走到\(j\)；
• 如果\(p_j < p_i\)， 随机一下，有\(\frac{p_j}{p_i}\)的概率走到\(j\)，否则待在原地不动。

这样\(p_j\)更大的\(j\)更有可能被走到，很符合直觉吧！

把上面的文字策略用公式表示一下，马尔科夫链中点\(i\)到点\(j\)的边权\(p_{ij}\)就是

\[p_{ij} = \frac1r \min(1, \frac{p_j}{p_i}), \]

而剩下的概率都分给了\(p_{ii}\)，

\[p_{ii} = 1 - \sum_{j\ne i} p_{ij}. \]

那么这个马尔科夫链到底能不能收敛到\(\mathbf p\)这个概率分布呢？

引理

令\(\mathbf \pi\)为一个概率分布，\(P\)为马尔科夫链的边权矩阵，如果对任意两点\(x,y\)，均有

\[\pi_x p_{xy} = \pi_y p_{yx}, \]

则\(\mathbf \pi\)就是马尔科夫链收敛到的概率分布。

证明：若\(\pi_x p_{xy} = \pi_y p_{yx}\)，那么枚举\(y\)，对等式两边分别求和，就有\(\pi_x = \sum_y \pi_y p_{yx}\)，即\(\mathbf \pi = \mathbf \pi P\)。根据上面的马尔科夫链基本定理可知收敛到\(\mathbf \pi\)。

把MH算法构造的马尔科夫链代入这个引理，可知确实收敛到\(\mathbf p\)！

（课上讲了一个应用MH算法破译密码的有趣例子，我来不及敲完了，可以看看这篇文章的开头。）

Gibbs 采样

Gibbs 采样(Gibbs Sampling)是另一种MCMC，它使用的是一种特殊的无向图：每个点对应\(d\)维空间中的一个格点\(\mathbf x = (x_1, x_2, \cdots, x_d), x_i\in\mathbb Z\)，如果\(\mathbf x\)和\(\mathbf y\)只有一个个坐标不同，则二者之间有连边。这样形成的就是一个类似\(d\)维网格（lattice）的结构，但每条（与坐标轴平行的）直线上的点都形成一个团。

（一个用来Gibbs取样的无向图，图片来自Blum, Hopcroft, Kannan, "Foundations of Data Science"）

对于相邻的两点\(\mathbf x\)和\(\mathbf y\)，不失一般性，假设它们第一维坐标不同（\(x_1 \ne y_1\)）而剩下的坐标都相同。那么就令

\[p_{\mathbf{xy}} = \frac1d p(y_1|x_2, x_3, \cdots, x_d). \]

而剩余的概率留给\(p_{\mathbf{xx}}\),

\[p_{\mathbf{xx}} = 1-\sum_{\mathbf y \ne x}p_{\mathbf{xy}}. \]

验证：

\[p(x_1|x_2, x_3, \cdots, x_d) p_{\mathbf{xy}} = p(y_1|x_2, x_3, \cdots, x_d) p_{\mathbf{yx}}, \]

故根据条件概率公式有

\[p(x_1) p_{\mathbf{xy}} = p(y_1) p_{\mathbf{yx}}, \]

所以概率分布确实收敛到\(\mathbf p.\)

ε-混合时间(ε-mixing time)

（我没查到这个ε-mixing time怎么翻译成中文……有大佬知道的话评论区告诉我一下qaq）

尽管我们知道马尔科夫链早晚会收敛的，但它到底早收敛还是晚收敛，对我们很重要。为了衡量收敛的快慢，定义ε-混合时间(ε-mixing time)为\(t\)的最小值，满足：对任意初始状态\(\mathbf p(0)\),\(\|\mathbf a(t) - \mathbf \pi(t)\|_1 < \epsilon\)（这个范数就是曼哈顿距离那种，各维度距离直接相加的和）。为了保证算法的效率，我们想知道这个ε-混合时间的一个上界。

一个上界是由导率(conductance)确定的。啥是导率呢？顾名思义，就和物理上的电导率差不多。对于节点集合\(S\)，令\(\pi(S) = \sum_{x\in S} \pi_x\)。令\(\bar S\)为\(S\)的补集，若\(\pi(S) \le \pi(\bar S)\)，那么导率\(\Phi(S)\)就是

\[\Phi(S) = \frac{\sum_{(x, y)\in(S, \bar S)}\pi_xp_{xy}}{\pi(S)} = \sum_{x\in S}\frac{\pi_x}{\pi(S)}\sum_{y\in\bar S}p_{xy}, \]

即“已知当前在\(S\)中，下一步跳出\(S\)的概率”。

直觉上，如果\(\Phi(S)\)都很大，那么我们可以在图上来去自如，收敛就比较快；反之，如果某个\(\Phi(S)\)很小，一旦进入\(S\)就被困住、很难逃脱，收敛就可能很慢。
因此，定义整个马尔科夫链的导率为

\[\Phi = \min_{S\subset V, S\ne \emptyset} \Phi(S). \]

定理

任何无向图随机游走的ε-混合时间有上界

\[O\left(\frac{\ln(1/\pi_{min})}{\Phi^2\epsilon^3}\right). \]

（证明……敲不完……下次也不一定）

例子

一条\(n\)个节点顺序组成的链，\(\mathbf \pi = (\frac1n, \frac1n, \cdots, \frac1n)\)。

可以构造马尔科夫链\(p(1, 1) = p(n, n) = \frac12, p(i, i+1) = p(i + 1, i) = \frac12, \forall 1 \le i \le n-1\)。

最难逃脱的\(S\)就由前半条链（后半条也行）构成的集合，\(\pi(S) = \frac12\),

\[\Phi = \Phi(S)= \frac{\frac1n \cdot \frac12}{\frac12} = \frac1n, \]

ε-混合时间的上界就是

\[O\left(\frac{\ln n}{\left(\frac1n\right)^2 \epsilon^3}\right) = O\left(\frac{n^2\ln n}{\epsilon^3}\right). \]

---

写作业去啦 未完待续_(:з」∠)_