常用采样方法

常用采样方法

最近在学习 MCMC，一种特殊的采样方法，顺便把其他常用的方法了解了一下。

为什么要采样？

很多问题，我们只需要使用数学解析的方法即可解决。例如对 f(x)做积分，如果 f(x) = x^2，那么直接积分就行，很简单。

若f(x)是标准正态分布的概率密度函数（pdf），求[a,b]之间的定积分，那么直接用数学解析方法就搞不定了，因为我们知道正态分布的积分是不可求的。既然无法用解析方法计算精确值，那么退而求其次，不妨寻找一种可以求得近似值的方法。例如：

f(x) = f(x) * g(x) * 1/g(x) = g(x) * (f(x) * 1/g(x))。这里我们找到一个容易采样的 pdf g(x)，并对其采样得到{X1…Xn}。连续函数的积分就可以转化为离散值的求和。这被称为蒙特卡洛积分。

从上面例子可以看到，采样方法可以对诸如积分这种无法用解析方法求解的方法进行近似求解。接下来就是问题的核心：怎样采样？

简单采样方法

最简单的抽样是均匀采样，也就是均匀产生[0,1]之间的随机数，编程语言中一般使用 rand()函数。之所以说是均匀采样，是因为这些样本服从均匀分布。随机数的产生，在计算机中一般通过线性同余的方法实现。

另外，一些简单的分布，如正态分布，也可以在均匀采样的基础上实现采样。

对于其他函数我们就无能为力了，只能通过一些高级点的方式进行采样。

接受-拒绝采样

我们需要对一个分布π(x)进行采样，但是却很难直接进行采样，所以我们想通过另外一个容易采样的分布q(x)的样本，用某种机制去除掉一些样本，从而使得剩下的样本就是来自与所求分布f(x)的样本。

条件：

• 对q(x)采样比较容易
• q(x)的轮廓接近π(x)，且有 π(x)≤Mq(x),∀x

过程：

• 产生样本X~q(x)，和U~Uniform[0,1]
• Y=U*Mq(x)，若Y≤π(X)，则接受X，否则拒绝。

解释：

• 根据 q(x)采样 X，得到 Mq(X)，在[0,1]之间产生随机数U，也就是在[0,Mq(X)]产生随机数Y=U*Mq(x)。如果Y 在π(x)曲线下方，那么就选择接受，否则拒绝。为什么呢？在上图两个曲线相隔越远的地方，随机Y在π(x)下方的概率越小，即接受这个采样X的概率越小。这是合理的，这时候两个曲线间隔远， Mq(x)的采样X不能直接用于π(X)。反之，如果两个曲线相隔近，那么U越可能在π(X)下方，越可能接受这个采样，既然曲线相隔近，那么对Mq(x)的采样就可以近似对π(X)的采样了嘛！

但是拒绝采样效率会比较低，因为有很多采样都拒绝了嘛！

重要性采样

如图，如果我们想求积分，也就是面积，且 f(x)不能求积分形式，一种方法就是在[a,b]间均匀采样N 个点，并用 f(Xi)乘以(b-a)/N 即宽度，累加求和就能得到近似的积分值。这里我们用的权重是相同的：(b-a)/N，就是说每个小矩形的宽度都是相等的。

但很多时候，曲线比较高的地方需要多采样并精确刻画，曲线低的地方可以少采样，这样能减小最后结果与真实值之间的误差。如下图：

我们采用与 f(x)类似的 g(x)来采样，g(x)如图中右上角所示。此时宽度怎么确定呢？宽度就是1/g(Xi)，这就能体现出不同点的权重不同了。

MCMC采样

上述都是独立性采样，采样的效率还不高。MCMC 是一种关联采样，当前采样有赖于前一个采样结果。MCMC 的全程是马尔可夫蒙特卡罗。马尔可夫，是说前后两个采样结果的关联性。蒙特卡罗是一种随机模拟方法，用采样的方法解决解析问题，如前面提到的蒙特卡洛积分。

若想对π(x)进行采样，首先构建一个马尔可夫链，该马尔可夫链的状态转移矩阵满足特定条件时，会存在一个稳定状态，稳定分布就是π(x)。根据某定律，我们从某个状态出发，在马尔可夫链每步得到的分布中采样，得到 M 个样本，这些样本就近似服从π(x)。要想构建符合条件的状态转移矩阵并抽样，有两种方法：Metropolis-Hastings 算法和Gibbs Sampling方法。具体介绍可参考文末给出的链接。

参考：

http://www.wuhaijia.com/wordpress/%E6%8E%A5%E5%8F%97-%E6%8B%92%E7%BB%9D%E9%87%87%E6%A0%B7%EF%BC%88acceptance-rejection-sampling%EF%BC%89/

http://blog.sina.com.cn/s/blog_4e5740460100cw5b.html

http://www.52nlp.cn/lda-math-mcmc-%E5%92%8C-gibbs-sampling2