马尔可夫-蒙特卡洛(MCMC)方法

发表于 | 阅读次数：
马尔可夫-蒙特卡洛(MCMC)方法的简单理解，翻译自外文

Markov Chain Monte Carlo Without all the Bullshit

文章翻译于Markov Chain Monte Carlo Without all the Bullshit

我有一个特点：我才不喜欢统计学中的术语，符号和写作风格。因为它复杂而没有必要。比如，当我想要了解马尔可夫链-蒙特卡洛(Markov Chain Monte Carlo)方法时候，Encyclopedia of Biostatistics.中摘要是这么说的。

马尔可夫-蒙特卡洛(MCMC)是一种通过模拟估计复杂模型中统计量的期望的技术。连续随机选择形成马尔可夫链，这个固定分布就是目标分布。他对贝叶斯模型的后验分布的评估非常有效。….

我只能模糊地理解作者在这里所说的内容（因为我提前知道MCMC是啥）。它肯定提到了很多比这篇博客要介绍的更高级的东西。但很难找到MCMC的多余解释。作者这些不明不白的“废话”是行话需要。这么做可能是为了解释更加高级的应用程度，但显然没有必要来定义和分析基本思想。

所以，这是我个人的MCMC描述，是受到John Hopcroft and Ravi Kannan的处理启发。

从分布中提取的问题

马尔可夫蒙特卡洛是用来解决从一个复杂分布中采样问题的技术。让我构建以下场景来解释。我有一个魔法箱，他能很好地估计婴儿名字地概率。我给他一个类似于“李二狗(Malcolm)”地字符串，它能告诉我这个名字准确的概率$P_{李二狗}$，然后你就可以给你下一个孩子取这个名字。所以，这里有一个在所有名字上的分布$D$，它非常了解你的偏好，为了继续研究，这个分布式固定的。

那么问题来了：假如我想有效描述(efficiently draw)在分布$D$下的名字。这个问题就是MCMC想要解决的。为什么是这个问题？因为我不知道你选取名字的步骤，所以我不能模仿你的步骤。你可能想到另一种方式：随时生成一个名字$x$，然后从魔法箱得到$P_x$，然后掷一枚概率为$P_x$硬币来决定要不要。这个方式问题是你有爆炸级别数量$n=|x|$的名字！所以概率$P_x$会非常非常的小，你需要非常长的时间才能得到同一个名字。或者只有几个名字概率非零，然后需要我指数级别的重复才能找到它们。效率低的我选择死亡。

所以这是一个非常严重的问题！我们重新规范下描述以表达清楚。

定义(采样问题)：$D$为一个有限集$X$上的分布。你有一个可以得到概率分布函数$p(x)$的黑盒，它的输出为根据$D$得到的$x \in X$的概率。请设计一个有效随机算法$A$，它能输出$X$中元素，使得输出的$x$概率近似$p_(x)$。更一般的，根据$p_(x)$输出$X$元素的采样。

假设算法$A$只能使用公平的随机硬币，尽管这允许我们有效地模拟任何所需概率的硬币。

需要注意，像这种算法我们能做类似估计随机变量$f:X\rightarrow \mathbb{R}$的期望。我们可以通过解决采样的问题得到一个大采样$S \in X$，然后在采样上计算$f$的平均。当采样很容易的时候， 这就是蒙特卡洛。事实上，马尔可夫是用来解决采样问题的，它允许我们在一次采样上估计$\mathbb{E}(f)$。

但核心问题实际上是一个采样问题，而“马尔可夫链蒙特卡罗”更准确地称为“马尔可夫链采样方法”。所以让我们看看为什么马尔可夫链可能会帮助我们。

随机游走，MCMC的’马尔可夫链’部分

马尔可夫链是本质上是在图上随机游走的术语。

给你一个图$G={V,E}$,和每一条边$e=(u,v) \in E$。给你一个数字$P_{u,v} \in [0,1]$。为了使随机游走有效，$p_{u,v}$的和必须为1。

如果满足我们就可以在$G$上根据概率进行随机游走：从顶点$x_0$开始，然后根据出度的边概率随机选择，然后走到顶点$v_1$。然后尽可能重复。

‘尽可能’是因为任意一个图并不是每个顶点都有出向边。我们需要添加一些限定条件，来在MCMC种应用随机游走。但无论如何，随机游走是明确定义的，我们将$(V,E,{P_e}_{e\in E})$整体对象称为马尔可夫链。

这是一个示例，其中图中的顶点对应于情绪状态。

An example Markov chain; image source http://www.mathcs.emory.edu/~cheung/

在统计学，它们非常重视随机游走的’状态’含义。它们呈边概率为’状态到状态转移’。

马尔可夫链的最有用的主要理论是静态分布定理(有时称为’马尔可夫链的基本定理’)。从直觉上说，对于很长的随机游走，你在顶点V结束的概率与你开始时候的地方无关！所有这些概率一起称为随机游走的静态分布，它由马尔可夫链唯一确定。

然而，由于以上原因，静态分布定理并不是对所有马尔可夫链都适用。重点是图G是强连接的。当你忽略图中边的方向时候，每个顶点都能连接到其它顶点的路径，这称为连通有向图。如果考虑方向，则称为强连接。如果边没有零概率，那么强连接等同于以下属性：

对于每一个顶点$v \in V(G)$，一个从顶点$v$无限的随机游走都会以1的概率返回顶点$v$。

事实上它经常会无限返回。这个属性称为统计下的状态$v$的持久性。我不喜欢这个术语，因为它似乎描述的是一个顶点的属性，对我来说，它描述的是包含该顶点的连通组件的属性。在任何情况下，由于马尔可夫蒙特卡洛，我们会通过设计来确保图是强连接的。

我们使用线性代数来描述静态分布，将转移概率矩阵记为$A$，其中$a_{j,i}=p(i,j)$，边$(i,j)\in E$。矩阵的行列和图$G$的顶点相关，每列$i$表示在随机游走中，从状态$i$到其它状态的概率分布。注意$A$是有向加权图$G$的加权邻接矩阵的转置，其中权重是转移概率(我这么做是因为矩阵向量乘法矩阵在左而不是右，见下)。

这个矩阵允许我用线性代数来很好描述一些东西。在实例中，如果给定一个基本顶点$e_i$作为’’在当前顶点$i$的随机游走’’，$Ae_i$表示了一个顶点在一步随机游走后能到达的$j$个顶点对应的概率。此外，如果有在所有顶点的概率分布$q$，那么$Aq$概率向量表示为：

如果一个随机游走以概率$q_i$在状态$i$,然后$Aq$的第$j$项是指在得到向量$j$后一步随机游走后到达顶点$j$的概率。

以这种方式解释，静态分布就是一个概率分布$\pi$，使得$A\pi=\pi$。换句话说，$\pi$就是具有特征值1的$A$的特征向量。

对于博客的读者来说，这是个一个快速的注意事项：对随机游走的分析正是我们在本博客早期研究用于排名的PageRank算法所作的。在PageRank中，我们称$A$为一个”网络矩阵”，通过随机游走，发现了一个特殊的特征值，其特征向量是我们用来对网页进行排名的的“静态分布”(这用了Perron-Frobenius定理，它说随机游走矩阵具有特殊的特征向量)。我们描述了一种通过迭代乘以$A$来真正找到特征向量的算法。以下定理本质上是该算法的变体，但条件限定更弱一些。对于网络矩阵，在给定需要更强的条件下添加额外的’假’边。

定理：图$G$表示强连接图，其中相关的概率为${p_e}_e \in E$边形成了马尔可夫链。对每个概率向量$x_0$，定义$x_{t+1}=Ax_t$，$t \geq 1$，并且$v_t$表示平均$v_t = \frac{1}{t}\sum_{s=1}^t x_s$。 然后：

  1. 有唯一确定概率向量$\pi$使得$A_\pi = \pi$
  2. 对所有的$x_0$，$\lim_{t\rightarrow \infin v_t = \pi}$

证明。因为$v_t$是概率向量，我们希望随着$|t\rightarrow \infin|$而$|Av_t-v_t|\rightarrow 0$。事实上，我们能做以下展开：
$$
\begin{aligned} A v_{t}-v_{t} &=\frac{1}{t}\left(A x_{0}+A x_{1}+\cdots+A x_{t-1}\right)-\frac{1}{t}\left(x_{0}+\cdots+x_{t-1}\right) \ &=\frac{1}{t}\left(x_{t}-x_{0}\right) \end{aligned}
$$
但是$x_t,x_0$是单位向量，所以它们的不同主要在2，$|Av_t-v_t|\leq \frac{2}{t} \rightarrow 0$。现在这很清晰表明它并不依赖$v_0$。对于唯一性，我们使用Perron-Frobenius定义，它表明这种形式的任何矩阵都具有唯一的(归一化)的特征向量。

另外一点，除了通过实际计算这个平均值或使用本征结算期来计算静态分布之外，还可以通过解析方法将其解析为特定矩阵的逆矩阵。定义$B=A-I_n$，$I_n$是一个$n\times n$的单位阵。假设$C$为删除第一行并添加最后一行的$B$。然后可以明白最后一列$C^{-1}$是$\pi$。此处练习留给读者，因为实践中没有人使用这种方法。

最后一句关于为什么需要上面定理中的所有的$x_t$取平均值。有一个额外的技术条件添加到强连接，称为非周期性。这允许我们加强定理来使$x_t$收敛到静态分布。严格说，非周期性使这样的属性，无论你在哪里开始随机游走，在经过一些足够多的步骤后，随机游走都有可能在每个后续步骤处于每个顶点。作为非周期性失败的图的样例：偶数个顶点上的五项循环。在那种情况下，每隔一个步骤只存在某些顶点的正概率，并且对这两个长期序列求平均给出时间的静态分布。

保证马尔可夫链使非周期性的一种方法是确保在任何顶点有正概率。即，图具有自循环。

构建一个图来游走

回想一下，我们的问题是从有限集$X$的分布得到概率函数$p(x)$。MCMC方法是构建一个马尔可夫链，它静态分布正好是$p$，即使你只能通过黑盒访问。也就是说，你(隐式的)选择了一个图$G$,并(隐式地)选择了边地转移概率来构成静态分布$p$。然后你在图$G$上进行足够长地随机游走，输出你所走过的状态对应的$x$。

如果一个图$G$有正确的竞态分布，这是很简单的(事实上，‘大部分’图都会有效)。难得是，给你一个图，你需要证明随机游走到静态分布的收敛速度要比$x$的增加要快。这超出了本文范围，但如何选择’正确’的图并不难理解。

本文将介绍Metropolis-Hastings算法。输入是对p(x)的黑盒访问，输出是一组隐式定义图上随机游走的规则。它的工作原理如下：选择一些方法来放置$X$格子上，这样每个状态对应${0,1,2,\cdots n}^d$中的一些顶点。然后你给所有的相邻的格子添加顶点(两种直接方式)。比如$n=5,d=2$，和$d=3,n\in {2,3}$为例：

你需要注意保证为$X$选择的顶点没有断开，但在很多应用中，$X$自然已经是一个格子。

然后我们必须描述转换概率。令$r$为格子中顶点的最大度数$r=2d$。假设我们在顶点$i$，想知道下一个步骤，可以按以下操作：

  1. 以$1/r$的概率选择一个邻居$j$。(有可能留在原地)
  2. 如果选择的邻居$j$概率大于$i$,$p(j)\geq p(i)$，就决定去$j$
  3. 其它情况，$p(j)<p(i)$，以$p(j)/p(i)$概率去$j$

我们能更简洁保存边$(i,j)$的概率权重$p_{i,j}$：
$$
\begin{aligned} p_{i, j} &=\frac{1}{r} \min (1, p(j) / p(i)) \ p_{i, i} &=1-\sum_{(i, j) \in E(G) ; j \neq i} p_{i, j} \end{aligned}
$$
对每个顶点$i$都很容易检查其实际上是一个概率分布。所以我们展示了$p(x)$是当前随机游走的静态分布。

需要注意的是，如果每个$x \in X$都有概率分布$v(x)$，对所有的$(x,y\in X)$都有$v(x)P_{x,y} = v(y)P_{y,x}$性质,那么$v$就是静态分布。证明它，固定$x$然后对两边在所有$y$求和。结果就是公式$v(x) = \sum_y v(y)P_{y,x}$，这等于$v=Av$。由于静态分布是满足这个公式的唯一向量，所以$v$是静态分布。

使用$p(i)$来选择是很简单的，因为$p(i)p_{i,j}$和$p(i)p_{j,i}$都等于$\frac{1}{r} \min(p(i),p(j))$。所以完成了，可以根据这些概率随机游走来获得样本。

最后

关于MCMC最后我想说的是，在Metropolis-Hastings图随机游走的过程中，你可以估计函数的期望值$\mathbb{E}(f)$(或者任意竞态分布为$p(x)$的图)。通过定义期望为$\sum_xf(x)p(x)$。

现在需要做的就是计算在随机游走中访问到的状态的$f(x)$的平均值。通过一些其它工作，可以证明这个数量会在随机游走中收敛到静态分布。论证懒得写，但重点是它有效。

我没有从估计函数的期望来描述MCMC是因为核心问题是一个抽样问题。此外MCMC的需要应用都只需要一个样本。比如，MCMC可以用于估计任意凸集的体积。可参考 these lecture notes