如何快速理解马尔科夫链蒙特卡洛法？

©PaperWeekly 原创 · 作者｜邓妍蕾

学校｜香港大学硕士

研究方向｜NLP、语音识别

概览

马尔科夫蒙特卡洛法（Markov Chain Monte Carlo, MCMC）经常用在贝叶斯概率模型的推理和学习中，主要是为了解决计算困难的问题。

日常中我们会用采样的方法采集样本，进行近似的数值计算，比如计算样本均值，方差，期望。虽然许多常见的概率密度函数（t 分布，均匀分布，正态分布），我们都可以直接在 Numpy，Scikit-Learn 中找到，但很多时候，一些概率密度函数，特别是高维的概率密度函数，并不是常见的分布，这个时候我们就需要用到 MCMC。 

在开始马尔科夫蒙特卡洛法之前，我们先简单介绍一些蒙特卡洛法和马尔科夫链。

蒙特卡洛法 Monte Carlo

蒙特卡洛法是比较宽泛的一系列算法的统称（想了解可自行 google），它的特点是假设概率分布已知，通过重复的随机采样来获得数值结果。比如根据大数定理，我们可以用采样得到的样本计算得到的样本均值来估计总体期望（例子 1）。又比如，积分的运算往往可以表示为随机变量在一个概率密度函数分布上的期望（例子 2）。 

例子 1：假设有随机变量 x，定义域 X，其概率密度函数为 p(x)，f(x) 为定义在 X 上的函数，目标是求函数 f(x) 关于密度函数 p(x) 的数学期望 。蒙特卡洛法根据概率分布 p(x) 独立地抽样 n 个样本 x1,x2,…..xn，得到近似的 f(x) 期望为：

例子 2： 假设我们想要求解 h(x) 在 X 上的积分：

我们将 h(x) 分解成一个函数 f(x) 和一个概率密度函数 p(x) 的乘积，进而又将问题转换为求解函数 f(x) 关于密度函数 p(x) 的数学期望 ：

这里，f(x) 表示为 ，则有：

更一般的，假设我们想要求解 ，熟悉积分的同学肯定已经知道答案为 ，那么如何用采样的方法来得到这个值呢？

令 ，0<x<10，那么 。

import random
numSamples = 10000
# 在 0-10的均匀分布内采样
samples = [random.uniform(0,10) for _ in range(num_samples)]
f_samples = [10 * sample ** 2 for sample in samples]
result = 1/10000.0 * sum(f_samples)
#>>> result
#333.7766822849899

对于复杂的 h(x)，这种方法计算起来显然就更加方便了（特别是忘记积分怎么算的同学）。

到这里为止，我们简单的介绍了蒙特卡洛方法，但是依旧没有提到要怎么利用复杂的概率密度函数进行采样。接下来我们来看一下接受-拒绝法（accept-reject sampling method），它也是蒙特卡洛法中的一种类型适用于不能直接抽样的情况。 

接受-拒绝法 

假设有一个非常复杂不常见的分布 p(x)，我们没有现成的工具包可以调用后进行采样，那么我们可以用我们已经有的采样分布比如高斯分布作为建议分布(proposal distribution)，用 q(x) 表示，来进行采样，再按照一定的方法来拒绝当前的样本，使得最后得到的样本尽可能的逼近于 p(x)。 

首先我们需要找到一个常数 k 使得 kq(x) 一定大于等于 p(x), 也就是如图所示（图摘自 [2]），p(x) 在 q(x) 下面。接着对 q(x) 进行采样，假设得到的样本为 。然后我们按照均匀分布在  中采样得到 u。如果 u 落到了图中的灰色区域，则拒绝这次采样，否则接受样本 。重复这个过程得到一系列的样本 。

在这个过程中，我们可以发现只有当建议分布 q(x) 和复杂分布 p(x) 重合的尽可能多的地方的样本更有可能被接受。那么相反的，如果他们重合部分非常少，那么就会导致拒绝的比例很高，抽样的效率就降低了。很多时候我们时候在高维空间进行采样，所以即使 p(x) 和 q(x) 很接近，两者的差异也可能很大。 

我们可以发现接受-拒绝法需要我们提出一个建议分布和常量，而且采样的效率也不高，那么我们就需要一个更一般的方法，就是我们的马尔科夫蒙特卡洛法。不过 MCMC 依旧有个问题, 它的抽样样本不是独立的。到了 Gibbs Sampling 的部分，我们可以看到它做了一些小 trick，来假装使得样本独立。

马尔科夫链

马尔科夫链的定义和概念相信大家都很熟悉了（不熟悉的请 google）。在这里，我们回顾一下几个重要的知识点：

平稳分布 

设有马尔科夫链 , 其状态空间为 S，转移概率矩阵为