马氏链模型及matlab程序

一、用法，用来干什么，什么时候用 二、步骤，前因后果，算法的步骤，公式 三、程序 四、举例
五、前面国赛用到此算法的备注一下
马氏链模型
用来干什么
马尔可夫预测法是应用概率论中马尔可夫链（Markov chain ）的理论和方法来研究分析时间序列的变化规律，并由此预测其未来变化趋势的一种预测技术。

什么时候用
应用马尔可夫链的计算方法进行马尔可夫分析， 主要目的是根据某些变量现在的情 况及其变动趋向，来预测它在未来某特定区间可能产生的变动，作为提供某种决策的依 据。

马尔可夫链的基本原理
我们知道，要描述某种特定时期的随机现象如某种药品在未来某时期的销售情况，比如说第n 季度是畅销还是滞销，用一个随机变量X n 便可以了，但要描述未来所有时期的情况，则需要一系列的随机变量 X 1，X 2，…，X n ，…．称{ X t ，t ∈T ，T 是参数集}为随机过程，{ X t }的取值集合称为状态空间．若随机过程{ X n }的参数为非负整数， X n 为离散随机变量，且{ X n }具有无后效性（或称马尔可夫性），则称这一随机过程为马尔可夫链（简称马氏链）．所谓无后效性，直观地说，就是如果把{ X n }的参数n 看作时间的话，那么它在将来取什么值只与它现在的取值有关，而与过去取什么值无关．
对具有N 个状态的马氏链，描述它的概率性质，最重要的是它在n 时刻处于状态i 下一时刻转移到状态j 的一步转移概率：
N j i n p i X j X P j i n n ,,2,1,)
()|(1 ====+
若假定上式与n 无关，即 ====)()1()0(n p p p j i j i j i ，则可记为j i p （此时，称过程是平稳的），并记
⎪⎪
⎪
⎪
⎪⎭
⎫
⎝
⎛=N N N N N N p p p p p p p p p P
2
12222111211
（1）
称为转移概率矩阵．
转移概率矩阵具有下述性质：
（1）N j i p j i ,,2,1,,
0 =≥．即每个元素非负．
（2）N i p N
j j i ,,2,1,
11
==∑=．即矩阵每行的元素和等于1．
如果我们考虑状态多次转移的情况，则有过程在n 时刻处于状态i ，n +k 时刻转移到状态j 的k 步转移概率：
N j i n p i X j X P k j i n k n ,,2,1,)
()|()
( ====+
同样由平稳性，上式概率与n 无关，可写成)
(k j i p ．记
⎪⎪⎪⎪
⎪⎭
⎫ ⎝⎛=)()(2
)(1
)(2)(22)(21)(1)(12)
(11)
(k N N k N k N k N k k k N k k k p p p p p p p p p P
（2）
称为k 步转移概率矩阵．其中)
(k j i p 具有性质：
N j i p k j
i ,,2,1,,0)
( =≥； N i p N
j k j i ,,2,1,11)
( ==∑=．
一般地有，若P 为一步转移矩阵，则k 步转移矩阵
⎪⎪⎪⎪
⎪⎭
⎫ ⎝⎛=)()(2
)(1
)(2)(22)(21)(1)(12)
(11)
(k N N k N k N k N k k k N k k k p p p p p p p p p P
（3） （2）状态转移概率的估算
在马尔可夫预测方法中，系统状态的转移概率的估算非常重要．估算的方法通常有两种：一是主观概率法，它是根据人们长期积累的经验以及对预测事件的了解，对事件发生的可能性大小的一种主观估计，这种方法一般是在缺乏历史统计资料或资料不全的情况下使用．二是统计估算法，现通过实例介绍如下．
例3 记录了某抗病毒药的6年24个季度的销售情况，得到表1．试求其销售状态的转移概率矩阵．
表1 某抗病毒药24个季度的销售情况
季度
销售状态
季度
销售状态
季度
销售状态
季度
销售状态
1 1 (畅销) 7 1(畅销) 13 1(畅销) 19 2(滞销)
2 1(畅销) 8 1(畅销) 14 1(畅销) 20 1(畅销)
3 2(滞销) 9 1(畅销) 15 2(滞销) 21 2(滞销)
4 1(畅销) 10 2(滞销) 16 2(滞销) 22 1(畅销)
5 2(滞销) 11 1(畅销) 17 1(畅销) 23 1(畅销) 6
2(滞销)
12
2(滞销)
18
1(畅销)
24
1(畅销)
分析表中的数据，其中有15个季度畅销，9个季度滞销，连续出现畅销和由畅销转入滞销以及由滞销转入畅销的次数均为7，连续滞销的次数为2．由此，可得到下面的市场状态转移情况表（表2）．
表2 市场状态转移情况表
现计算转移概率．以频率代替概率，可得连续畅销的概率：
117
0.5151
p =
=
=-连续出现畅销的次数出现畅销的次数
分母中的数为15减1是因为第24季度是畅销，无后续记录，需减1．
同样得由畅销转入滞销的概率：
127
0.5151
p =
=
=-畅销转入滞销的次数出现畅销的次数
滞销转入畅销的概率：
217
0.789
p =
=
=滞销转入畅销的次数出现滞销的次数
连续滞销的概率：
222
0.229
p =
=
=连续滞销的次数出现滞销的次数
综上，得销售状态转移概率矩阵为：
⎪⎪⎭⎫ ⎝
⎛=⎪⎪⎭⎫ ⎝⎛=22.078.05.05.022211211p p
p p P 从上面的计算过程知，所求转移概率矩阵P 的元素其实可以直接通过表2中的数字计算而得到，即将表中数分别除以该数所在行的数字和便可：
777
11+=
p 777
12+=
p 27721+=
p 772
22+=
p Matlab 程序：
format rat clc
a=[ 1 1 2 1 2 2 1 1 1 2 1 2,1 1 2 2 1 1 2 1 2 1 1 1]; for i=1:2 for j=1:2
f(i,j)=length(findstr([i j],a)); end end f
ni=(sum(f'))' for i=1:2
p(i,:)=f(i,:)/ni(i); end p
由此，推广到一般情况，我们得到估计转移概率的方法：假定系统有m 种状态S 1，S 2，…，
S m ，根据系统的状态转移的历史记录，得到表3的统计表格，以j i p
ˆ表示系统从状态i 转移到状态j 的转移概率估计值，则由表3的数据计算估计值的公式如下：
表3 系统状态转移情况表
（3）带利润的马氏链
在马氏链模型中，随着时间的推移，系统的状态可能发生转移，这种转移常常会引起某种经济指标的变化．如抗病毒药的销售状态有畅销和滞销两种，在时间变化过程中，有时呈连续畅销或连续滞销，有时由畅销转为滞销或由滞销转为畅销，每次转移不是盈利就是亏本．假定连续畅销时盈r 11元，连续滞销时亏本r 22元，由畅销转为滞销盈利r 12元，由滞销转为畅销盈利r 21元，这种随着系统的状态转移，赋予一定利润的马氏链，称为有利润的马氏链．对于一般的具有转移矩阵
⎪⎪⎪
⎪
⎪⎭⎫
⎝
⎛=N N N N N N p p p p p p p p p P
2
12222111211
的马氏链，当系统由i 转移到j 时，赋予利润r ij （i ，j =1，2，…，N ），则称
⎪⎪⎪⎪
⎪⎭
⎫
⎝
⎛=N N N N N N r r r r r r r r r R
2
122221112
11 （5） 为系统的利润矩阵，r ij ＞0称为盈利，r ij ＜0称为亏本，r ij = 0称为不亏不盈．
随着时间的变化，系统的状态不断地转移，从而可得到一系列利润，由于状态的转移是随机的，因而一系列的利润是随机变量，其概率关系由马氏链的转移概率决定．例如从抗病毒药的销售状态的转移矩阵，得到一步利润随机变量)
1(1x 、)
1(2x 的概率分布分别为：
其中 p 11+ p 12 = 1 ，p 21+ p 22 = 1．
如果药品处于畅销阶段，即销售状态为i =1，我们想知道，经过n 个季度以后，期望获得的利润是多少？为此，引入一些计算公式．
首先，定义)
(n i v 为抗病毒药现在处于)2,1(=i i ，经过n 步转移之后的总期望利润，
则
一步转移的期望利润为：
∑==+==
2
1
2211)1()
1()
(j j i j i i i i i i i p r p r p r x E v
其中)()
1(i x E 是随机变量)
1(i x 的数学期望．
二步转移的期望利润为：
∑=+=
+
++
==
2
1
)1(2
)
1(221
)
1(11)2()
2(][][][)
(j j i j j i i i i i i i p v r p v r p v r x E v
其中随机变量)
2(i
x （称为二步利润随机变量）的分布为：
2,1,
)()1()2(==+=j p v r x P j i j j i i
例如，若
⎪⎪⎭⎫ ⎝⎛=6.04.05.05.0P ， ⎪⎪
⎭
⎫
⎝⎛-=7339R
则抗病毒药销售的一步利润随机变量：
抗病毒药畅销和滞销时的一步转移的期望利润分别为：
65.035.09)(12121111)
1(1)1(1=⨯+⨯=+==p r p r x E v 36.074.03)(22222121)1(2)1(2-=⨯-⨯=+==p r p r x E v
二步利润随机变量为：
抗病毒药畅销和滞销时的二步转移的期望利润分别为：
12)
1(21211)1(111)2(1)2(1][][)(p v r p v r x E v +++==
5.75.0)33(5.0)69(=⨯-+⨯+=
22)1(22221)1(121)2(2)2(2][][)(p v r p v r x E v +++==
4.26.0)37(4.0)63(-=⨯--+⨯+=
一般地定义k 步转移利润随机变量),2,1()
(N i x k i
=的分布为：
N j p v r x P j
i k j j i k i ,2,1)()1()(==+=-
则系统处于状态i 经过k 步转移后所得的期望利润)
(k i
v 的递推计算式为：
j i k j N
j j i k i k i p v r x E v )()
()1(1)()
(-=+==
∑
∑∑
∑=-=-=+=
+=
N
j j i k j i N
j j i k j N
j j i j i p v v p v p r 1
)1()1(1
)
1(1
（6）
当k =1时，规定边界条件0)
0(=i
v ．
称一步转移的期望利润为即时的期望利润，并记
N i q v i i ,2,1,)1(==．
可能的应用题型
题型一、市场占有率预测
例题1在购买该药的总共1000家对象（购买力相当的医院、药店等）中，买A 、B 、C 三药厂的各有400家、300家、300家，预测A 、B 、C 三个厂家生产的某种抗病毒药在未来的市场占有情况。

顾客订货情况如下表5：
表5 顾客订货情况表
下季度订货情况
合计
来自
A B C
A 160 120 120 400
B 180 90 30 300
C 180 30 90 300
合计520 240 240 1000
模型建立与求解
一、问题分析
目前的市场占有情况为：在购买该药的总共1000家对象（购买力相当的医院、药店等）中，买A、B、C三药厂的各有400家、300家、300家，那么A、B、C三药厂目前的市场占有份额分别为：40%、30%、30%．称（0.4，0.3，0.3）为目前市场的占有分布或称初始分布．此外，我们需要查清使用对象的流动情况。

流动情况的调查可通过发放信息调查表来了解顾客以往的资料或将来的购买意向，也可从下一时期的订货单得出。

由题已知顾客订货情况如下表5
表5 顾客订货情况表
下季度订货情况合计
来自
A B C
A 160 120 120 400
B 180 90 30 300
C 180 30 90 300
合计520 240 240 1000
二、模型的建立
2.1模型构建
假定在未来的时期内，顾客相同间隔时间的流动情况不因时期的不同而发生变化，以1、2、3分别表示顾客买A、B、C三厂家的药这三个状态，以季度为模型的步长（即转移一步所需的时间），那么根据表5，我们可以得模型的转移概率矩阵：
⎪⎪⎪⎭
⎫
⎝⎛=⎪⎪⎪⎪⎪⎪
⎭
⎫
⎝⎛=⎪⎪⎪⎭⎫ ⎝⎛=3.01
.06.01.03.06.03.03.04.030090300
30300
18030030300903001804001204001204001603332
31
232221
131211p p p p p p p p p P 矩阵中的第一行（0.4，0.3，0.3）表示目前是A 厂的顾客下季度有40%仍买A 厂的药，转为买B 厂和C 厂的各有30%．同样，第二行、第三行分别表示目前是B 厂和C 厂的顾客下季度的流向．
由P 我们可以计算任意的k 步转移矩阵，如三步转移矩阵：
⎪⎪⎪⎭
⎫
⎝
⎛=⎪⎪⎪⎭⎫ ⎝
⎛==252.0244
.0504.0244.0252.0504
.0252.0252
.0496.03.01
.06.01.03.06
.03.03.04
.03
3
)
3(P P 从这个矩阵的各行可知三个季度以后各厂家顾客的流动情况．如从第二行（0.504，0.252，0.244）知，B 厂的顾客三个季度后有50.4%转向买A 厂的药，25.2%仍买B 厂的，24.4%转向买C 厂的药．
设),,()
(3)(2)(1)
(k k k k p p p S
=表示预测对象k 季度以后的市场占有率，初始分布则为
),,()
0(3)0(2)0(1)0(p p p S =，市场占有率的预测模型为
P S P S S k k k ⋅=⋅=-)1()0()( （7）
已知)3.0,3.0,4.0()
0(=S
，由此，我们可预测任意时期A 、B 、C 三厂家的市场占有
率．例如，三个季度以后的预测值为：
⎪⎪⎪⎭
⎫
⎝
⎛=⋅==252.0244
.0504.0244.0252.0504
.0252.0252.0496
.0)3.03.04.0(),,(3
)0()
3(3
)3(2)3(1)3(P S p p p S )2496.02496.05008.0(=
大致上，A 厂占有一半的市场，B 厂、C 厂各占四分之一．
模型（7）可推广到N 个状态的情形：
k
N N N N N N N
k
k k p p p p p p p p p p
p p P S
P S
S
⎪⎪⎪⎪
⎪⎭
⎫
⎝⎛===-
2
1
2222111211)0()0(2)
0(1)
0()
1()
(),,(（8）
如果我们按公式（7）继续逐步求A 、B 、C 三家的市场占有率，会发现，当k 大到一定的程度，S (k )
将不会有多少改变，即有稳定的市场占有率，设其稳定值为),,(321p p p S =，满足1321=++p p p ．
事实上，如果市场的顾客流动趋向长期稳定下去，则经过一段时期以后的市场占有率将会出现稳定的平衡状态，即顾客的流动，不会影响市场的占有率，而且这种占有率与初始分布无关．如何求出这种稳定的市场占有率呢？
2.2模型求解
以A 、B 、C 三家的情况为例，当市场出现平衡状态时，从公式（7）可得方程S = S P ，即
⎪⎪⎪⎭
⎫
⎝
⎛=3.01
.06.01.03.06
.03.03.04
.0),,(),,(321321p p p p p p 由此得
⎪⎩⎪
⎨⎧++=++=++=3213
32123
2113.01.03.01.03.03.06.06.04.0p
p p p p p p p p p p p 经整理，并加上条件1321=++p p p ，得
⎪⎪⎩⎪⎪⎨
⎧=++=-+=+-=++-1
p p p 00.7p 0.1p 0.3p 00.1p 0.7p 0.3p 0
0.6p 0.6p 0.6p 3213213
21321 0.4 0.3 0.3，0.6 0.3 0.1,0.6 0.1 0.3
上方程组是三个变量四个方程的方程组，在前三个方程中只有二个是独立的，任意删去一个，从剩下的三个方程中，可求出唯一解：
5.01=p ， 25.02=p ， 25.03=p
这就是A 、B 、C 三家的最终市场占有率．
一般N 个状态的稳定市场占有率（稳态概率）),,(21N p p p S =可通过解方程组
⎪⎪⎪⎩⎪⎪⎪
⎨⎧=⎪⎪⎪
⎪⎪⎭
⎫
⎝⎛=∑=1),,(),,(1
2
1
22221112112121N
k k N N N N N N N
N p p p p p p p p p p p
p p p p p
（9） 求得，而（9）的前N 个方程中只有N -1个是独立的，可任意删去一个。

MATLAB 程序：
format rat
p=[0.4 0.3 0.3，0.6 0.3 0.1,0.6 0.1 0.3]; a=[p'-eye(3);ones(1,3)]; b=[zeros(3,1);1]; p_limit=a\b
题型二、期望利润预测
企业追逐市场占有率的真正目的是使利润增加，因此，竞争各方无论是为了夺回市场份额，还是为了保住或者提高市场份额，在制订对策时都必须对期望利润进行预测．
预测主要分两步进行：①市场统计调查．首先调查销路的变化情况，即查清由畅销到滞销或由滞销到畅销，连续畅销或连续滞销的可能性是多少．其次统计出由于销路的变化，获得的利润和亏损情况．②建立数学模型，列出预测公式进行预测．
例如，通过市场调查，我们得到如下的销路转移表（表6）和利润变化表（表7）．由此，我们来建立数学模型．
表6 销路转移表
销路转移表说明连续畅销的可能性为50%，由畅销转入滞销的可能性也是50%，由滞销到畅销为40%，连续滞销的可能性为60%．利润表说明的是连续畅销获利900万元，由畅销
到滞销或由滞销到畅销均获利300万元，连续滞销则亏损700万元．从而得到销售状态的转移矩阵P 和利润矩阵R 分别为：
表7 利润变化表（单位：百万元）
⎪⎪⎭
⎫
⎝⎛=⎪⎪⎭⎫ ⎝⎛=6.04.05.05.022211211p p
p p P ⎪⎪⎭
⎫
⎝⎛-=⎪⎪⎭⎫ ⎝⎛=733922211211
r r
r r R P 和R 便构成一个有利润的马氏链．由前面所述的基本原理及公式（6）得下面的预测公式：
即时期利润：2,12
1
)
1(==
=
∑=i p r v q j j
i j i i i
k
步以后的期望利润：2,12
1
)1(2
1
)1(2
1
)
(=+=+=
∑∑
∑=-=-=i p v q p v p r v j j
i k j i j j i k j j j i j i k i
将调查数据代入上公式则可预测各时期的期望利润值．如：
65.035.091=⨯+⨯=q 36.074.032-=⨯-⨯=q
由此可知，当本季度处于畅销时，在下一季度可以期望获得利润600万元；当本季度处于滞销时，下一季度将期望亏损300万元．
同样算得： 5.7)
2(1
=v ， 4.2)
2(2-=v
55.8)3(1=v ， 44.1)3(2-=v
由此可预测本季度处于畅销时，两个季度后可期望获利750万元，三个季度后可期望获利855万元；当本季度处于滞销时，两个季度后将亏损240万元，三个季度后亏损144万元。

MATLAB 程序：。