【详解+推导！！】蒙特卡洛方法、接受拒绝采样、重要性采样、MCMC方法

蒙特卡洛方法、接受拒绝采样、重要性采样、MCMC方法

一、蒙特卡洛方法

1. 什么是蒙特卡洛方法

蒙特卡罗方法也称统计模拟方法，以概率为基础的方法，与它对应的是确定性算法，其核心思想是使用采样+平均(随机近似)的方式去估计出无法计算的值。

蒙特卡洛方法的一个重要应用就是求定积分。来看下面的一个例子。

2. 蒙特卡洛积分推导

我们先说明一般情况再举具体的例子。假设要计算的积分如下：

$$I={\int }_a^{b}g(x)dx$$
其中被积函数$g(x)$在$[a,b]$内可积,我们选择一个概率密度函数为$f_X(x)$的方式进行采样，并且保证在$[a,b]$内该采样函数的概率和为1，即${\int }_a^b{f}_X(x)dx=1$，那么原积分可写成：
$$\begin{array}{rl}I& ={\int }_a^b\frac{g(x)}{f_X(x)}{f}_X(x)dx\\ & ={\int }_a^b\frac{g(x)}{f_X(x)}dx\cdot {\int }_a^b{f}_X(x)dx\\ & ={\int }_a^b\frac{g(x)}{f_X(x)}dx\end{array}$$
我们以概率密度函数$f_X(x)$的方式进行采样，得到$N$个样本，则可以将上面的式子写成离散形式：
$$I^{\prime }=\frac{1}{N}\sum _{i=1}^N\frac{g(x_i)}{f_X(x_i)}$$
则有 $I^{\prime }\approx I$

举例子：

上面的例子其实就是取$f_X$为均匀分布：
$$f_X(x)=\frac{1}{b-a},a\le x\le b$$
带入我们最终推导出来的式子即：
$$I^{\prime }=\frac{b-a}{N}\sum _{i=1}^{N}g(x_i)$$

3.python实例

首先看一个经典的用蒙特卡洛方法求$\pi$值。

import random

def calpai():
    n = 1000000
    r = 1.0
    a, b = (0.0, 0.0)
    x_neg, x_pos = a - r, a + r
    y_neg, y_pos = b - r, b + r

    count = 0
    for i in range(0, n):
        x = random.uniform(x_neg, x_pos)
        y = random.uniform(y_neg, y_pos)
        if x*x + y*y <= 1.0:
            count += 1

    print (count / float(n)) * 4

正方形内部有一个相切的圆，它们的面积之比是$\pi /4$。现在，在这个正方形内部，随机产生$n$个点，计算它们与中心点的距离，并且判断是否落在圆的内部。若这些点均匀分布，则圆周率 $\pi =4\ast count/n$, 其中$count$表示落到圆内投点数 $n$表示总的投点数。

承上启下

从上面的介绍中我们可以看出蒙特卡洛方法非常直观，但是实际中样本的概率分布$f_X$[ 后面我们以$p(x)$表示样本点的概率分布 ]往往十分复杂，我们很难采样出符合这个概率分布的一组样本。

为了解决概率分布$p(x)$复杂，难以采样的问题派生出了两种解决方案，分别为：接受-拒绝采样和重要性采样。这两种方法都是基于一个事实，目标分布$p(x)$无法直接采样，于是引入一个提议分布$q(x)$来辅助采样。所谓提议分布就是一个我们已知的，很好采样的分布，例如：均分分布、高斯分布等。

二、接受-拒绝采样

1. 核心思想介绍

这里我们介绍最基本的接受-拒绝采样方法，目标是根据给定的概率密度函数$p(x)$产生服从该分布的样本集$X$。

基本思想是利用一个易于采样的提议分布$q(x)$，例如高斯分布，进行采样，而后再对利用$q(x)$得到的样本进行加工，接受符合$p(x)$分布的样本，拒绝不符合$p(x)$分布的样本。

2. 例子

首先我们已知一个复杂分布$p(x)$；再取一个提议分布$q(x)$(这里取高斯分分布，图中的变量为$z$，我文字描述中为$x$)；随后确定一个常量$k$，保证满足$kq(x)\ge p(x)$，在保证上式成立的情况下常量$k$尽可能取小一些。绘制出图像如下图所示，蓝色为$kq(x)$，红色为$p(x)$。

采样过程：

① 首先利用$q(x)$采样得到一个样本点$x_0$；

② 再在$[0,kq(x_0)]$区间内进行均匀采样，得到${\mu }_0$；

③ 判断：若${\mu }_0<p(x)$，则接受该样本，否则拒绝。

采样结果分布图：

由上图可知，该方法可以采样出已知概率分布的样本集。

三、重要性采样

1. 普通重要性采样

回顾一下：我们已知$f(x)$，其中$x～p(x)$，我们计算$f(x)$的期望如下：
$$E[f(x)]=\int f(x)p(x)dx\approx \frac{1}{N}\sum _{i=1}^{N}f(x_i)$$
上式中的$x_i$按照$x～p(x)$的方式进行采样，即可计算出$E[f(x)]$的近似值。

如果$p(x)$比较复杂无法完成采样，这时候我们借助一个简单的提议分布$q(x)$来完成采样，则期望的计算变成如下的样子：
$$E[f(x)]=\int f(x)\frac{p(x)}{q(x)}q(x)dx\approx \frac{1}{N}\sum _{i=1}^{N}f(x_i)\frac{p(x_i)}{q(x_i)}$$
上面式子中的$p(x)$和$q(x)$都是已知的，所以便可以利用这个公式计算出$f(x)$的期望，并且由推导可知该积分估计是无偏估计。

我们定义重要性权重：${\omega }_i=\frac{p(x_i)}{q(x_i)}$，则可以写出普通重要性权重的积分估计方程：
$$E[f(x)]\approx \frac{1}{N}\sum _{i=1}^N{\omega }_{i}f(x_i)$$

2. 方差分析

所谓无偏估计是指重要性采样计算得到的期望估计值等于真实的期望，这一点从推导过程中可以看出来，这也是我们可以用重要性采样的依据。下面我们来分析一下方差，方差的计算公式为：
$$D(x)=E(x^2)-E(x)$$
这里的$x$是对应的是$\frac{p(x)}{q(x)}f(x)$，由此可以看出如果$\frac{p(x)}{q(x)}$很大，那么方差也会很大，理论上方差可以达到无穷大。所以当$p(x)$和$q(x)$差异很大的时候，方差会变大，二者比较接近的时候，方差会比较接近真正的方差。下面我们来举个实际的例子。

我们定义：$f(x)=\frac{1}{1+exp(-x)}$，其函数图像如下：

• 实验一：$p(x)$和$q(x)$比较接近

  

  

  

• 实验二：$p(x)$和$q(x)$差异较大

  

由此我们可以看出来，如果$p(x)$和$q(x)$差异较大，我们就需要做更多的采样才能获得比较准确的期望估计，并且会导致估计的方差较大。

3. 加权重要性采样

加权重要性采样是一种减少重要性采样积分方差的方法，公式如下：
$$E[f]\approx \sum _{n=1}^N\frac{{\omega }^n}{{\sum }_{m=1}^N{\omega }^m}f(x^n)$$

四、MCMC方法

MCMC是指马尔科夫链蒙特卡洛(Markov Chain Monte Carlo)方法，被视为二十世纪Top10的算法，正是MCMC方法的提出使得许多贝叶斯统计问题的求解成为可能。

当采样空间的维数比较高时，“接受-拒绝采样”和“重要性采样”都不适用，原因是当维数升高时很难找到适合的提议分布，采样效率比较差。而MCMC算法不需要提议分布，只需要一个初始样本点，后面的样本点会依据当前样本点随机产生。

这里需要马尔科夫链中的一些术语和知识可以参考：https://blog.csdn.net/qq_33302004/article/details/115027798

本节会先介绍一些定理和知识，最后再讲解MCMC的采样过程，如果不能理解基础定理的话很难理解MCMC方法，大家耐心理解，也可以看一下后面我的引用文章。

1. 马氏链平稳分布

先举例子再讲公式：

由此我们可以看出马尔科夫链最终的稳定分布和初始状态无关，只与状态转移概率矩阵$P$有关。

我们的目标分布$\pi$这个就是最终的稳定分布，专业一点表述就是我们本小节的小标题马氏链平稳分布。我们也可以知道，目标分布$\pi$是由状态转移概率矩阵$P$决定的。

用数学方式描述马氏链平稳分布：

其背后蕴含的意义：

• 任意两个状态连通并非指状态i可以一步转移到j，而是指状态i可以通过有限的次数转移到j。
• 状态转移矩阵自乘多次转移后，会得到一个稳定值。
• 马氏链稳定后，所有状态i转移到状态j的概率之和是稳定的。
• 一个状态转移矩阵$P$只有唯一一个对应的稳定分布$\pi (x)$

2. 平稳分布判定：细致平稳条件

假设我们已知一个状态概率分布$\pi$，如何判定它是否是一个状态转移概率矩阵$P$的稳定分布呢？这里就用到了细致平稳条件。

如果状态转移矩阵$P$和分布$\pi$满足细致平稳条件，则这个分布$\pi$就是该状态转移矩阵$P$对应的稳定分布。细致平稳条件定义如下：
$$\pi (i)p_{ij}=\pi (j)p_{ji},\text{ for all i, j}$$

证明如下：

3. 构造状态转移概率矩阵$P$

我们实际应用中的场景是这样的：已知目标分布$\pi$，完成按照目标分布的采样。所以我们的目的是构造出这个目标分布对应的状态转移概率矩阵$P$，这样我们就可以从一个初始状态不断使用MCMC进行采样了。

我们能够利用的工具就是上面介绍的细致平稳条件。

我们任取一个状态转移矩阵$Q$，一般这个任取的$Q$是步满足细致平稳条件的：
$$p(i)q(i,j)\ne p(j)q(j,i)$$
我们可以构造一个$\alpha$使上式满足细致平稳条件：
$$p(i)q(i,j)\alpha (i,j)=p(j)q(j,i)\alpha (j,i)$$
其中：
$$\alpha (i,j)=p(j)q(j,i)，\alpha (j,i)=p(i)q(i,j)$$
这样我们就构造出了目标马尔科夫状态转移矩阵，其中每一项都满足：
$$P(i,j)=Q(i,j)\alpha (i,j)$$
继承“接受-拒绝采样”和“重要性采样”中的说法，我们可以把$Q$称作提议马尔科夫链，其中状态$i$转移到状态 $j$的概率是$Q(i,j)$；而我们的目标马尔科夫链中，从状态$i$转移到状态 $j$的概率是$Q(i,j)\alpha (i,j)$。

继续继承“接受-拒绝采样”中的思想，我们可以很容易的把$\alpha (i,j)$看做接受概率或者接受率。

4. MCMC的算法步骤

紧接上一小节，$\alpha (i,j)$表示的是采样样本的接受概率，所以MCMC的核心思想就是按照提议马尔科夫链$Q$进行采样，而后按照$\alpha (i,j)$确定采样是否被接受。MCMC的主要步骤如下：

1. 初始化马尔科夫链的初始状态 $X_0=x_0$;
2. 循环执行以下步骤（3~6），迭代若干次，$t=0,1,2,...$
3. 在$t$时刻状态$X_t=x_t$，按照$Q$进行采样，即$y～q(x|x_t)$
4. 从均匀分布中进行采样$u～Uniform[0,1]$
5. 若 $u<\alpha (x_t,y)=p(y)q(x_t|y)$ 则接受该采样，状态$X_{t+1}=y$；否则不接受转移$X_{t+1}=x_t$

5. Metropolis-Hastings采样方法

M-H采样方法是对MCMC方法的一点改进，改进之后采样的接受度升高，加快了采样效率。

从第3小结中我们知道
$$\alpha (i,j)=p(j)q(j,i)$$

所以通常$\alpha (i,j)$是一个很小的值，因此在一次采样中很容易被拒绝，从而导致采样的时间成本、计算成本很高。M-H(Metropolis-Hastings)采样，通过放大接受率$\alpha (i,j)$提高了跳转概率，缩短了采样过程，推导过程如下：
$$\begin{array}{rl}p(i)q(i,j)p(j)q(j,i)& =p(j)q(j,i)p(i)q(i,j)\\ p(i)q(i,j)\frac{p(j)q(j,i)}{p(i)q(i,j)}& =p(j)q(j,i)\end{array}$$
有：
$$\alpha (i,j)=\frac{p(j)q(j,i)}{p(i)q(i,j)}，\alpha (j,i)=1$$
从公式中我们可以看出，这样计算的$\alpha (i,j)$是可能大于1的，所以我们定义：
$$\alpha (i,j)=min(\frac{p(j)q(j,i)}{p(i)q(i,j)},1)$$
这样我们就将接受概率放大了$\frac{1}{p(i)q(i,j)}$倍，最多可以实现满概率跳转。

最后推荐几篇相关文章：

蒙特卡洛方法：https://blog.csdn.net/bitcarmanlee/article/details/82716641
接受拒绝采样：https://www.zhihu.com/question/38056285/answer/1803920100
重要性采样：https://mathpretty.com/12375.html
马尔科夫链的平稳分布：http://www.360doc.com/content/20/0313/14/10724725_898892977.shtml
细致平稳条件、Metropolis-Hastings采样方法：https://www.zhihu.com/question/63305712/answer/1804780073
MCMC方法：https://blog.csdn.net/arcers/article/details/88732639