MCMC(一)：蒙特卡罗方法和马尔科夫链

作为一种随机采样方法，马尔科夫链蒙特卡罗（Markov Chain Monte Carlo，以下简称MCMC）在机器学习,深度学习以及自然语言处理等领域都有广泛的应用，是很多复杂算法求解的基础。比如分解机(Factorization Machines)推荐算法或者受限玻尔兹曼机（RBM）的原理，都用到了MCMC来做一些复杂运算的近似求解。下面我们就对MCMC的原理做一个总结。

从名字我们可以看出，MCMC由两个MC组成，即蒙特卡罗方法（Monte Carlo Simulation，简称MC）和马尔科夫链（Markov Chain ，也简称MC）。要弄懂MCMC的原理我们首先得搞清楚蒙特卡罗方法和马尔科夫链的原理。我们将用两篇来完整学习MCMC。在本篇，我们分别讲解蒙特卡罗方法和马尔科夫链。

本篇博客主要转自参考文献【1】【2】，在原文的基础上，为了更容易增加理解，略有删改增。

一、蒙特卡罗方法

1.1 蒙特卡罗方法引入

蒙特卡罗原来是一个赌场的名称，用它作为名字大概是因为蒙特卡罗方法是一种随机模拟的方法，这很像赌博场里面的扔骰子的过程。最早的蒙特卡罗方法都是为了求解一些不太好求解的求和或者积分问题。比如积分：
$$\theta ={\int }_a^{b}f(x)dx$$

如果我们很难求解出$f(x)$的原函数，那么这个积分比较难求解。当然我们可以通过蒙特卡罗方法来模拟求解近似值。如何模拟呢？假设我们函数图像如下图:

则一个简单的近似求解方法是在$[a,b]$之间随机的采样一个点。比如$x_0$,然后用$f(x_0)$代表在$[a,b]$区间上所有的$f(x)$的值。那么上面的定积分的近似求解为:
$$(b-a)f(x_0)$$

当然，用一个值代表$[a,b]$区间上所有的$f(x)$的值，这个假设太粗糙。那么我们可以采样$[a,b]$区间的$n$个值：$x_0,x_1,...x_{n-1}$,用它们的均值来代表$[a,b]$区间上所有的$f(x)$的值。这样我们上面的定积分的近似求解为:
$$\frac{b-a}{n}\sum _{i=0}^{n-1}f(x_i)$$

虽然上面的方法可以一定程度上求解出近似的解，但是它隐含了一个假定，即$x$在$[a,b]$之间是均匀分布的，而绝大部分情况，$x$在$[a,b]$之间不是均匀分布的。如果我们用上面的方法，则模拟求出的结果很可能和真实值相差甚远。

怎么解决这个问题呢？ 如果我们可以得到$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 _{i=0}^{n-1}\frac{f(x_i)}{p(x_i)}\text{(1)}$$

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

可以看出，最上面我们假设$x$在$[a,b]$之间是均匀分布的时候，$p(x_i)=1/(b-a)$，带入我们有概率分布的蒙特卡罗积分的(1)式，可以得到：
$$\frac{1}{n}\sum _{i=0}^{n-1}\frac{f(x_i)}{1/(b-a)}=\frac{b-a}{n}\sum _{i=0}^{n-1}f(x_i)$$

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

---

注：关于（1）式，我们可以这么理解最后一步转换：由于${\int }_a^b\frac{f(x)}{p(x)}p(x)dx$可以看做是$\frac{f(x)}{p(x)}$基于概率分布$p(x)$的期望，那么我们可以用期望的方法来求这个式子的值。而计算期望的一个近似方法是取$\frac{f(x)}{p(x)}$的若干个基于分布$p(x)$的采样点，然后求平均值得到。

对于（1）式，我们再举一个例子：$f(x)$的取值只有2个，$x=1，2$。对应的$y$值分别是$f(1)=1,f(2)=4$。其中$x$的取值不是平均的，取1的概率$p(1)=0.25$,取2的概率是$p(2)=0.75$。

那么严格来说，根据（1）式，对应的$f(x)$的积分等于$\frac{f(x)}{p(x)}$基于概率分布$p(x)$的期望，根据期望的公式，则有$\frac{1}{0.25}\ast 0.25+\frac{4}{0.75}\ast 0.75=5$。

此时我们去采样三次，期望求近似结果。发现第一次采样到1，第二次和第三次采样到2，那么最后的近似结果是$\frac{1}{3}(\frac{1}{0.25}+(\frac{4}{0.75}+(\frac{4}{0.75})=4.89$。

这个4.89就是我们5的近似。虽然有些距离，但是是是由于采样太少的原因。

假设我们采样100次，得到26次1，74次2，那么最后的近似结果是$\frac{1}{100}(\frac{1}{0.25}\times 26+(\frac{4}{0.75}\times 74)=4.99$

可见越来越接近的。接近的原因是随着采样数的增多，采样的样本的分布越来越接近于$x$本来的分布。

---

1.2 概率分布采样和拒绝采样

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

对于常见的分布，如均匀分布，高斯分布，指数分布，t分布，F分布，Beta分布，Gamma分布等，可以采用逆采样的方法进行采样；关于逆采样的知识请参看这篇博客《逆采样(Inverse Sampling)和拒绝采样(Reject Sampling)原理详解》

不过很多时候，我们的$x$的概率分布不是常见的分布，这些分布的概率分布函数CDF不可逆，因此没有办法用逆采样来采样，这意味着我们没法方便的得到这些非常见的概率分布的样本集。拒绝采样就是用来解决这个问题。关于拒绝采样的知识也请参看这篇博客《逆采样(Inverse Sampling)和拒绝采样(Reject Sampling)原理详解》

这样。我们得到采样的$n$个样本$z_0,z_1,...z_{n-1}$后，最后的蒙特卡罗方法求解结果为：
$$\frac{1}{n}\sum _{i=0}^{n-1}\frac{f(z_i)}{p(z_i)}$$

1.3 蒙特卡罗方法小结

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

1）对于一些二维分布$p(x,y)$，有时候我们只能得到条件分布$p(x|y)$和$p(y|x)$，却很难得到二维分布的概率密度函数$p(x,y)$的一般形式，这时我们无法用拒绝采样得到其样本集。

2）对于一些高维的复杂非常见分布$p(x_1,x_2,...,x_n)$，我们要找到一个合适的$q(x)$和$k$非常困难。

因此，实际上，我们仍然要找到一种方法可以解决如何方便得到各种复杂概率分布的对应的采样样本集的问题。马尔科夫链有能力帮助找到这些复杂概率分布的对应的采样样本集。因此，下一章我们会讲解马尔科夫链。

二、马尔科夫链

在上一章中，我们讲到了如何用蒙特卡罗方法来随机模拟求解一些复杂的连续积分或者离散求和的方法，但是这个方法需要得到对应的概率分布的样本集，而想得到这样的样本集很困难。因此我们需要本篇讲到的马尔科夫链来帮忙。

2.1 马尔科夫链概述

马尔科夫链定义本身比较简单，它假设某一时刻状态转移的概率只依赖于它的前一个状态。举个形象的比喻，假如每天的天气是一个状态的话，那个今天是不是晴天只依赖于昨天的天气，而和前天的天气没有任何关系。当然这么说可能有些武断，但是这样做可以大大简化模型的复杂度，因此马尔科夫链在很多时间序列模型中得到广泛的应用，比如循环神经网络RNN，隐式马尔科夫模型HMM等，当然MCMC也需要它。

如果用精确的数学定义来描述，则假设我们的序列状态是$...X_{t-2},X_{t-1},X_t,X_{t+1},...$，那么我们的在时刻$X_{t+1}$的状态的条件概率仅仅依赖于时刻$X_t$，即：
$$P(X_{t+1}|...X_{t-2},X_{t-1},X_t)=P(X_{t+1}|X_t)$$

既然某一时刻状态转移的概率只依赖于它的前一个状态，那么我们只要能求出系统中任意两个状态之间的转换概率，这个马尔科夫链的模型就定了。我们来看看下图这个马尔科夫链模型的具体的例子(来源于维基百科)。

这个马尔科夫链是表示股市模型的，共有三种状态：牛市（Bull market）, 熊市（Bear market）和横盘（Stagnant market）。每一个状态都以一定的概率转化到下一个状态。比如，牛市以0.025的概率转化到横盘的状态。这个状态概率转化图可以以矩阵的形式表示。如果我们定义矩阵$P$某一位置$P(i,j)$的值为$P(j|i)$,即从状态$i$转化到状态$j$的概率，并定义牛市为状态0， 熊市为状态1, 横盘为状态2。这样我们得到了马尔科夫链模型的状态转移矩阵为：
$$P=\left(\begin{array}{ccc}0.9& 0.075& 0.025\\ 0.15& 0.8& 0.05\\ 0.25& 0.25& 0.5\end{array}\right)$$

讲了这么多，那么马尔科夫链模型的状态转移矩阵和我们蒙特卡罗方法需要的概率分布样本集有什么关系呢？这需要从马尔科夫链模型的状态转移矩阵的性质讲起。

2.2 马尔科夫链模型状态转移矩阵P的性质

得到了马尔科夫链模型的状态转移矩阵，我们来看看马尔科夫链模型的状态转移矩阵的性质。

2.2.1 观察1

仍然以上面的这个状态转移矩阵为例。假设我们当前股市的概率分布为：$[0.3,0.4,0.3]$，即30%概率的牛市，40%概率的熊盘与30%的横盘。然后这个状态作为序列概率分布的初始状态$t_0$，将其带入这个状态转移矩阵计算$t_1,t_2,t_3...$的状态。代码如下：

import numpy as np
matrix = np.matrix([[0.9,0.075,0.025],[0.15,0.8,0.05],[0.25,0.25,0.5]], dtype=float)
vector1 = np.matrix([[0.3,0.4,0.3]], dtype=float)
for i in range(100):
    vector1 = vector1*matrix
    print "Current round:" , i+1
    print vector1

部分输出结果如下：

Current round: 1
[[ 0.405   0.4175  0.1775]]
Current round: 2
[[ 0.4715   0.40875  0.11975]]
Current round: 3
[[ 0.5156  0.3923  0.0921]]
Current round: 4
[[ 0.54591   0.375535  0.078555]]
。。。。。。
Current round: 58
[[ 0.62499999  0.31250001  0.0625    ]]
Current round: 59
[[ 0.62499999  0.3125      0.0625    ]]
Current round: 60
[[ 0.625   0.3125  0.0625]]
。。。。。。
Current round: 99
[[ 0.625   0.3125  0.0625]]
Current round: 100
[[ 0.625   0.3125  0.0625]]

可以发现，从第60轮开始，我们的状态概率分布就不变了，一直保持在$[0.6250.31250.0625]$，即62.5%的牛市，31.25%的熊市与6.25%的横盘。那么这个是巧合吗？

我们现在换一个初始概率分布试一试，现在我们用$[0.7,0.1,0.2]$作为初始概率分布，然后这个状态作为序列概率分布的初始状态$t_0$，将其带入这个状态转移矩阵计算$t_1,t_2,t_3...$的状态。代码如下：

matrix = np.matrix([[0.9,0.075,0.025],[0.15,0.8,0.05],[0.25,0.25,0.5]], dtype=float)
vector1 = np.matrix([[0.7,0.1,0.2]], dtype=float)
for i in range(100):
    vector1 = vector1*matrix
    print "Current round:" , i+1
    print vector1

部分输出结果如下：

Current round: 1
[[ 0.695   0.1825  0.1225]]
Current round: 2
[[ 0.6835   0.22875  0.08775]]
Current round: 3
[[ 0.6714  0.2562  0.0724]]
Current round: 4
[[ 0.66079   0.273415  0.065795]]
。。。。。。。
Current round: 55
[[ 0.62500001  0.31249999  0.0625    ]]
Current round: 56
[[ 0.62500001  0.31249999  0.0625    ]]
Current round: 57
[[ 0.625   0.3125  0.0625]]
。。。。。。。
Current round: 99
[[ 0.625   0.3125  0.0625]]
Current round: 100
[[ 0.625   0.3125  0.0625]]

可以看出，尽管这次我们采用了不同初始概率分布，最终状态的概率分布趋于同一个稳定的概率分布$[0.6250.31250.0625]$， 也就是说我们的马尔科夫链模型的状态转移矩阵收敛到的稳定概率分布与我们的初始状态概率分布无关。这是一个非常好的性质，也就是说，如果我们得到了这个稳定概率分布对应的马尔科夫链模型的状态转移矩阵，则我们可以用任意的概率分布样本开始，带入马尔科夫链模型的状态转移矩阵，这样经过一些序列的转换，最终就可以得到符合对应稳定概率分布的样本。

这个性质不光对我们上面的状态转移矩阵有效，对于绝大多数的其他的马尔科夫链模型的状态转移矩阵也有效。同时不光是离散状态，连续状态时也成立。

2.2.2 观察2

同时，对于一个确定的状态转移矩阵$P$，它的$n$次幂$P^n$在当$n$大于一定的值的时候也可以发现是确定的，我们还是以上面的例子为例，计算代码如下：

matrix = np.matrix([[0.9,0.075,0.025],[0.15,0.8,0.05],[0.25,0.25,0.5]], dtype=float)
for i in range(10):
    matrix = matrix*matrix
    print "Current round:" , i+1
    print matrix

输出结果如下：

Current round: 1
[[ 0.8275   0.13375  0.03875]
 [ 0.2675   0.66375  0.06875]
 [ 0.3875   0.34375  0.26875]]
Current round: 2
[[ 0.73555   0.212775  0.051675]
 [ 0.42555   0.499975  0.074475]
 [ 0.51675   0.372375  0.110875]]
。。。。。。
Current round: 5
[[ 0.62502532  0.31247685  0.06249783]
 [ 0.6249537   0.31254233  0.06250397]
 [ 0.62497828  0.31251986  0.06250186]]
Current round: 6
[[ 0.625   0.3125  0.0625]
 [ 0.625   0.3125  0.0625]
 [ 0.625   0.3125  0.0625]]
Current round: 7
[[ 0.625   0.3125  0.0625]
 [ 0.625   0.3125  0.0625]
 [ 0.625   0.3125  0.0625]]
。。。。。。
Current round: 9
[[ 0.625   0.3125  0.0625]
 [ 0.625   0.3125  0.0625]
 [ 0.625   0.3125  0.0625]]
Current round: 10
[[ 0.625   0.3125  0.0625]
 [ 0.625   0.3125  0.0625]
 [ 0.625   0.3125  0.0625]]　

我们可以发现，在$n\ge 6$以后，$P^n$的值稳定不再变化，而且每一行都为$[0.6250.31250.0625]$，这和我们前面的稳定分布是一致的。这个性质同样不光是离散状态，连续状态时也成立。

2.2.3 性质总结

好了，现在我们可以用数学语言总结下马尔科夫链的收敛性质了：

如果一个非周期的马尔科夫链有状态转移矩阵$P$, 并且它的任何两个状态是连通的，那么${\lim }_{n\to \infty }{P}_{ij}^n$与$i$无关，我们有：
1）
$$\underset{n\to \infty }{\lim }{P}_{ij}^n=\pi (j)$$
$n\to \infty$即无数个$P$相乘的含义，由2.2.2观察2，可以得到性质1。
2）
$$\underset{n\to \infty }{\lim }{P}^n=\left(\begin{array}{ccccc}\pi (1)& \pi (2)& \dots & \pi (j)& \dots \\ \pi (1)& \pi (2)& \dots & \pi (j)& \dots \\ \dots & \dots & \dots & \dots & \dots \\ \pi (1)& \pi (2)& \dots & \pi (j)& \dots \\ \dots & \dots & \dots & \dots & \dots \end{array}\right)$$
由2.2.2观察2，可以得到性质2。
3）
$$\pi (j)=\sum _{i=0}^{\infty }\pi (i)P_{ij}$$

求和 $i=0$至无穷，就是不停的计算$\pi (i)P_{ij}$得到的结果再去乘$P_{ij}$，直到收敛为止。这样第$i$个分量会收敛到平稳分布对应序号的分量。由2.2.1观察1，可以得到性质3。

4） $\pi$是方程$\pi P=\pi$的唯一非负解，其中：
$$\pi =[\pi (1),\pi (2),...,\pi (j),...]\sum _{i=0}^{\infty }\pi (i)=1$$

由2.2.1观察1，可以得到性质4。

上面的性质中需要解释的有：

1）非周期的马尔科夫链：这个主要是指马尔科夫链的状态转化不是循环的，如果是循环的则永远不会收敛。幸运的是我们遇到的马尔科夫链一般都是非周期性的。用数学方式表述则是：对于任意某一状态$i$，$d$为集合 { n ∣ n ≥ 1 , P i i n > 0 } \{n \mid n \geq 1,P_{ii}^n>0 \} {n∣n≥1,Piin​>0}的最大公约数，如果$d=1$，则该状态为非周期的。

2）任何两个状态是连通的：这个指的是从任意一个状态可以通过有限步到达其他的任意一个状态，不会出现条件概率一直为0导致不可达的情况。

3）马尔科夫链的状态数可以是有限的，也可以是无限的。因此可以用于连续概率分布和离散概率分布。

4）$\pi$通常称为马尔科夫链的平稳分布。

2.3 基于马尔科夫链采样

如果我们得到了某个平稳分布所对应的马尔科夫链状态转移矩阵$P$，我们就很容易采用出这个平稳分布的样本集。

假设我们任意初始的概率分布是${\pi }_0(x)$, 经过第一轮马尔科夫链状态转移后的概率分布是${\pi }_1(x)$，。。。第$i$轮的概率分布是${\pi }_i(x)$。假设经过 $n$ 轮后马尔科夫链收敛到我们的平稳分布$\pi (x)$，即：
$${\pi }_n(x)={\pi }_{n+1}(x)={\pi }_{n+2}(x)=...=\pi (x)$$

对于每个分布${\pi }_i(x)$，我们有：
$${\pi }_i(x)={\pi }_{i-1}(x)P={\pi }_{i-2}(x)P^2={\pi }_0(x)P^i$$

现在我们可以开始采样了，首先，基于初始任意简单概率分布比如高斯分布${\pi }_0(x)$采样得到状态值$x_0$，基于条件概率分布$P(x|x_0)$采样状态值$x_1$，一直进行下去，当状态转移进行到一定的次数时，比如到$n$次时，我们认为此时的采样集$(x_n,x_{n+1},x_{n+2},...)$即是符合我们的平稳分布的对应样本集，可以用来做蒙特卡罗模拟求和了。

总结下基于马尔科夫链的采样过程：

1）输入马尔科夫链状态转移矩阵$P$，设定状态转移次数阈值$n_1$，需要的样本个数$n_2$

2）从任意简单概率分布采样得到初始状态值$x_0$

3）for $t=0$ to $n_1+n_2-1$: 从条件概率分布$P(x|x_t)$中采样得到样本$x_{t+1}$。

4）样本集$(x_{n_1},x_{n_1+1},...,x_{n_1+n_2-1})$即为我们需要的平稳分布对应的样本集。

对采样过程的几点再说明一下：

• 对第3）步，每对$P(x|x_t)$采样一次就意味着模拟了一次${\pi }_i(x)={\pi }_{i-1}(x)P$操作。在此过程中，状态转化矩阵$P$没有变化。

• 对第4）步，注意的是采样的样本集的下标是$n_1,n_1+1,...,n_1+n_2+1$，即状态稳定之后采样的样本才是我们需要的样本集。即我们需要从真正趋于稳定时候的分布抽样样本。

2.4 马尔科夫链采样小结

如果假定我们可以得到我们需要采样样本的平稳分布所对应的马尔科夫链状态转移矩阵$P$，那么我们就可以用马尔科夫链采样得到我们需要的样本集，进而进行蒙特卡罗模拟。但是一个重要的问题是，随意给定一个平稳分布$\pi$,如何得到它所对应的马尔科夫链状态转移矩阵$P$呢？这是个大问题。我们绕了一圈似乎还是没有解决任意概率分布采样样本集的问题。

幸运的是，MCMC采样通过迂回的方式解决了上面这个大问题，我们接下来的文章来讨论MCMC的采样，以及它的使用改进版采样: M-H采样和Gibbs采样。链接如下：
MCMC(二)：MCMC采样和M-H采样
MCMC(三)：Gibbs采样

参考文献

【1】MCMC(一)蒙特卡罗方法

【2】MCMC(二)马尔科夫链

【2】 马尔可夫链蒙特卡罗算法（MCMC）