机器学习：Gibbs 采样

1. 什么是Gibbs采样

Gibbs采样是MH算法的一种特例（α==1)，因此可以保证Gibbs抽取的样本，也构成一个非周期不可约稳定收敛的马氏链；Gibbs采样适用于样本是两维或以上的情况；通过积分去除掉相关但是不感兴趣的变量，称为“collapsed”的Gibbs采样；并且个人的一个感觉是，观测量所直接依赖的那些变量是不能被积分掉的，否则无法有效的进行抽样…

gibbs采样需要知道样本中一个属性在其它所有属性下的条件概率，然后利用这个条件概率来分布产生各个属性的样本值。gibbs采样属于随机模拟抽样算法中的一种（一类近似求解的方法）。随机模拟的核心是对一个分布进行抽样，常用的抽样算法包括：1. 接受-拒绝抽样；2）重要性抽样；3）MCMC（马尔科夫链蒙特卡洛方法）方法，它包括两个非常著名的采样算法（metropolis-hasting算法和它的特例Gibbs采样算法）（补充：MCMC方法最早由Metropolis（1954）给出，后来Metropolis的算法由Hastings改进，合称为M-H算法。M-H算法是MCMC的基础方法。由M-H算法演化出了许多新的抽样方法，包括目前在MCMC中最常用的Gibbs抽样也可以看做M-H算法的一个特例）。

Gibbs算法，就是用条件分布的抽样来替代全概率分布的抽样。例如，X={x1,x2,...xn}满足分布p(X)，如何对p(X)进行抽样呢？如果我们知道它的条件分布p(x1|X_{-1}),...,p(xi|X_{-i}),....，其中X_{-i}表示除了xi之外X的所有变量。如果这些条件分布都是很容易抽样的，那么我们就可以通过对条件分布的抽样来对全概率分布p(X)进行抽样。

Gibbs采样算法的步骤：

1. 给定一个初始样本X0={x10,x20,...,xn0}

2.已知一个样本Xi={x1i,x2i,...,xni}，对于x1_{i+1}进行抽样，x1_{i+1} ~ p(x1|Xi_{-1})

3. 对于x2_{i+1}进行抽样，x2_{i+1} ~ p(x2|x1_{i+1}, x3i,...xni)

4.对于xn_{i+1}进行抽样，xn_{i+1} ~ p(xn|x1_{i+1}, x2_{i+1},...x_{n-1}_{i+1})

5.步骤2~4可以得到X的一个样本，然后重复步骤2~4可以不断地得到X的样本。

当然无论是metropolis-hasting算法还是gibbs算法，都有一个burn in的过程，所谓burn in的过程就是因为这个两个算法本身都是markov chain的算法，要达到稳定状态需要一定的步骤才能达到，所以需要一个burn in过程，只有在达到平衡状态时候得到的样本才能是平衡状态时候的目标分布的样本，因此，在burn in过程中产生的样本都需要被舍弃。如何判断一个过程是否达到了平衡状态还没有一个成熟的方法来解决，目前常见的方法是看是否状态已经平稳（例如画一个图，如果在较长的过程中，变化已经不大，说明很有可能已经平衡）当然这个方法并不能肯定一个状态是否平衡，你可以举出反例，但是却是实际中没有办法的办法。Gibbs采样的目的是获得一个样本，不是计算概率，但可以通过其他方法来统计概率。

2. 常见的采样方法

2.1 直接采样（简单）

2.2 接受-拒绝抽样（Acceptance-Rejection sampling)下面内容来源

很多实际问题中，p(x)是很难直接采样的的，因此，我们需要求助其他的手段来采样。既然 p(x) 太复杂在程序中没法直接采样，那么我设定一个程序可抽样的分布 q(x) 比如高斯分布，然后按照一定的方法拒绝某些样本，达到接近 p(x) 分布的目的，其中q(x)叫做 proposal distribution（建议分布）。

具体操作如下，设定一个方便抽样的函数 q(x)，以及一个常量 k，使得 p(x) 总在 kq(x) 的下方。

1）.x 轴方向：从 q(x) 分布抽样得到 a。(如果是高斯，就用之前说过的 tricky and faster 的算法更快）

2）.y 轴方向：从均匀分布（0, kq(a)) 中抽样得到 u。

3）.如果刚好落到灰色区域： u > p(a), 拒绝，否则接受这次抽样

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

3. 马尔科夫链，马尔科夫稳态

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

也就是说前一个状态只与当前状态有关，而与其他状态无关，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)。如何能做到这一点呢？我们主要使用如下的定理。