蒙特卡罗方法概述

前言

作为一种随机采样方法，马尔科夫链蒙特卡罗（Markov Chain Monte Carlo，以下简称MCMC）在机器学习,深度学习以及自然语言处理等领域都有广泛的应用，是很多复杂算法求解的基础。

从名字我们可以看出，MCMC由两个MC组成，即蒙特卡罗方法（Monte Carlo Simulation，简称MC）和马尔科夫链（Markov Chain ，也简称MC）。要弄懂MCMC的原理我们首先得搞清楚蒙特卡罗方法和马尔科夫链的原理。

蒙特卡罗方法

蒙特卡罗原来是一个赌场的名称，用它作为名字大概是因为蒙特卡罗方法是一种随机模拟的方法，这很像赌博场里面的扔骰子的过程。原理是通过大量随机样本，去了解一个系统，进而得到所要计算的值。

它非常强大和灵活，又相当简单易懂，很容易实现。对于许多问题来说，它往往是最简单的计算方法，有时甚至是唯一可行的方法。

具体的适合理解的例子参考这篇文章

给定统计样本集，如何估计产生这个样本集的随机变量概率密度函数,是我们比较熟悉的概率密度估计问题。
求解概率密度估计问题的常用方法是最大似然估计、最大后验估计等。但是，我们思考概率密度估计问题的逆问题:给定一个概率分布p(x)，如何让计算机生成满足这个概率分布的样本。
这个问题就是统计模拟中研究的重要问题–采样(Sampling)。

如果有一个我们很难求解出$f(x)$的原函数的函数, 现要求其在定义域 [a, b] 上的积分, 如果这个函数是均匀分布, 那么我们可以采样 [a,b] 区间的 n 个值：${x_0,x_1,…x_{n-1}}$,用它们的均值来代表 [a,b] 区间上所有的 f(x) 的值。这样我们上面的定积分的近似求解为:

$$\theta = \int_a^b f(x)dx = \frac{b-a}{n}\sum\limits_{i=0}^{n-1}f(x_i)$$

如果不是均匀分布, 并假设我们可以得到 $x$ 在$[a,b]$的概率分布函数$p(x)$，那么我们的定积分求和可以这样进行：

$$\theta = \int_a^b f(x)dx = \int_a^b \frac{f(x)}{p(x)}p(x)dx \approx \frac{1}{n}\sum\limits_{i=0}^{n-1}\frac{f(x_i)}{p(x_i)}$$

上式最右边的这个形式就是蒙特卡罗方法的一般形式。当然这里是连续函数形式的蒙特卡罗方法，但是在离散时一样成立。

可以看出，最上面我们假设x在[a,b]之间是均匀分布的时候，p(xi)=1/(b−a)，带入我们有概率分布的蒙特卡罗积分的上式，可以得到：

$$\frac{1}{n}\sum\limits_{i=0}^{n-1}\frac{f(x_i)}{1/(b-a)} = \frac{b-a}{n}\sum\limits_{i=0}^{n-1}f(x_i)$$

也就是说，我们最上面的均匀分布也可以作为一般概率分布函数$p(x)$在均匀分布时候的特例。那么我们现在的问题转到了如何在已知分布求出 x 的分布$p(x)$对应的若干个样本上来。

采样方法

主要有概率分布采样及接受-拒绝采样方法.

概率分布采样

如果求出了$x$的概率分布，我们可以基于概率分布去采样基于这个概率分布的 n 个$x$的样本集，带入蒙特卡罗求和的式子即可求解。但是还有一个关键的问题需要解决，即如何基于概率分布去采样基于这个概率分布的 n 个$x$的样本集。

一般而言均匀分布Uniform(0,1)的样本是相对容易生成的。 通过线性同余发生器可以生成伪随机数，我们用确定性算法生成[0,1]之间的伪随机数序列后，
这些序列的各种统计指标和均匀分布Uniform(0,1)的理论计算结果非常接近。这样的伪随机序列就有比较好的统计性质，可以被当成真实的随机数使用。线性同余随机数生成器如下:

$$x_{n+1}=(ax_n+c)~\textrm{mod}~m$$

式中$a，c，m$是数学推导出的合适的常数。这种算法产生的下一个随机数完全依赖当前的随机数，当随机数序列足够大的时候，随机数会出现重复子序列的情况。
当然，也有很多更加先进的随机数产生算法出现，比如numpy用的是 Mersenne Twister 等.

而其他常见的概率分布，无论是离散的分布还是连续的分布，它们的样本都可以通过Uniform(0,1)的样本转换而得, 但是如何产生满足其他分布下的随机数呢？

比如二维正态分布的样本$(Z1,Z2)$可以通过通过独立采样得到的uniform(0,1)样本对$(X1,X2)$通过如下的式子转换而得：

$$Z_1 = \sqrt{-2 ln X_1}cos(2\pi X_2) \\ Z_2 = \sqrt{-2 ln X_1}sin(2\pi X_2)$$

其他一些常见的连续分布，比如t分布，F分布，Beta分布，Gamma分布等，都可以通过类似的方式从uniform(0,1)得到的采样样本转化得到。在python的numpy，scikit-learn等类库中，都有生成这些常用分布样本的函数可以使用。

不过很多时候，我们的x的概率分布不是常见的分布，这意味着我们没法方便的得到这些非常见的概率分布的样本集。那这个问题怎么解决呢？

接受-拒绝采样

对于概率分布不是常见的分布，一个可行的办法是采用接受-拒绝采样来得到该分布的样本。既然 $p(x)$ 太复杂在程序中没法直接采样，那么我设定一个程序可采样的分布 $q(x$) 比如高斯分布，然后按照一定的方法拒绝某些样本，以达到接近 $p(x)$ 分布的目的，其中$q(x)$叫做 proposal distribution。

具体采用过程如下，设定一个方便采样的常用概率分布函数 q(x)，以及一个常量 k，使得 p(x) 总在 kq(x) 的下方。

首先，采样得到q(x)的一个样本$z_0$，采样方法如上，使用uniform(0,1)转换得到。然后，从均匀分布$(0,kq(z_0))$中采样得到一个值$u$。如果$u$落在了上图中的灰色区域，则拒绝这次抽样，否则接受这个样本$z_0$。重复以上过程得到 n 个接受的样本 $z_0,z_1,…z_{n−1}$,则最后的蒙特卡罗方法求解结果为：

$$\frac{1}{n}\sum\limits_{i=0}^{n-1}\frac{f(z_i)}{p(z_i)}$$

整个过程中，我们通过一系列的接受拒绝决策来达到用q(x)模拟p(x)概率分布的目的。

小结

使用接受-拒绝采样，我们可以解决一些概率分布不是常见的分布的时候，得到其采样集并用蒙特卡罗方法求和的目的。但是接受-拒绝采样也只能部分满足我们的需求，在很多时候我们还是很难得到我们的概率分布的样本集。比如：

1. 对于一些二维分布$p(x,y)$，有时候我们只能得到条件分布$p(x|y)$和$p(y|x)$,却很难得到二维分布$p(x,y)$一般形式，这时我们无法用接受-拒绝采样得到其样本集。
2. 对于一些高维的复杂非常见分布$p(x_1,x_2,…,x_n)$，我们要找到一个合适的$q(x)和$k$非常困难。

从上面可以看出，要想将蒙特卡罗方法作为一个通用的采样模拟求和的方法，必须解决如何方便得到各种复杂概率分布的对应的采样样本集的问题。

此时就需要使用一些更加复杂的随机模拟的方法来生成样本。比如马尔科夫链蒙特卡罗方法，了解这个算法我们首先要对马尔科夫链的平稳分布的性质有基本的认识。