随机采样方法整理与讲解（MCMC、Gibbs Sampling等）

本文是对参考资料中多篇关于sampling的内容进行总结+搬运，方便之后本身翻阅。其实参考资料中的资料写的比我好，你们能够看一下！好东西多分享！PRML的第11章也是sampling，有时间后面写到PRML的笔记中去：)

背景

随机模拟也能够叫作蒙特卡罗模拟(Monte Carlo Simulation)。这个方法的发展始于20世纪40年代，和原子弹制造的曼哈顿计划密切相关，当时的几个大牛，包括乌拉姆、冯.诺依曼、费米、费曼、Nicholas Metropolis， 在美国洛斯阿拉莫斯国家实验室研究裂变物质的中子连锁反应的时候，开始使用统计模拟的方法,并在最先的计算机上进行编程实现。[3]

随机模拟中有一个重要的问题就是给定一个几率分布p(x)，咱们如何在计算机中生成它的样本。通常而言均匀分布 Uniform(0,1)的样本是相对容易生成的。 经过线性同余发生器能够生成伪随机数，咱们用肯定性算法生成[0,1]之间的伪随机数序列后，这些序列的各类统计指标和均匀分布 Uniform(0,1) 的理论计算结果很是接近。这样的伪随机序列就有比较好的统计性质，能够被当成真实的随机数使用。

 

下面总结这么几点：

一、蒙特卡洛数值积分

二、均匀分布，Box-Muller 变换

三、Monte Carlo principle

四、接受-拒绝抽样（Acceptance-Rejection sampling)

五、重要性抽样(Importance sampling)

六、马尔科夫链，马尔科夫稳态

七、MCMC——Metropolis-Hasting算法

八、MCMC——Gibbs Sampling算法

 

一、蒙特卡洛数值积分

若是咱们要求f(x)的积分，如

而f(x)的形式比较复杂积分很差求，则能够经过数值解法来求近似的结果。经常使用的方法是蒙特卡洛积分：

这样把q(x)看作是x在区间内的几率分布，而把前面的分数部门看作一个函数，而后在q(x)下抽取n个样本，当n足够大时，能够用采用均值来近似：

所以只要q（x）比较容易采到数据样本就好了。随机模拟方法的核心就是如何对一个几率分布获得样本，即抽样（sampling）。下面咱们将介绍经常使用的抽样方法。

 

二、均匀分布，Box-Muller 变换

在计算机中生成[0,1]之间的伪随机数序列，就能够当作是一种均匀分布。而随机数生成方法有不少，最简单的如：

固然计算机产生的随机数都是伪随机数，不过通常也就够用了。

 

[Box-Muller 变换]  若是随机变量 U1,U2 独立且U1,U2∼Uniform[0,1]，

则 Z0,Z1 独立且服从标准正态分布。

 

三、Monte Carlo principle

Monte Carlo 抽样计算随即变量的指望值是接下来内容的重点：X 表示随即变量，服从几率分布 p(x), 那么要计算 f(x) 的指望，只须要咱们不停从 p(x) 中抽样xi，而后对这些f（xi）取平均便可近似f(x)的指望。

 

 

四、接受-拒绝抽样（Acceptance-Rejection sampling)[2]

不少实际问题中，p(x)是很难直接采样的的，所以，咱们须要求助其余的手段来采样。既然 p(x) 太复杂在程序中无法直接采样，那么我设定一个程序可抽样的分布 q(x) 好比高斯分布，而后按照必定的方法拒绝某些样本，达到接近 p(x) 分布的目的，其中q(x)叫作 proposal distribution 。

具体操做以下，设定一个方便抽样的函数 q(x)，以及一个常量 k，使得 p(x) 总在 kq(x) 的下方。（参考上图）

• x 轴方向：从 q(x) 分布抽样获得 a。(若是是高斯，就用以前说过的 tricky and faster 的算法更快）
• y 轴方向：从均匀分布（0, kq(a)) 中抽样获得 u。
• 若是恰好落到灰色区域： u > p(a), 拒绝， 不然接受此次抽样
• 重复以上过程

在高维的状况下，Rejection Sampling 会出现两个问题，第一是合适的 q 分布比较难以找到，第二是很难肯定一个合理的 k 值。这两个问题会致使拒绝率很高，无用计算增长。

 

五、重要性抽样(Importance sampling)[2]

Importance Sampling 也是借助了容易抽样的分布 q (proposal distribution)来解决这个问题，直接从公式出发：

其中，p(z) / q(z) 能够看作 importance weight。咱们来考察一下上面的式子，p 和 f 是肯定的，咱们要肯定的是 q。要肯定一个什么样的分布才会让采样的效果比较好呢？直观的感受是，样本的方差越小指望收敛速率越快。好比一次采样是 0, 一次采样是 1000, 平均值是 500,这样采样效果不好，若是一次采样是 499, 一次采样是 501, 你说指望是 500,可信度还比较高。在上式中，咱们目标是 p×f/q 方差越小越好，因此 |p×f| 大的地方，proposal distribution q(z) 也应该大。举个稍微极端的例子：

第一个图表示 p 分布， 第二个图的阴影区域 f = 1，非阴影区域 f = 0, 那么一个良好的 q 分布应该在左边箭头所指的区域有很高的分布几率，由于在其余区域的采样计算实际上都是无效的。这代表 Importance Sampling 有可能比用原来的 p 分布抽样更加有效。

可是惋惜的是，在高维空间里找到一个这样合适的 q 很是难。即便有 Adaptive importance sampling 和 Sampling-Importance-Resampling(SIR) 的出现，要找到一个同时知足 easy to sample 而且 good approximations 的 proposal distribution, it is often impossible！

 

六、马尔科夫链，马尔科夫稳态

在讲蒙特卡洛方法以前，必需要先讲一下马尔科夫链；马氏链的数学定义：

也就是说前一个状态只与当前状态有关，而与其余状态无关，Markov Chain 体现的是状态空间的转换关系，下一个状态只决定与当前的状态(能够联想网页爬虫原理，根据当前页面的超连接访问下一个网页)。以下图：

举一个例子，若是当前状态为 u(x) = (0.5, 0.2, 0.3), 那么下一个矩阵的状态就是 u(x)T = (0.18, 0.64, 0.18), 依照这个转换矩阵一直转换下去，最后的系统就趋近于一个稳定状态 (0.22, 0.41, 0.37) (此处只保留了两位有效数字)。而事实证实不管你从那个点出发，通过很长的 Markov Chain 以后都会聚集到这一点。[2]

再举一个例子，社会学家常常把人按其经济情况分红3类：下层(lower-class)、中层(middle-class)、上层(upper-class)，咱们用1,2,3 分别表明这三个阶层。社会学家们发现决定一我的的收入阶层的最重要的因素就是其父母的收入阶层。若是一我的的收入属于下层类别，那么他的孩子属于下层收入的几率是 0.65, 属于中层收入的几率是 0.28, 属于上层收入的几率是 0.07。事实上，从父代到子代，收入阶层的变化的转移几率以下

使用矩阵的表示方式，转移几率矩阵记为

咱们发现从第7代人开始，这个分布就稳定不变了，事实上，在这个问题中，从任意初始几率分布开始都会收敛到这个上面这个稳定的结果。

注：要求图是联通的（没有孤立点），同时不存在一个联通的子图是没有对外的出边的（就像黑洞同样）。

这个马氏链的收敛定理很是重要，全部的 MCMC(Markov Chain Monte Carlo) 方法都是以这个定理做为理论基础的。

对于给定的几率分布p(x),咱们但愿能有便捷的方式生成它对应的样本。因为马氏链能收敛到平稳分布， 因而一个很的漂亮想法是：若是咱们能构造一个转移矩阵为P的马氏链，使得该马氏链的平稳分布刚好是p(x), 那么咱们从任何一个初始状态x0出发沿着马氏链转移, 获得一个转移序列 x0,x1,x2,⋯xn,xn+1⋯,， 若是马氏链在第n步已经收敛了，因而咱们就获得了 π(x) 的样本xn,xn+1⋯。

这个绝妙的想法在1953年被 Metropolis想到了，为了研究粒子系统的平稳性质， Metropolis 考虑了物理学中常见的波尔兹曼分布的采样问题，首次提出了基于马氏链的蒙特卡罗方法，即Metropolis算法，并在最先的计算机上编程实现。Metropolis 算法是首个普适的采样方法，并启发了一系列 MCMC方法，因此人们把它视为随机模拟技术腾飞的起点。 Metropolis的这篇论文被收录在《统计学中的重大突破》中， Metropolis算法也被遴选为二十世纪的十个最重要的算法之一。

咱们接下来介绍的MCMC 算法是 Metropolis 算法的一个改进变种，即经常使用的 Metropolis-Hastings 算法。由上一节的例子和定理咱们看到了，马氏链的收敛性质主要由转移矩阵P 决定, 因此基于马氏链作采样的关键问题是如何构造转移矩阵P,使得平稳分布刚好是咱们要的分布p(x)。如何能作到这一点呢？咱们主要使用以下的定理。

 

 

                                    马氏链转移和接受几率

 

假设咱们已经有一个转移矩阵Q(对应元素为q(i,j)), 把以上的过程整理一下，咱们就获得了以下的用于采样几率分布p(x)的算法。

 

 

八、MCMC——Gibbs Sampling算法

          平面上马氏链转移矩阵的构造

 

 

以上算法收敛后，获得的就是几率分布p(x1,x2,⋯,xn)的样本，固然这些样本并不独立，可是咱们此处要求的是采样获得的样本符合给定的几率分布，并不要求独立。一样的，在以上算法中，坐标轴轮换采样不是必须的，能够在坐标轴轮换中引入随机性，这时候转移矩阵 Q 中任何两个点的转移几率中就会包含坐标轴选择的几率，而在一般的 Gibbs Sampling 算法中，坐标轴轮换是一个肯定性的过程，也就是在给定时刻t，在一根固定的坐标轴上转移的几率是1。

 

参考资料

[1] http://blog.csdn.net/xianlingmao/article/details/7768833

[2] http://www.cnblogs.com/daniel-D/p/3388724.html

[3] http://cos.name/2013/01/lda-math-mcmc-and-gibbs-sampling/

[4] An Introduction to MCMC for Machine Learning，2003

[5] Introduction to Monte Carlo Methods