马尔科夫链

Posted 上帝不玩骰子

tags:

篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了马尔科夫链相关的知识,希望对你有一定的参考价值。

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

1. 马尔科夫链概述

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

    如果用精确的数学定义来描述,则假设我们的序列状态是...Xt2,Xt1,Xt,Xt+1,......Xt−2,Xt−1,Xt,Xt+1,...,那么我们的在时刻Xt+1Xt+1的状态的条件概率仅仅依赖于时刻XtXt,即:

P(Xt+1|...Xt2,Xt1,Xt)=P(Xt+1|Xt)P(Xt+1|...Xt−2,Xt−1,Xt)=P(Xt+1|Xt)

 

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

    这个马尔科夫链是表示股市模型的,共有三种状态:牛市(Bull market), 熊市(Bear market)和横盘(Stagnant market)。每一个状态都以一定的概率转化到下一个状态。比如,牛市以0.025的概率转化到横盘的状态。这个状态概率转化图可以以矩阵的形式表示。如果我们定义矩阵阵PP某一位置P(i,j)P(i,j)的值为P(j|i)P(j|i),即从状态i转化到状态j的概率,并定义牛市为状态0, 熊市为状态1, 横盘为状态2. 这样我们得到了马尔科夫链模型的状态转移矩阵为:

P=⎛⎝⎜0.90.150.250.0750.80.250.0250.050.5⎞⎠⎟P=(0.90.0750.0250.150.80.050.250.250.5)

 

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

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

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

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

复制代码
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.625   0.3125  0.0625],即62.5%的牛市,31.25%的熊市与6.25%的横盘。那么这个是巧合吗?

    我们现在换一个初始概率分布试一试,现在我们用[0.7,0.1,0.2]作为初始概率分布,然后这个状态作为序列概率分布的初始状态t0t0,将其带入这个状态转移矩阵计算t1,t2,t3...t1,t2,t3...的状态。代码如下:

复制代码
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.625   0.3125  0.0625], 也就是说我们的马尔科夫链模型的状态转移矩阵收敛到的稳定概率分布与我们的初始状态概率分布无关。这是一个非常好的性质,也就是说,如果我们得到了这个稳定概率分布对应的马尔科夫链模型的状态转移矩阵,则我们可以用任意的概率分布样本开始,带入马尔科夫链模型的状态转移矩阵,这样经过一些序列的转换,最终就可以得到符合对应稳定概率分布的样本。

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

    同时,对于一个确定的状态转移矩阵PP,它的n次幂PnPn在当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]]  

    我们可以发现,在n6n≥6以后,PnPn的值稳定不再变化,而且每一行都为[0.625   0.3125  0.0625],这和我们前面的稳定分布是一致的。这个性质同样不光是离散状态,连续状态时也成立。

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

    如果一个非周期的马尔科夫链有状态转移矩阵PP, 并且它的任何两个状态是连通的,那么limnPnijlimn→∞Pijn与i无关,我们有:

    1)

limnPnij=π(j)limn→∞Pijn=π(j)

 

    2)

limnPn=⎛⎝⎜⎜⎜⎜⎜⎜π(1)π(1)π(1)π(2)π(2)π(2)π(j)π(j)π(j)⎞⎠⎟⎟⎟⎟⎟⎟limn→∞Pn=(π(1)π(2)…π(j)…π(1)π(2)…π(j)………………π(1)π(2)…π(j)………………)

 

    3)

π(j)=i=0π(i)Pijπ(j)=∑i=0∞π(i)Pij

 

    4) ππ是方程πP=ππP=π的唯一非负解,其中:

π=[π(1),π(2),...,π(j),...]i=0π(i)=1π=[π(1),π(2),...,π(j),...]∑i=0∞π(i)=1

 

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

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

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

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

    4)ππ通常称为马尔科夫链的平稳分布。

3. 基于马尔科夫链采样

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

    假设我们任意初始的概率分布是π0(x)π0(x), 经过第一轮马尔科夫链状态转移后的概率分布是π1(x)π1(x),。。。第i轮的概率分布是πi(x)πi(x)。假设经过n轮后马尔科夫链收敛到我们的平稳分布π(x)π(x),即:

πn(x)=πn+1(x)=πn+2(x)=...=π(x)πn(x)=πn+1(x)=πn+2(x)=...=π(x)

 

    对于每个分布πi(x)πi(x),我们有:

πi(x)=πi1(x)P=πi2(x)P2=π0(x)Piπi(x)=πi−1(x)P=πi−2(x)P2=π0(x)Pi

 

    现在我们可以开始采样了,首先,基于初始任意简单概率分布比如高斯分布π0(x)π0(x)采样得到状态值x0x0,基于条件概率分布P(以上是关于马尔科夫链的主要内容,如果未能解决你的问题,请参考以下文章

深度剖析HMM(附Python代码)1.前言及隐马尔科夫链HMM的背景

深度剖析HMM(附Python代码)2.隐马尔科夫链HMM的EM训练过程

马尔可夫链模型实际运用(以金融领域为例)(超详细代码)

蒙特卡罗方法生成指定状态空间下对应长度的马尔可夫链--MATLAB源程序

蒙特卡罗方法生成指定状态空间下对应长度的马尔可夫链--MATLAB源程序

机器学习算法之——隐马尔可夫模型(Hidden Markov Models,HMM) 代码实现