白话马尔科夫链蒙特卡罗方法(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 ightarrow 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 e 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_iinmathbb Z),如果(mathbf x)(mathbf y)只有一个个坐标不同,则二者之间有连边。这样形成的就是一个类似(d)维网格(lattice)的结构,但每条(与坐标轴平行的)直线上的点都形成一个


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

对于相邻的两点(mathbf x)(mathbf y),不失一般性,假设它们第一维坐标不同((x_1 e 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 e 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_{xin S} pi_x)。令(ar S)(S)的补集,若(pi(S) le pi(ar S)),那么导率(Phi(S))就是

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

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

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

[Phi = min_{Ssubset V, S e emptyset} Phi(S). ]

定理

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

[Oleft(frac{ln(1/pi_{min})}{Phi^2epsilon^3} ight). ]

(证明……敲不完……下次也不一定)

例子

一条(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, ]

ε-混合时间的上界就是

[Oleft(frac{ln n}{left(frac1n ight)^2 epsilon^3} ight) = Oleft(frac{n^2ln n}{epsilon^3} ight). ]


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

原文地址:https://www.cnblogs.com/RabbitHu/p/MCMC.html