MCMC(二)蒙特卡洛之马尔可夫链

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

马尔科夫链简述

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

马尔可夫链的数学定义如下：
$P(X_{t+1} |...X_{t-2}, X_{t-1}, X_{t} ) = P(X_{t+1} | X_{t})$P(Xt+1​∣...Xt−2​,Xt−1​,Xt​)=P(Xt+1​∣Xt​)
也就是某一时刻状态转移的概率只依赖于它的前一个状态。下面来看一个例子(来自LDA数学八卦)。
社会学家经常把人按其经济状况分成3类:下层(lower-claaa}、中层(middle-claaa}、上层(upper-class，我们用1,2,3分别代表这三个阶层。社会学家们发现决定一个人的收入阶层的最重要的因素就是其父母的收入阶层。如果一个人的收入属于下层类别，那么他的孩子属于下层收入的概率是0.65,属于中层收入的概率是0.28,属于上层收入的概率是0.07。事实上，从父代到子代，收入阶层的变化的转移概率如下：

使用矩阵的表示方式，转移概率矩阵记为：
$P=\left( \begin{array}{ccc} 0.65&0.28&0.07\\ 0.15&0.67& 0.18\\ 0.12&0.36&0.52 \end{array} \right)$P=⎝⎛​0.650.150.12​0.280.670.36​0.070.180.52​⎠⎞​
假设当前这一代人处在下层、中层、上层的人的比例是概率分布向量$\pi_0=[\pi_0(1),\pi_0(2),\pi_0(3)]$π0​=[π0​(1),π0​(2),π0​(3)],那么他们的子女的分布比例将是:
$\pi_1=\pi_0P$π1​=π0​P,他们的孙子代的分布比例将是$\pi_2=\pi_1P=\pi_0P^2$π2​=π1​P=π0​P2,…第n代子孙的收入分布比例将是：$\pi_n=\pi_{n-1}P=\pi_0P^n$πn​=πn−1​P=π0​Pn.
假设初始概率分布为$\pi_0=[0.21,0.68,0.11]$π0​=[0.21,0.68,0.11],那么我们可以计算前n代人的收入分布情况如下：

我们发现从第七代开始，这个分布就稳定不变了，这个是偶然的么？我们换一个初始概率分布$\pi_0=[0.75,0.15,0.1]$π0​=[0.75,0.15,0.1]试试看，继续计算前n代人的分布情况：

我们发现从第九代人开始，分布又收敛了。并且，两次初始概率分布不同，最终收敛得到的分布却是一致的。在计算一下$P^n$Pn:
$P^{20}=P^{21}=.....=P^{100}=\left( \begin{array}{ccc} 0.286&0.489&0.225\\ 0.286&0.489&0.225\\ 0.286&0.489&0.225 \end{array} \right)$P20=P21=.....=P100=⎝⎛​0.2860.2860.286​0.4890.4890.489​0.2250.2250.225​⎠⎞​
上面的情况足以说明，最终收敛得到的稳定分布与初始概率分布无关，只与状态转移矩阵有关。

马尔可夫链的收敛性质

如果一个非周期的马尔科夫链有状态转移矩阵P, 并且它的任何两个状态是连通的，那么$\lim_{n \to \infty}P_{ij}^n与i$limn→∞​Pijn​与i无关，我们有：
1）
$\lim_{n \to \infty}P_{ij}^n = \pi(j)$n→∞lim​Pijn​=π(j)
2)
$\lim_{n \to \infty}P^n = \left( \begin{array}{ccc} \pi(1)&\pi(2)&\ldots&\pi(j)&\ldots \\ \pi(1)&\pi(2)&\ldots&\pi(j)&\ldots \\ \ldots&\ldots&\ldots&\ldots&\ldots \\ \pi(1)&\pi(2)&\ldots&\pi(j)&\ldots \\ \ldots&\ldots&\ldots&\ldots&\ldots \end{array} \right)$n→∞lim​Pn=⎝⎜⎜⎜⎜⎛​π(1)π(1)…π(1)…​π(2)π(2)…π(2)…​……………​π(j)π(j)…π(j)…​……………​⎠⎟⎟⎟⎟⎞​
3)
$\pi(j) = \sum\limits_{i=0}^{\infty}\pi(i)P_{ij}$π(j)=i=0∑∞​π(i)Pij​
4)$\pi$π是方程$\pi=\pi$π=πP的唯一非负解,其中：
$\pi = [\pi(1),\pi(2),...,\pi(j),...]\;\; \sum\limits_{i=0}^{\infty}\pi(i) = 1$π=[π(1),π(2),...,π(j),...]i=0∑∞​π(i)=1
$\pi$π通常称为马尔科夫链的平稳分布,这个马尔可夫链的收敛定理十分重要，所有的MCMC方法都是以这个定理为理论基础的。

假设我们用$X_i$Xi​表示在马氏链上跳转到第$i$i步后所处的状态，从初始概率分布$\pi_0$π0​出发，我们在马氏链上做状态转移，记$X_i$Xi​的概率分布为$\pi_i$πi​，则有：$X_0\sim\pi_0(x) ,X_i\sim\pi_i(x)$X0​∼π0​(x),Xi​∼πi​(x),其中，$\pi_i(x)=\pi_{i-1}(x)P=\pi_0(x)P^n$πi​(x)=πi−1​(x)P=π0​(x)Pn，由马氏链的收敛定理，概率分布$\pi_i(x)$πi​(x)将收敛到平稳分布$\pi(x)$π(x)。假设第n步的时候马氏链收敛，则有：
$\left( \begin{array}{ccc} X_0\sim\pi_0(x)\\ X_1\sim\pi_1(x)\\ \ldots\\ X_n\sim\pi_n(x)=\pi(x)\\ X_{n+1}\sim\pi_{n+1}(x)\\ X_{n+2}\sim\pi_{n+2}(x)\\ \ldots\end{array} \right)$⎝⎜⎜⎜⎜⎜⎜⎜⎜⎛​X0​∼π0​(x)X1​∼π1​(x)…Xn​∼πn​(x)=π(x)Xn+1​∼πn+1​(x)Xn+2​∼πn+2​(x)…​⎠⎟⎟⎟⎟⎟⎟⎟⎟⎞​
所以，$X_n,X_{n+1},X_{n+2},\ldots\sim\pi(x)$Xn​,Xn+1​,Xn+2​,…∼π(x)都是同分布的随机变量，当然他们并不独立。如果我们从一个具体的初始状态$x_0$x0​开始，沿着马氏链按照状态转移矩阵做跳转，那么我们得到的一个转移序列$x_0,x_1,x_2,\ldots x_n,x_{n+1} \ldots$x0​,x1​,x2​,…xn​,xn+1​…,由于马氏链的收敛行为，$x_n,x_{n+1} \ldots$xn​,xn+1​…都将是平稳分布$\pi(x)$π(x)的样本。