一.基本思路:马尔可夫链收敛到平稳态

《12_05_PGM_马尔科夫链_初探及代码实现》这一节,我们首次探索了马尔可夫链,并在本末讨论并说明了它的一个重要性质:马尔可夫链的平稳态,而且在后面的PageRank算法中使用了它的这个性质,这部分内容其实是要求这样一个问题:

**已知马尔可夫模型的状态转移概率矩阵$P$,求它的平稳态$\pi^*$,使得$P\pi^*=\pi^*$**

将这个问题反过来,那就是MCMC(马尔科夫蒙特卡洛抽样法)的主要思想咯:

**已知某分布$\pi$,求一个马尔可夫模型的状态转移概率矩阵$P^*$,使得$P^*\pi=\pi$**

于是,我们求一个$P^*$满足上面的条件即可,但是要用于MCMC,对$P^*$的要求更加苛刻一些:要保证在$P^*$的条件下,平稳态必须是唯一的,换言之,给定一个概率转移矩阵,它的平稳态可能有多个,比如下面的例子

In [1]:
import numpy as np
import pandas as pd
import os
os.chdir('../')
from ml_models.pgm import SimpleMarkovModel
In [2]:
smm=SimpleMarkovModel(status_num=3)
smm.P=np.asarray([
    [1,1/3,0],
    [0,1/3,0],
    [0,1/3,1]
])
In [3]:
#定义初始状态
init_prob=np.asarray([[0.2],[0.7],[0.1]])
pd.DataFrame({50:smm.predict_prob_distribution(1,set_init_prob=init_prob).reshape(-1).tolist(),
              100:smm.predict_prob_distribution(10,set_init_prob=init_prob).reshape(-1).tolist(),
              500:smm.predict_prob_distribution(50,set_init_prob=init_prob).reshape(-1).tolist(),
              1000:smm.predict_prob_distribution(100,set_init_prob=init_prob).reshape(-1).tolist(),
              2000:smm.predict_prob_distribution(200,set_init_prob=init_prob).reshape(-1).tolist()})
Out[3]:
50 100 500 1000 2000
0 0.433333 0.549994 5.500000e-01 5.500000e-01 5.500000e-01
1 0.233333 0.000012 9.750689e-25 1.358228e-48 2.635403e-96
2 0.333333 0.449994 4.500000e-01 4.500000e-01 4.500000e-01
In [4]:
#定义另一组初始状态
init_prob=np.asarray([[0.8],[0.1],[0.1]])
pd.DataFrame({50:smm.predict_prob_distribution(1,set_init_prob=init_prob).reshape(-1).tolist(),
              100:smm.predict_prob_distribution(10,set_init_prob=init_prob).reshape(-1).tolist(),
              500:smm.predict_prob_distribution(50,set_init_prob=init_prob).reshape(-1).tolist(),
              1000:smm.predict_prob_distribution(100,set_init_prob=init_prob).reshape(-1).tolist(),
              2000:smm.predict_prob_distribution(200,set_init_prob=init_prob).reshape(-1).tolist()})
Out[4]:
50 100 500 1000 2000
0 0.833333 0.849999 8.500000e-01 8.500000e-01 8.500000e-01
1 0.033333 0.000002 1.392956e-25 1.940325e-49 3.764862e-97
2 0.133333 0.149999 1.500000e-01 1.500000e-01 1.500000e-01
In [5]:
#再定义另一组初始状态
init_prob=np.asarray([[0.4],[0.5],[0.1]])
pd.DataFrame({50:smm.predict_prob_distribution(1,set_init_prob=init_prob).reshape(-1).tolist(),
              100:smm.predict_prob_distribution(10,set_init_prob=init_prob).reshape(-1).tolist(),
              500:smm.predict_prob_distribution(50,set_init_prob=init_prob).reshape(-1).tolist(),
              1000:smm.predict_prob_distribution(100,set_init_prob=init_prob).reshape(-1).tolist(),
              2000:smm.predict_prob_distribution(200,set_init_prob=init_prob).reshape(-1).tolist()})
Out[5]:
50 100 500 1000 2000
0 0.566667 0.649996 6.500000e-01 6.500000e-01 6.500000e-01
1 0.166667 0.000008 6.964778e-25 9.701626e-49 1.882431e-96
2 0.266667 0.349996 3.500000e-01 3.500000e-01 3.500000e-01

由上面的例子可以发现,由于初始状态的不同,平稳态可能有多个,所以我们构造的$P$必须要保证平稳态能唯一的收敛到我们的目标分布$\pi$,不然,如上面的例子,假如我们的目标分布是$\pi=[0.65,0,0.35]$,然后我们构造了一个概率转移矩阵$P=\begin{bmatrix} 1 & 1/3 & 0\\ 0 & 1/3 &0 \\ 0 & 1/3 & 1 \end{bmatrix}$ ,满足$P\pi=\pi$,然后,我们随便定义了一个初始点$x_0=[0.8,0.1,0.1]$,在$P$上随机游走采样$\{x_1,x_2,...,x_n\}$,最后发现它成功收敛到了$x_n\rightarrow [0.85,0,0.15]$,哈哈哈哈哈~~~~

保证马尔科夫模型仅收敛到唯一平稳态的充分条件是有的!那就是细致平衡方程

二.保证平稳态唯一:细致平衡方程

细致平衡方程的定义如下,设有马尔科夫链$X=\{X_0,X_1,...,X_t,...\}$,状态空间为$S$,转移概率矩阵为$P$,如果有状态分布$\pi=(\pi_1,\pi_2,...)^T$,对任意时刻状态$i,j\in S$,对任意时刻$t$均满足:

$$ P(X_t=i\mid X_{t-1}=j)\pi_j=P(X_{t-1}=j\mid X_t=i)\pi_i,i,j=1,2,.... $$

或者简写为:

$$ p_{ij}\pi_j=p_{ji}\pi_i,i,j=1,2,.... $$

为了更加直观,也可以写作:

$$ p(j\rightarrow i)\pi_j=p(i\rightarrow j)\pi_i $$

这便是细致平衡方程,此时的马尔科夫链称为可逆马尔科夫链,显然对$P$矩阵按列求和为1,即$\sum_{i}p_{ij}=1$或者$\sum_ip(j\rightarrow i)=1$ ,接下来简单证明一下,满足细致平衡条件,则有$P\pi=\pi$成立:

$$ (P\pi)_i=\sum_jp_{ij}\pi_j=\sum_jp_{ji}\pi_i=\pi_i\sum_{j}p_{ji}=\pi_i,i=1,2,... $$

高维/连续状态空间:上面只是对一维离散状态空间做的说明,而细致平衡方程对高维空间或者连续状态空间一样是成立的,只需对相应符号做修改即可,比如将$p_{ij}$修改为一个函数的形式$p(状态j,状态i)$,求和符号$\sum$需要替换为积分符号$\int$

三.小结

最后小结一下,这一节的重点引出细致平衡方程:

$$ p_{ij}\pi_j=p_{ji}\pi_i $$

因为在满足该方程约束的条件下,我们构造的$P^*$,必然有$P^*\pi=\pi$,且$\pi$唯一(这里$\pi$是我们的目标分布),这也是检验MCMC是否有效的一条重要标准,下一节将要介绍的Metropolis-Hastings算法(MH算法)便是满足该约束条件的一套算法

采样步骤

再对MCMC的采样步骤做个整理

(1)首先在随机变量$x$的状态空间上构造一个马尔科夫链$P$,使它能平稳且唯一的收敛到我们的目标分布$p(x)$;

(2)从状态空间某一点$x_0$出发,用构造的马尔科夫链$P$进行随机游走,产生样本序列$x_0,x_1,...,x_t,....$;

(3)确定一个正整数$m,n$($m<n$),得到样本集合$\{x_m,x_{m+1},...,x_n\}$,对目标指标进行计算,比如某函数$f(x)$的期望:

$$ \hat{E}f=\frac{1}{n-m}\sum_{i=m+1}^nf(x_i) $$
补充

(1)上面提到的随机游走采样简单解释一下,它其实和《12_06_PGM_马尔科夫链_语言模型及文本生成》中的文本生成类似,只是每一步不会做贪婪搜索或者beam search,而是按照当前步的概率分布进行采样,例如当前样本在状态$i$,则按照$p_{ji},j=1,2,...$的分布对下一个样本进行采样;

(2)另外由于马尔科夫链的性质,当$t\rightarrow \infty$时,采样的样本分布才会更加近似于目标分布,所以通常我们会从所有的采样样本$\{x_0,x_1,x_2,...,x_m,...,x_n\}$中截取尾部的一部分样本$\{x_m,x_{m+1},...,x_n\}$来做相应的统计指标计算,比如计算期望,积分等...

(3)另外另外,有时为了方便直观理解,转移概率我有时会采用不同的写法,但请注意一下它们是等价的:$p_{ij}=p(i\mid j)=p(j\rightarrow i)=p(j,i)$

In [ ]: