数学建模算法第17章 马氏链模型Word格式.docx
《数学建模算法第17章 马氏链模型Word格式.docx》由会员分享，可在线阅读，更多相关《数学建模算法第17章 马氏链模型Word格式.docx（27页珍藏版）》请在冰点文库上搜索。
j|ξn=i}
（2）
则称{ξn,n=1,2,L}为一个马尔可夫链（简称马氏链），
（2）式称为马氏性。
事实上，可以证明若等式
（2）对于m=1成立，则它对于任意的正整数m也成立。
因此，只要当m=1时
（2）式成立，就可以称随机序列{ξn,n=1,2,L}具有马氏性，即{ξn,n=1,2,L}是一个马尔可夫链。
定义3设{ξn,n=1,2,L}是一个马氏链。
如果等式
（2）右边的条件概率与n无关，即
P{ξn+m=
j|ξn=i}=pij（m）
（3）
则称{ξn,n=1,2,L}为时齐的马氏链。
称pij（m）为系统由状态i经过m个时间间隔（或m步）转移到状态j的转移概率。
（3）称为时齐性。
它的含义是：
系统由状态i到状态j的转移概率只依赖于时间间隔的长短，与起始的时刻无关。
本章介绍的马氏链假定都是时齐的，因此省略“时齐”二字。
2.2转移概率矩阵及柯尔莫哥洛夫定理
对于一个马尔可夫链{ξn,n=1,2,L}，称以m步转移概率pij（m）为元素的矩阵P（m）=（pij（m））为马尔可夫链的m步转移矩阵。
当m=1时，记P
（1）=P称为马尔可夫链的一步转移矩阵，或简称转移矩阵。
它们具有下列三个基本性质：
（i）对一切i,j∈E，0≤pij（m）≤1；
（ii）对一切i∈E，∑pij（m）=1；
j∈E
⎧
1,
（iii）对一切i,j∈E，pij（0）=δij=⎨
当i=j时
。
⎩0，当i≠j时
当实际问题可以用马尔可夫链来描述时，首先要确定它的状态空间及参数集合，然后确定它的一步转移概率。
关于这一概率的确定，可以由问题的内在规律得到，也可以由过去经验给出，还可以根据观测数据来估计。
例4某计算机机房的一台计算机经常出故障，研究者每隔15分钟观察一次计算
机的运行状态，收集了24小时的数据（共作97次观察）。
用1表示正常状态，用0表示不正常状态，所得的数据序列如下：
111001*********0011110111111001111111110001101101
111011*********101110111101111110011011111100111
解设Xn（n=1,L,97）为第n个时段的计算机状态，可以认为它是一个时齐马氏链，状态空间E={0,1}，编写如下Matlab程序：
a1='
1110010011111110011110111111001111111110001101101'
;
a2='
111011011010111101110111101111110011011111100111'
a=[a1a2];
f00=length（findstr（'
00'
a））f01=length（findstr（'
01'
a））f10=length（findstr（'
10'
a））f11=length（findstr（'
11'
a））
或者把上述数据序列保存到纯文本文件data1.txt中，存放在Matlab下的work子目录中，编写程序如下：
clc,clear
formatratfid=fopen（'
data1.txt'
'
r'
）;
a=[];
while（~feof（fid））a=[afgetl（fid）];
end
fori=0:
1
forj=0:
s=[int2str（i）,int2str（j）];
f（i+1,j+1）=length（findstr（s,a））;
fs=sum（f'
fori=1:
2
f（i,:
）=f（i,:
）/fs（i）;
endf
求得96次状态转移的情况是：
0→0，8次；
1→0，18次；
0→1，18次；
1→1，52次，
因此，一步转移概率可用频率近似地表示为
p00
p
=P{X
n+1
=0|Xn
=1|X
=0}≈
8=4
8+1813
18=9
01n+1n
p10
=0|Xn
=1}≈
18+5235
p11
=1|Xn
52
18+52
=26
35
例5设一随机系统状态空间E={1,2,3,4}，记录观测系统所处状态如下：
4321431123
2123443311
1332122244
2323112431
若该系统可用马氏模型描述，估计转移概率pij。
解首先将不同类型的转移数nij统计出来分类记入下表
i→j转移数nij
3
4
行和ni
10
11
7
各类转移总和∑∑nij等于观测数据中马氏链处于各种状态次数总和减1，而行和ni是
ij
系统从状态i转移到其它状态的次数，nij是由状态i到状态j的转移次数，则pij的估
nij
n
计值pij=。
计算得
i
⎡2/52/51/101/10⎤
ˆ⎢3/112/114/112/11⎥
P=⎢
⎢4/11
4/11
⎦
2/11
⎥
1/11⎥
⎣
⎢1/7
4/7
2/7⎥
Matlab计算程序如下：
formatratclc
a=[4
3...
1...
4...
1];
fori=1:
4forj=1:
f（i,j）=length（findstr（[ij],a））;
end
ni=（sum（f'
））'
p（i,:
）/ni（i）;
例6（带有反射壁的随机徘徊）如果在原点右边距离原点一个单位及距原点s（s>
1）个单位处各立一个弹性壁。
一个质点在数轴右半部从距原点两个单位处开始随机徘徊。
每次分别以概率p（0<
p<
1）和q（q=1-p）向右和向左移动一个单位；
若在+1处，则
以概率p反射到2，以概率q停在原处；
在s处，则以概率q反射到s-1，以概率p停在原处。
设ξn表示徘徊n步后的质点位置。
{ξn,n=1,2,L}是一个马尔可夫链，其状态空间E={1,2,L,s}，写出转移矩阵P。
⎨0
解P{ξ=i}=⎧1,
当i=2时
p1j
psj
⎧q,
⎨
=⎪p,
⎪⎩0,
⎧p,
=⎪q,
⎩0，当i≠2时当j=1时
当j=2时其它
当j=s时当j=s-1时其它
pij
当j-i=1时
当j-i=-1时（i=2,3,L,s-1）
其它
因此，P为一个s阶方阵，即
L
q
⎡q0⎤
⎢⎥
⎢00⎥
P=⎢⎥。
⎢LL⎥
⎢0p⎥
⎣⎦
定理1（柯尔莫哥洛夫—开普曼定理）设{ξn,n=1,2,L}是一个马尔可夫链，其
状态空间E={1,2,L}，则对任意正整数m,n有
pij（n+m）=∑pik（n）pkj（m）
k∈E
其中的i,j∈E。
定理2设P是一个马氏链转移矩阵（P的行向量是概率向量），P（0）是初始分布行向量，则第n步的概率分布为
P（n）=P（0）Pn。
例7若顾客的购买是无记忆的，即已知现在顾客购买情况，未来顾客的购买情况不受过去购买历史的影响，而只与现在购买情况有关。
现在市场上供应A、B、C三个不同厂家生产的50克袋状味精，用“ξn=1”、“ξn=2”、“ξn=3”分别表示“顾客第n次购买A、B、C厂的味精”。
显然，{ξn,n=1,2,L}是一个马氏链。
若已知第一次顾客购买三个厂味精的概率依次为0.2，0.4，0.4。
又知道一般顾客购买的倾向由表2给出。
求顾客第四次购买各家味精的概率。
表2
下次购买
A
B
C
上次购买
AB
0.8
0.5
0.1
0.3
0.4
0.2
解第一次购买的概率分布为
P
（1）=[0.20.40.4]
⎡0.8
⎢
转移矩阵P=⎢0.5
⎢⎣0.5
0.1⎤
0.4⎥
0.2⎥⎦
则顾客第四次购买各家味精的概率为
P（4）=P
（1）P3=[0.7004
0.136
0.1636]。
2.3转移概率的渐近性质—极限概率分布
现在我们考虑，随n的增大，Pn是否会趋于某一固定向量？
先考虑一个简单例子：
⎡0.5
转移矩阵P=⎢0.7
0.5⎤
0.3⎥，当n→+∞时，
⎡7
Pn→⎢12
⎣12
5⎤
12⎥
5⎥
12⎦
又若取u=⎡7
5⎤，则uP=u，uT为矩阵PT的对应于特征值λ=1的特征（概
率）向量，u也称为P的不动点向量。
哪些转移矩阵具有不动点向量？
为此我们给出正则矩阵的概念。
定义４一个马氏链的转移矩阵P是正则的，当且仅当存在正整数k，使Pk的每一元素都是正数。
定理３若P是一个马氏链的正则阵，那么：
（i）P有唯一的不动点向量W，W的每个分量为正。
（ii）P的n次幂Pn（n为正整数）随n的增加趋于矩阵W，W的每一行向量均等于不动点向量W。
例8信息的传播一条新闻在a1,a2,L,an,L等人中间传播，传播的方式是a1传
给a2，a2传给a3，…如此继续下去，每次传播都是由ai传给ai+1。
每次传播消息的失真概率是p，0<
1，即ai将消息传给ai+1时，传错的概率是p，这样经过长时间传播，第n个人得知消息时，消息的真实程度如何？
设整个传播过程为随机转移过程，消息经过一次传播失真的概率为p，转移矩阵
假真
假⎡1-pp⎤
P=⎢⎥
真⎢⎣p
1-p⎥⎦
P是正则矩阵。
又设V是初始分布，则消息经过n次传播后，其可靠程度的概率分布为V⋅Pn。
一般地，设时齐马氏链的状态空间为E，如果对于所有i,j∈E，转移概率pij（n）
存在极限
→∞
limpij（n）=πj，（不依赖于i）
或
⎡π1π2
π
⎢12
P（n）=Pn⎯⎯⎯→⎢LL
Lπj
LL
L⎤
L⎥
L⎥，
（n→∞）
LπjL⎥
⎢⎣L
LLL
L⎥⎦
则称此链具有遍历性。
又若∑πj=1，则同时称π=（π1,π2,L）为链的极限分布。
j
下面就有限链的遍历性给出一个充分条件。
定理4设时齐（齐次）马氏链{ξn,n=1,2,L}的状态空间为E={a1,L,aN}，P=（pij）是它的一步转移概率矩阵，如果存在正整数m，使对任意的ai,aj∈E，都有
pij（m）>
0，i,j=1,2,L,N
则此链具有遍历性；
且有极限分布π=（π1,L,πN），它是方程组
N
的满足条件
π=πP
或即πj=∑πipij，j=1,L,N
i=1
πj>
0，∑πj=1
j=1
的唯一解。
例9根据例7中给出的一般顾客购买三种味精倾向的转移矩阵，预测经过长期的多次购买之后，顾客的购买倾向如何？
解这个马氏链的转移矩阵满足定理4的条件，可以求出其极限概率分布。
为此，解下列方程组：
⎧p1=0.8p1+0.5p2+0.5p3
⎪
⎪2=0.1p1+0.1p2+0.3p3
⎪3=0.1p1+0.4p2+0.2p3
⎪⎩p1+p2+p3=1
编写如下的Matlab程序：
formatrat
p=[0.80.10.1;
0.50.10.4;
0.50.30.2];
a=[p'
-eye（3）;
ones（1,3）];
b=[zeros（3,1）;
p_limit=a\b
或者利用求转移矩阵P的转置矩阵PT的特征值1对应的特征（概率）向量，求得极限概率。
编写程序如下：
p=sym（p'
[x,y]=eig（p）fori=1:
x（:
i）=x（:
i）/sum（x（:
i））;
x
51113
求得p1=7，p2=84，p3=84。
这说明，无论第一次顾客购买的情况如何，经过长期多次购买以后，A厂产的味
5
精占有市场的
，B,C两厂产品分别占有市场的11，13。
8484
2.4吸收链
马氏链还有一种重要类型—吸收链。
若马氏链的转移矩阵为
1⎡0.3
23
0.30
0.4⎤
2⎢0.2
P=3⎢0
4⎢0
0.3⎥，
P的最后一行表示的是，当转移到状态4时，将停留在状态4，状态4称为吸收状态。
如果马氏链至少含有一个吸收状态，并且从每一个非吸收状态出发，都可以到达某
个吸收状态，那么这个马氏链被称为吸收链。
具有r个吸收状态，s（s=n-r）个非吸收状态的吸收链，它的n⨯n转移矩阵的标准形式为
⎢R
S
P=⎡IrO⎤
（4）
其中Ir为r阶单位阵，O为r⨯s零阵，R为s⨯r矩阵，S为s⨯s矩阵。
从（4）得
⎢Q
Pn=⎡Ir
O⎤
Sn⎥
（5）
（5）式中的子阵Sn表示以任何非吸收状态作为初始状态，经过n步转移后，处于s个非吸收状态的概率。
在吸收链中，令F=（I-S）-1，则F称为基矩阵。
对于具有标准形式（即（4）式）转移矩阵的吸收链，可以证明以下定理：
定理5吸收链的基矩阵F中的每个元素，表示从一个非吸收状态出发，过程到达每个非吸收状态的平均转移次数。
定理6设N=FC，F为吸收链的基矩阵，C=[11L
1]T，则N的每个
元素表示从非吸收状态出发，到达某个吸收状态被吸收之前的平均转移次数。
定理7设B=FR=（bij），其中F为吸收链的基矩阵，R为（4）式中的子阵，则bij表示从非吸收状态i出发，被吸收状态j吸收的概率。
例10智力竞赛问题甲、乙两队进行智力竞赛。
竞赛规则规定：
竞赛开始时，甲、乙两队各记2分，在抢答问题时，如果甲队赢得1分，那么甲队的总分将增加1
分，同时乙队总分将减少1分。
当甲（或乙）队总分达到4分时，竞赛结束，甲（或乙）获胜。
根据队员的智力水平，知道甲队赢得1分的概率为p，失去1分的概率为1-p，
求：
（i）甲队获胜的概率是多少？
（ii）竞赛从开始到结束，分数转移的平均次数是多少？
（iii）甲队获得1、2、3分的平均次数是多少？
分析甲队得分有5种可能，即0、1、2、3、4，分别记为状态a0,a1,a2,a3,a4，其中a0和a4是吸收状态，a1,a2和a3是非吸收状态。
过程是以a2作为初始状态。
根据
甲队赢得1分的概率为p，建立转移矩阵：
a0a1a2a3a4
a0⎡10000⎤
a⎢1-p0p00⎥
1⎢
P=a2⎢0
1-p0
p0⎥
（6）
a3⎢001-p0
a4⎢⎣00
将（6）式改记为标准形式：
P=⎡I2O⎤
001⎥⎦
其中
⎡1-p0⎤
⎡0p0⎤
R=⎢0
0⎥，S=⎢1-p0
p⎥，
⎣⎢0
计算
p⎥⎦
1-p
0⎥⎦
⎡1-pqpp2⎤
F=（I
-S）-1=
1⎢q1p⎥
其中q=1-p。
1-2pq⎢
⎣q2
q1-pq⎥⎦
因为a2是初始状态，根据定理5，甲队获得1，2，3分的平均次数为1-2pq，
1-2pq
，。
又
N=FC=
q1-pq⎥⎦
=1
[1+2p2
21+2p2]
根据定理6，以a2为初始状态，甲队最终获胜的分数转移的平均次数为
又因为
⎡（1-pq）p
p3⎤
B=FR=
1⎢q2
p2⎥
⎣q3
（1-pq）p⎥⎦
p2
根据定理7，甲队最后获胜的概率b22=1-2pq。
Matlab程序如下：
symspqr=[q,0;
0,0;
0,p];
s=[0,p,0;
q,0,p;
0,q,0];
f=（eye（3）-s）^（-1）;
f=simple（f）n=f*ones（3,1）;
n=simple（n）b=f*r;
b=simple（b）
3马尔可夫链的应用
应用马尔可夫链的计算方法进行马尔可夫分析，主要目的是根据某些变量现在的情况及其变动趋向，来预测它在未来某特定区间可能产生的变动，作为提供某种决策的依据。
例11（服务网点的设置问题）为适应日益扩大的旅
升级会员 