马尔可夫链蒙特卡罗算法 MCMC_mcmc算法

马尔可夫链蒙特卡罗算法（MCMC）是贝叶斯推断中的明星算法，困惑笔者颇久，阅读了刘建平大佬的博客及知乎上的一些阅读笔记后，终于有了些自己的理解。本文基于刘建平大佬的博客进行梳理，复制粘贴较多（甚至截图），权且当成读书笔记。

蒙特卡罗方法

原文链接：MCMC(一)蒙特卡罗方法

要理解好MCMC，得先从第二个MC——Monte Carlo开始说起。

引入

求解积分问题（即面积计算类问题，尤其是难以直接数学推导时），采用Monte Carlo算法是个不错的选择。比如计算圆周率，可以在边长为2的正方形内画一个内切圆，然后再正方形内进行N次随机采样，求出圆内样本点的个数所占比例即可近似求出圆面积，从而间接推导圆周率的值。再比如下图，求蝙蝠侠的大小，也可以用类似的思想求解。这便是蒙特卡洛（Monte Carlo）的基本哲学。

      

当然，实际情况运用这个思想，得稍微灵活一点。以下用原文的例子一步步贴近实际需求进行说明。

问题：求解如下函数的积分：

        

      

这个例子很好，只是笔者在第一次读的时候没理解为什么求积分需要用到概率分布，总感觉均匀采样就行了，使用概率分布的必要性在哪？仔细一想，发现其实这才是现实的需求，这种做法很有必要！以吃饭为例（暴露了吃货的本质）：f(x)表示学生饭量随身高的变化趋势（曲线应当是逐渐上升），如果求食堂需要准备多少米饭，此时我们还有必要了解学生的身高分布情况（频数分布）；如果求单个学生饭量的期望值，则要用身高的概率分布。这样类似的例子还有很多很多。

接受-拒绝采样

如果已知概率分布函数，如何基于这个概率分布进行采样呢？

我们知道，对于正态分布、Beta分布、Gamma分布等等经典分布，程序包都有相应的随机采样方法可以直接调用。其实这些分布之所以能直接使用，是基于一些变换得到的（比如正态分布是在均匀分布的基础上经过Box−Muller变换得到）。

但是对于自定义的复杂概率分布，一般无法直接进行随机采样。这时，我们采取“接受-拒绝采样”的策略。

      

总的来说，就是我们通过怎样的采样策略，让获得的样本的分布形态与概率分布函数的形态接近。如何接地气地理解上述策略，笔者思考如下。最自然的想法是，我们找一个矩形（定义域两端开发时 矩形无限长），让分布函数被包含在这个矩形内，然后在矩形内均匀采样，若样本点落在了目标概率函数下方，则保留该样本。这其实就对应着上述方法中的proposal distribution为均匀分布的情况。

接下来，我们考虑优化这个自然的方法。矩形的高度和长度如何确定？选择最小的能覆盖目标概率函数的矩形即可！如果目标概率函数的两端无限延长，中间高耸，类似于正态分布呢？此时矩形两端无限延长，均匀采样的效率极低。将proposal distribution改成正态分布，采样效率明显提升。实际情况采用怎样的proposal distribution，需要根据目标概率函数的特点而定。有时均匀分布不均合适，有时状态分布合适，有时Beta分布更合适；合适的proposal distribution可以带来更高的采用效率。

      

从上面可以看出，要想将蒙特卡罗方法作为一个通用的采样模拟求和的方法，必须解决如何方便得到各种复杂概率分布的对应的采样样本集的问题。接下来要讲到的马尔科夫链就是帮助找到这些复杂概率分布的对应的采样样本集的白衣骑士。

 

马尔科夫链

原文链接：MCMC(二)马尔科夫链

概述

马尔科夫链（Markov chain）：假设某一时刻状态转移的概率只依赖于它的前一个状态。举个形象的比喻，假如每天的天气是一个状态的话，那个今天是不是晴天只依赖于昨天的天气，而和前天的天气没有任何关系。当然这么说可能有些武断，但是这样做可以大大简化模型的复杂度，因此马尔科夫链在很多时间序列模型中得到广泛的应用，比如循环神经网络RNN，隐式马尔科夫模型HMM等，当然MCMC也需要它。

      

注意：状态转移矩阵P的每一行的和为1。在连续概率分布中，马尔科夫链的状态数是无限的，也就是说，每一行的元素个数有无限个。马尔科夫链转移（以AB两个状态为例），不是指A状态转移到B状态，而是：当前处于A或B状态的概率，转移到下一个阶段（时刻）后，处于A或B状态的概率分别变成了多少。

马尔科夫链模型的状态转移矩阵的性质

原文中，作者通过模拟，引出了马尔科夫链的收敛性质，此处直接抛出结论（模拟及推导此处略去）。

马尔科夫链的收敛性质

      

      

 

基于马尔科夫链采样

      

马尔科夫链和蒙特卡洛是怎么联系起来的？先重新回顾之前蒙特卡洛采样：给定概率分布函数，如何获得有代表性的采样点。在初始的几个采样点，形成的概率分布可能不符合概率分布函数，但随着采样次数逐渐增加，越来越接近概率分布函数形态，并且分布的形态将达到稳态，这对也就应于马尔科夫链的稳态（马尔科夫链的节点是当前阶段的概率分布）。我们在达到平稳分布后进行采样，这些样本就可以用来进行蒙特卡洛模拟了。也就是说，我们已知马尔科夫链的稳态，需要重现马尔科夫链的状态转移过程，也就是需要状态转移矩阵（初始状态并不重要）；状态转移矩阵如何求得，接下来的内容将进行阐述。

马尔科夫链的细致平稳条件

      

MCMC采样

原文链接：MCMC(三)MCMC采样和M-H采样

      

      

仔细观察第3c步，这里就和朴素的蒙特卡洛方法很像了。由此理解接受率，其实这个接受率和朴素蒙特卡洛中的接受率是类似的意思，都是动态确定的，根据当前概率分布（此处的，蒙特卡洛中的p分布）及使用的辅助分布（此处的Q，蒙特卡洛中的q分布）进行计算确定。总的来说，上述MCMC采样过程，就是Markov chain采样过程和Monte Carlo采样过程的综合版。

上面这个过程基本上就是MCMC采样的完整采样理论了，但是这个采样算法还是比较难在实际中应用，为什么呢？问题在上面第三步的c步骤，接受率这儿。由于 可能非常的小，比如0.1，导致我们大部分的采样值都被拒绝转移，采样效率很低。有可能我们采样了上百万次马尔可夫链还没有收敛，也就是上面这个n1要非常非常的大，这让人难以接受，怎么办呢？这时就轮到我们的M-H采样出场了。

M-H采样

M-H采样是Metropolis-Hastings采样的简称，这个算法首先由Metropolis提出，被Hastings改进，因此被称之为Metropolis-Hastings采样或M-H采样。M-H采样解决了我们上一节MCMC采样接受率过低的问题。

      

 

很多时候，我们选择的马尔科夫链状态转移矩阵QQ如果是对称的，即满足Q(i,j)=Q(j,i),这时我们的接受率可以进一步简化为：

      

      

Gibbs采样解决了上面两个问题，因此在大数据时代，MCMC采样基本是Gibbs采样的天下，接下来我们就来讨论Gibbs采样。

Gibbs采样

原文链接：MCMC(四)Gibbs采样

Gibbs采样系列，笔者尚未理解透，暂时把原文截图存于此。

重新寻找合适的细致平稳条件

在M-H采样中我们通过引入接受率使细致平稳条件满足。现在我们换一个思路。

      

二维Gibbs采样

      

      

多维Gibbs采样

      

       由于Gibbs采样在高维特征时的优势，目前我们通常意义上的MCMC采样都是用的Gibbs采样。当然Gibbs采样是从M-H采样的基础上的进化而来的，同时Gibbs采样要求数据至少有两个维度，一维概率分布的采样是没法用Gibbs采样的,这时M-H采样仍然成立。

　　有了Gibbs采样来获取概率分布的样本集，有了蒙特卡罗方法来用样本集模拟求和，他们一起就奠定了MCMC算法在大数据时代高维数据模拟求和时的作用。

 

参考资料：

MCMC(一)蒙特卡罗方法

MCMC(二)马尔科夫链

MCMC(三)MCMC采样和M-H采样

MCMC(四)Gibbs采样

马尔可夫链蒙特卡罗算法（MCMC）

告别数学公式，图文解读什么是马尔可夫链蒙特卡罗方法（MCMC）