0%

【随机过程】马氏链的理论与仿真

2014年终总结中,我提到要对这学期学过的数学课中的部分算法进行仿真实现。《数值分析》和《工程优化》这两门数学课里面还有些专门讲算法的,可以用来仿真。在《随机过程》这门课中,几乎全都是公式推导,定理证明,实在难以仿真实现。最后发现,马尔科夫链这一章比较适合仿真,况且先前也写过类似的程序,更重要的是之前有人也问过关于马氏链的Matlab实现问题。关于马氏链的理论原理在这就不作描述,下面直接用程序来实现具体问题的求解。

假设有$9$个状态,其状态转移图如下所示:

图1

根据状态转移图,我们可以得到其一步转移概率矩阵为$P$:

图2


问题一:从状态$1$在$1000$次内转移到状态$9$的概率为多少?

理论值:

由一步转移概率矩阵$P$的性质知,$P$的$n$次幂矩阵$P^n$中的第$i$行第$j$列的元素表示的是从状态$i$经过$n$步转移到状态$j$的概率。注意到,状态$9$为吸收态。因此,在$1000$步以内的任意一步到达状态$9$后,就永远停留在了状态$9$。所以$P$的$1000$次幂矩阵$P^{1000}$中的第$1$行第$9$列就是从状态$1$在$1000$次内转移到状态$9$的概率。

仿真值:

在仿真中,我们设置仿真次数为$100000$次(越大越好),设置初始状态为$1$,在$1000$步中,如果到达状态$9$就记录,然后跳出循环,开始新一轮循环。具体代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
clear all
clc

%一步转移概率矩阵
P=[0.7 0.3 0 0 0 0 0 0 0
0.7 0 0.3 0 0 0 0 0 0
0 0.7 0 0.3 0 0 0 0 0
0 0 0.7 0 0.3 0 0 0 0
0 0 0 0.7 0 0.3 0 0 0
0 0 0 0 0.7 0 0.3 0 0
0 0 0 0 0 0.7 0 0.3 0
0 0 0 0 0 0 0.7 0 0.3
0 0 0 0 0 0 0 0 1];

p=P^1000;
p(1,9)%理论值
num=0;%记录到达状态9的次数
for i=1:100000%实验次数
state=1;%初始状态为1
for j=1:1000%1000步
if state==1
if rand()<0.3
state=2;
end
elseif state==9
state=9;

else
if rand()<0.3
state=state+1;
else
state=state-1;
end
end
if state==9
num=num+1;
break;
end
end
end
num/100000 %仿真值

问题二:从状态$1$到状态$9$平均需要几步?

理论值:

最直观的想法就是:假设从状态$1$到达状态$9$需要$n$步的概率为$p^n$,则平均步数可以表示为$1p^1+2p^2+\cdots+np^n+\cdots$。注意,这里的$p^n$是首达概率,即从状态$1$经过$k$步首次到达状态$9$的概率,也就是前$k-1$次没有到达状态$9$。显然这里的首达概率而不是转移概率矩阵中的概率。这里的首达概率不好求,不过可以矩阵论(可惜没有选这门课)的相关知识来理论求解,最终可以得到的表达式如下:

其中,$x$,$y$表示状态,在本例中,$x=1$,$y=9$,其它符号表示如下:

图3

仿真值:

有时候,理论比较复杂的问题,仿真起来就相对简单,和问题一中的仿真类似,假定仿真次数为$10000$次,设置最大步数为$100000000$,当到达状态$9$时,就记录所需步数并跳出循环,具体的程序实现如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
clear all
clc

P=[0.7 0.3 0 0 0 0 0 0 0
0.7 0 0.3 0 0 0 0 0 0
0 0.7 0 0.3 0 0 0 0 0
0 0 0.7 0 0.3 0 0 0 0
0 0 0 0.7 0 0.3 0 0 0
0 0 0 0 0.7 0 0.3 0 0
0 0 0 0 0 0.7 0 0.3 0
0 0 0 0 0 0 0.7 0 0.3
0 0 0 0 0 0 0 0 1];

sum=0;
p=P^1000;

x=1;
y=9;
P_size=9;%状态个数
Ix=zeros(1,P_size);%行向量
Ix(x)=1;
Iy=zeros(P_size,1);%列向量
Iy(y)=1;
I=eye(P_size);
Ey=I;
Ey(y,y)=0;

average=Ix*(I-P*Ey)^-2*P*Iy%理论值

state=1;%初始状态
num=0;%达到次数
result=0;%需要经过的步数
for k=1:100000%仿真10000次
state=1;
for i=1:100000000%设定最大步数(越大越好)
if state==1
if rand()<0.3
state=2;
end
elseif state==9
state=9;
else
if rand()<0.3
state=state+1;
else
state=state-1;
end
end
if state==9
num=num+1;
result=result+i;%记录总的步数
break;
end
end
end

result/num%计算平均步数
坚持原创技术分享,您的支持将鼓励我继续创作!