马尔可夫链例子

“随机过程随机过，实变函数学十遍，微机原理闹危机，汇编语言不会编”

1. 唯一让我彻底蒙圈的课程

这些课程真的太难了，大学里无数人为此伤透了脑筋，挂科率杠杠的。我当初也是的，特别是随机过程这门课，上完了一学期的课，只记住了几个公式，问我干嘛的？不知道！

像其他的高等数学啊，电磁场电磁波啊，通信原理啊，我都能大体知道是干嘛的，用在什么地方。讲真的，唯独就随机过程，感觉这门课太变态了，学的我云里雾里的，尤其是我当时那本只有公式，别无其他的影印教材，看了让人直蒙圈。

此后的很长一段时间，随机过程都是我的噩梦，是一段不忍回忆的历史，不过庆幸还好工作中不会接触到它，学不会就罢了。可是---直到我最近买了一本机器学习的书，翻开目录，看到了马尔可夫随机场时。。。哎，还是乖乖的去把马尔可夫链这个随机过程里很重要的概念搞搞清楚吧。

这几天又把之前的教材翻开复习了一遍，这才对随机过程和马尔可夫链有了更为深入的了解。

接下来，我用一个简单的例子，来给大家讲讲马尔可夫链到底是什么东东。

2. 随机过程是啥玩意儿

讲马尔可夫链不得不提到随机过程，因为它是随机过程课本中的，啊啊。

顾名思义，它其实就是个过程，比如今天下雨，那么明天下不下雨呢？后天下不下雨呢？从今天下雨到明天不下雨再到后天下雨，这就是个过程。那么怎么预测N天后到底下不下雨呢？这其实是可以利用公式进行计算的，随机过程就是这样一个工具，把整个过程进行量化处理，用公式就可以推导出来N天后的天气状况，下雨的概率是多少，不下雨的概率是多少。说白了，随机过程就是一些统计模型，利用这些统计模型可以对自然界的一些事物进行预测和处理，比如天气预报，比如股票，比如市场分析，比如人工智能。它的应用还真是多了去了。

话说回来，还真是佩服能构造出这些统计模型的大牛，简直脑洞大开啊。

好了，终于可以来看看马尔可夫链 （Markov Chain）到底是什么了。

它是随机过程中的一种过程，到底是哪一种过程呢？好像一两句话也说不清楚，还是先看个例子吧。

先说说我们村智商为0的王二狗，人傻不拉几的，见人就傻笑，每天中午12点的标配，仨状态：吃，玩，睡。这就是传说中的状态分布。

你想知道他n天后中午12点的状态么？是在吃，还是在玩，还是在睡？这些状态发生的概率分别都是多少？ （知道你不想，就假装想知道吧~~学习真的好累~~）

先看个假设，他每个状态的转移都是有概率的，比如今天玩，明天睡的概率是几，今天玩，明天也玩的概率是几几，看图更清楚一点。

这个矩阵就是转移概率矩阵P，并且它是保持不变的，就是说第一天到第二天的转移概率矩阵跟第二天到第三天的转移概率矩阵是一样的。（这个叫时齐，不细说了，有兴趣的同学自行百度）。

有了这个矩阵，再加上已知的第一天的状态分布，就可以计算出第N天的状态分布了。

S1 是4月1号中午12点的的状态分布矩阵 [0.6, 0.2, 0.2]，里面的数字分别代表吃的概率，玩的概率，睡的概率。

那么

4月2号的状态分布矩阵 S2 = S1 * P (俩矩阵相乘)。

4月3号的状态分布矩阵 S3 = S2 * P (看见没，跟S1无关，只跟S2有关)。

4月4号的状态分布矩阵 S4 = S3 * P (看见没，跟S1，S2无关，只跟S3有关)。

...

4月n号的状态分布矩阵 Sn = Sn-1 * P (看见没，只跟它前面一个状态Sn-1有关)。

总结：马尔可夫链就是这样一个任性的过程，它将来的状态分布只取决于现在，跟过去无关！

就把下面这幅图想象成是一个马尔可夫链吧。实际上就是一个随机变量随时间按照Markov性进行变化的过程。

=================================更新============

附：S2 的计算过程 (没兴趣的同学自行略过)

欢迎关注我的微信公众号：红猴子

这是一个工科生涨知识的号，公众号的内容有CS\EE技术, 职场以及经验谈，知乎专栏文章会首发于我的微信公众号，希望能给迷茫和困惑中的朋友一些启发与帮助，欢迎围观

什么是马尔可夫过程(Markov Process)

要说什么是马尔可夫过程，首先必须讲讲什么是随机过程(Stochastic Process)。

設 (Ω,,P) 為一概率空間，另設集合 T為一指標集合。如果對於所有t∈T，均有一隨機變量 ξt(ω)定義於概率空間 (Ω,,P) ，則集合 {ξt(ω)|t∈T}為一隨機過程。—Wikipedia

其实简单来说，随机过程就是对一连串随机变量(或事件)变迁或者说动态关系的描述。

马尔可夫过程其实就是随机过程的一种：

A stochastic process st:t∈N is said to be Markov if

prob(st+1|st)=prob(st+1|st)

where st is the history up to period t and st is the realization of the state at period t.

换句话说，马尔可夫过程说的就是在随机过程中，明天的状态只取决于今天的状态而和今天之前的状态无关。也就是不管你是怎样达到 st 的，只要你今天的状态是 st, 你明天状态的可能性就决定好了。

我们经常还会看到马尔可夫链。其实马尔可夫链就是状态空间为可数集的马尔可夫过程：

Let st∈S and if S is finite (countable), we call the Markov Process with this countable state space Markov Chain.

转移矩阵(Transition Matrix)

既然当知道 st 时我们就能知道 st+1 发生的概率，那么我们自然可以把有限可数的状态空间的转移法(transistion law)写成一个矩阵。举个例子，如果状态空间只包含了两个状态A和B。如果今天是A,那么明天还是A的概率为0.8。如果今日是B,那么明天还是B的概率为0.7。转移矩阵可写为：

P=∣∣∣0.8 0.30.20.7∣∣∣

稳态(steady－state)

首先我们先定义一个横向量 π；

π＝∣∣prob(A)prob(B)∣∣

假设我们今天的状态是A而不是B,那么我们明天是A还是B的可能性为：

∣∣10∣∣×∣∣∣0.8 0.30.20.7∣∣∣=∣∣0.80.2∣∣

后天是A还是B的可能性为：

∣∣0.80.2∣∣×∣∣∣0.8 0.30.20.7∣∣∣=∣∣0.70.3∣∣

大后天是A还是B的可能性为：

∣∣0.70.3∣∣×∣∣∣0.8 0.30.20.7∣∣∣=∣∣0.650.35∣∣

那么会不会某一时刻后我们的横向量 π能到达一个不动点：

π∗×P=π∗

这么繁杂的工作还是要我们的小奴Matlab来完成吧。

%initial pi%
pi=[1,0];

%transistion matrix%
P=[0.8,0.2;0.3,0.7];

%critical value%
cv=1;

while cv>0.0000000001;
    pi2=pi*P;
    cv=abs(min(pi2-pi));
    pi=pi2;
end

得出的结果为:

π∗=∣∣0.60.4∣∣

让我们换一个initial pi来试试看:

%initial pi%
pi=[0.35,0.65];

%transistion matrix%
P=[0.8,0.2;0.3,0.7];

%critical value%
cv=1;

while cv>10e-06;
    pi2=pi*P;
    cv=abs(min(pi2-pi));
    pi=pi2;
end

得出的结果仍然为:

π∗=∣∣0.60.4∣∣

其实仔细观察一下就可以看出这个 π 其实是转移矩阵 P 的转置的对应特征值(Eigenvalue)等于1的特征向量(Eigenvector)：

   π=πPπ(I−P)=0(I−P′)π′=0(P′−I)π′=0

既然知道 π 其实是 P‘ 的特征向量，我们就可以直接通过Matlab找特征向量而不必要用繁琐的迭代法啦。

%%%%%Eigenvector%%%%%
[n,v]=eig(P');
ppi=n(:,1);
scale=sum(ppi); %ensure the summation of the distribution vector pi is 1
ppi=ppi/scale;

得出的结果绝对是当然地为:

π∗=∣∣0.60.4∣∣

一个简单的模型

假设我们在研究一个粒子的运动。这个粒子随机地在A门与B门之间跳动。如果这个粒子跳入A门，那么它下一次跳入A门的概率为0.8。如果这个粒子跳入B门，那么它下次跳入B门的概率为0.7。我们的问题是，现在有100000个这样的粒子让它们跳跃，让它们跳跃1000次,我们会观察到多少个粒子在A门多少个在B门？

思路：

首先随机生成100000个粒子。1代表该粒子在A门，0代表在B门。

 size=100000;
 a=round(rand(1,size));

生成10000个以概率0.8为1，概率0.2为0的随机数。同时在生成10000个以概率0.7为0，概率0.3为1随机数。

 AA=rand(1,size)<=0.8;
 BB=rand(1,size)<=0.3;

用2a－AA。如果值为1说明该粒子从A门跳转到A门；如果值为2则说明该粒子从A门跳转到B门。另外值为0，－1则说明该粒子初始状态是在B门。

用2a－BB。如果值为0说明该粒子从B门跳转到B门；如果值为-1则说明该粒子从B门跳转到A门。另外值为2，1则说明该粒子初始状态是在A门。

A=2*a-AA;
B=2*a-BB;

找到A中大于等于1的元素与B中小于等于0的元素。同时新建一个向量a2用于存储结果。

indxA=find(A>=1);
indxB=find(B<=0);

a2=zeros(1,size);

将A中大于等于1的元素与B中小于等于0的元素按照对应的位置放入新的向量a中。

for i=1:length(indxA);
    a2(1,indxA(1,i))=A(1,indxA(1,i));
end

for i=1:length(indxB);
    a2(1,indxB(1,i))=B(1,indxB(1,i));
end

这时a2中1和－1代表粒子在A门；2和0代表粒子在B门。将－1换成1，2换成0。

indxA=find(a2==-1);
indxB=find(a2==2);

for i=1:length(indxA);
    a2(1,indxA(1,i))=1;
end

for i=1:length(indxB);
    a2(1,indxB(1,i))=0;
end

通过while loop让粒子跳跃1000次。

%%%%%Example%%%%%

%generate observations% 
size=100000;
a=round(rand(1,size)); %each column is one particle%

j=1000;
while j>1
AA=rand(1,size)<=0.8;
BB=rand(1,size)<=0.3;

A=2*a-AA;
B=2*a-BB;

indxA=find(A>=1);
indxB=find(B<=0);

a2=zeros(1,size);

for i=1:length(indxA);
    a2(1,indxA(1,i))=A(1,indxA(1,i));
end

for i=1:length(indxB);
    a2(1,indxB(1,i))=B(1,indxB(1,i));
end

indxA=find(a2==-1);
indxB=find(a2==2);

for i=1:length(indxA);
    a2(1,indxA(1,i))=1;
end

for i=1:length(indxB);
    a2(1,indxB(1,i))=0;
end

a=a2;
j=j-1;
end

sum(a2)

得出的结果大概在60000左右。也就是说大概能观察到60000个粒子在A门，40000个在B门。

正如前面得出的 π∗=∣∣0.60.4∣∣ ：粒子出现在A门的概率为0.6，出现在B门的概率为0.4。