马尔科夫模型应用

普大帝

原发表于数学中国微信公众号，关注数学中国微信公众号可查看更多。请用电脑查看，手机版图片可能无法正常显示。

​

清晨，你醒了，睁开双眼的你面临3个选择：继续睡一会，躺在床上听听音乐、看看手机，起床准备学习.勤奋的你选择了去学习.到了下午，你可能决定放松身心，娱乐一会.晚上你又有其他的选择.如果把每天分时段进行的活动看作是一个随机过程，而你做出的决定不取决于你之前干过什么，比如你决定下午玩耍一下和早上有没有看过书没有关系，那么这就是一个特殊的随机过程——马尔科夫过程.

把所有可能进行的活动总和记为状态空间E，从某时刻开始的时间段记为n，则上述例子可以描述为在状态空间E下的马尔科夫过程{Xn,n≥0}.

马尔科夫过程的特性在于未来的演变不依赖于它过去的演变.例如，明天是否会下雨不依赖于昨天是否下雨，这种性质被称作无后效性.

现实中很多问题都可以看作马尔科夫过程，如布朗运动、传染病爆发过程、车站候车人流量等.马尔科夫模型也在网站流量分析、教学质量评估、股票期权等方面得到了广泛的应用.下面给出几个马尔科夫模型适用的典型问题.

1 疾病健康问题

问题描述 人的健康状态随着时间的推移会随机地发生转变，保险公司要对投保人未来的健康状态进行估计，以制订保险金和理赔金的数额.人的健康状况分为健康和患病两种状态，设对特定年龄段的人，今年健康、明年保持健康状态的概率为0.8，今年患病、明年转为健康状态的概率为0.7.若某人投保时健康，问10年后他仍处于健康状态的概率是多少?

思路分析 首先需要求解出状态转移概率矩阵，由马尔科夫链的基本方程，根据给定的a（0），可以预测 a（n），这样就可以求出任意时间点的状态概率了.

本模型满足马尔科夫链的基本要求，即该人某年健康或患病的概率只与其前一年健康或患病状态有关，与再前面各年份健康情况无关.

模型建立 设Xn表示第n年投保人身体所处状态，记状态[插图],则Xn是时间状态均离散的马尔科夫链，其中，状态概率ai（n）=P（Xn=i）,i=1,2; n=0,1,…，转移概率pij=P（Xn+1=j|Xn=i）,i,j=1,2; n=0,1,….

健康状态与状态转移示意如图2.5所示.

图2.5 健康状态与状态转移示意图

转移方程如下

给定a（0）→预测 a（n）,n=1,2…

显然，转移概率矩阵为模型求解Xn+1只取决于Xn和Pij，与Xn-1,…无关，状态转移具有无后效性，由马尔科夫链的基本方程，根据给定a（0），可以预测 a（n）,n=1,2…

（1）设投保时健康，即a1（0）=1,a2（0）=0，由此得到状态概率a（n）如表2.3所示.

a1（1）=a1（0）p11+a2（0）p21=1×0.8+0×0.7=0.8

a2（1）=a1（0）p12+a2（0）p22=1×0.2+0×0.3=0.2

……

a1（2）=a1（1）p11+a2（1）p21=0.8×0.8+0.2×0.7=0.78

a2（2）=a1（1）p12+a2（1）p22=0.8×0.2+0.2×0.3=0.22

……

表2.3 投保时健康得到的状态概率

（2）设投保时患病，即a1（0）=0,a2（0）=1，由此类似得到状态概率a（n）如表2.4所示.

表2.4 投保时患病得到的状态概率

由此可知，当n→∞时，状态概率趋于稳定值，稳定值与初始状态无关.

事实上，由前述一般例子，a=0.2,b=0.7，从而取状态0的极限概率为

，取状态1的极限概率为，与系统所处的初始状态无关.

也可用如下方式计算极限分布w满足wP=w，有

，从而得w=（7/9,2/9）.

若患病病人平均理赔金额为2 000元，投保期为10年，投保人数为10人，则由稳态概率可近似计算每年每人的投保金额x（未考虑利息理论），由10×10x=2 000×2/9×10，从而得x=44.4（元）.

2 疾病健康死亡问题

问题描述 如果人的状态分为健康、疾病和死亡3种状态，记Xn=1表示n年后投保人身体健康，Xn=2表示投保人患病，Xn=3表示投保人因疾病死亡.若三种状态的转换概率如图2.6所示，问n年后，该人处于3种状态的概率分别为多少?

思路分析 注意，在这个模型中多了一个状态Xn=3，即死亡，一旦投保人从状态1或状态2转移到此状态，则其就永久停在这个状态了.这是和前一问题的最大不同.因而状态3转到自身的转移概率永远是1，不会再变回其他状态了.再根据图2.6，可以知道整个状态转移矩阵了.

图2.6 3种状态的转换概率

模型建立 由于这里的状态3是吸收状态，也就是说，到这个状态后，将不再会改变为其他状态，所以模型和前面有一定的不同.其转移矩阵为

由如下转移公式即可求得n+1个时刻，该人在健康、患病或者死亡状态的概率.

a1（n+1）=a1（n）p11+a2（n）p21+a3（n）p31，

a2（n+1）=a1（n）p12+a2（n）p22+a3（n）p32，

a3（n+1）=a1（n）p13+a2（n）p23+a3（n）p33.

模型求解 用MATLAB编程的一个程序如下.

n=input（＇n=＇）

A=zeros（3,n+1）；

A（1,1）=input（＇a01=＇）；

A（2,1）=input（＇a02=＇）；

A（3,1）=1-A（1,1）-A（2,1）；

for i=1:n

　　A（1,i+1）=0.8*A（1,i）+0.65*A（2,i）+0*A（3,i）；

　　A（2,i+1）=0.18*A（1,i）+0.25*A（2,i）+0*A（3,i）；

　　A（3,i+1）=0.02*A（1,i）+0.1*A（2,i）+1*A（3,i）；

end

A

从而可得每一时刻的状态与状态转移概率，如表2.5所示.

表2.5 同一时刻状态与状态转移概率

从表2.5看出，无论初始状态是哪种情形，当投保年份越来越多时，最终投保人都会转到状态3，即投保人死亡.一旦a1（k）= a2（k）=0,a3（k）=1，则对于n＞k，恒有a1（n）=0,a2（n）=0,a3（n）=1，即从状态3不再会转移到其他状态.该马尔科夫链被称为吸收链.

3 汽车工况问题

马尔科夫链在经济、社会、生态、遗传等许多领域中有着广泛的应用.值得指出的是，虽然它是解决随机转移过程问题的工具，但是一些确定性系统的状态转移问题也能用马尔科夫链模型处理.下面的实例就是马尔科夫链在汽车工况研究中的应用.

问题描述汽车工厂要了解某一类重型汽车的行驶状况（工况），来研究如何降低汽车油耗.但重型汽车在公路上行驶时间一般很长，且因为路况复杂，速度变化很不均匀，所以很难对重型汽车实际行驶状况进行分析.需要模拟并在实验室重现能够代替实际汽车道路行驶的工况.如何构造一定时间的汽车行驶工况，且能代表重型汽车长时间的实际行驶状况，从而可以在实验室对汽车发动机进行各种配置以找到最优配置?

思路分析 重型汽车在公路上行驶常常会因为路况、避让等原因造成速度间歇性地不均匀变化，需要消除这些不必要的行驶状态，为行驶实验提供稳定的行驶状态，但同时还需要让重组的行驶状态能代替汽车实际行驶状况，这是非常关键的一点.

重型汽车行驶的速度变化图是不规则的曲线.首先将这条曲线按照一定规则切割成无数小段，再通过提取和重组，构建出具有代表性的一段光滑曲线，以此代表汽车的工况，这是重型汽车行驶状况构建的基本思想.

模型建立 以上述思想为理念，可以构建具体的模型.首先将重型汽车长时间的行驶划分为各个片段，再用聚类分析的方法将各片段分为多个大类，片段会随时间变化在各个大类中转移，可以将这个过程看成马尔科夫链.设定片段拼接优化指标D，根据马尔科夫链的性质和指标D确定片段重组标准，最后检验重组的片段是否能代表实际行驶状况.模型步骤框图如图2.7所示.

图2.7 模型步骤框图

模型求解 重型汽车在公路上行驶是不规律的，建立片段切割标准将重型汽

车行驶时间进行分段切割.常用的切割方法有速度变化切割、最大似然估计和突变点切割.

速度变化切割是根据速度的变化情况进行切割.根据时刻-速度图像，以及加速时段、减速时段和匀速时段这三大类情况进行切割.例如，重型汽车加速到顶点的时段记为片段1，随后匀速行驶一段时间记为片段2，之后一直减速的时段记为片段3，如此切割下去，如图2.8所示.

图2.8 速度变化切割示意图

还可以根据有序聚类分析法找出序列的突变点，根据公式

以此来进行切割，示例如图2.9所示。

选择一种切割方式可以得到多个片段，随后采用系统聚类法对片段进行分类，样例采用计算类间距离中的离差平方和法（Ward法）进行分类，得到图2.10所示的树状图.

图2.9 根据序列的突变点切割示意图

图2.10 采用离差平方和法分类得到的树状图

根据需要得到分类线，图中的分类线将所有片段分为5大类，根据片段在大类之间转移的频数可以写出状态转移概率矩阵.

由于片段间转移的概率是不会相互影响的，从片段开始到结束，然后转移到新的片段，其转移过程可以看作马尔科夫链{Xt,t≥1}.片段一共分5大类，所以该马尔科夫Xt的状态空间为E={1,2,3,4,5}.

由于无法知道概率，所以用频数代替概率来写出转移概率矩阵，因为数据量很大，这样的假定通常是合理的.设大类i到大类j的转移概率为

其中nij为大类i转移到大类j的频数.根据计算公式和数据，可以得到相应的一步转移概率矩阵为

很明显，转移概率矩阵行的和为1.同时可以看出，大类1停留在大类1内的概率超过0.5，大类1转移到大类2的概率大约在0.2，大类1转移到大类3、4、5的概率都很小，但不为0；大类2转移到大类1的概率为0.6，大类2停留到大类2的概率为0.2，大类2转移到大类3的概率为0.2，大类2不会转移到大类4或大类5；其余大类的转移情况可以根据一步转移概率矩阵得到.

接下来对片段进行重组，要求重组的片段能够代表全部的运行状态，所以需要制定标准来选择片段进行重组.速度和加速度是重型汽车运行过程中两个重要的指标，可以根据速度和加速度的实际情况和重组情况的差别来制定标准.

首先根据速度和加速度可以写出联合概率分布.试验数据出现的加速度有1,0,-1,-2,-3,-4,-5；试验数据出现的速度在40～90之间，所以可以将速度分成[40,50）,[50,60）,[60,70）,[70,80）,[80,90）5个区间.用相应区间内的频数除以总频数即得到相应的频率分布，近似作为速度的概率分布，经过计算得到的速度-加速度联合概率分布如表2.6所示.

表2.6 速度-加速度联合概率分布

同理可以得到重组片段的速度-加速度联合概率分布.根据实际的速度-加速度联合概率分布和重组片段的速度-加速度联合概率分布，可以得到速度-加速度联合概率分布偏差

其中Pij为试验数据的速度-加速度分布概率，Qij为提取片段后建立的工况的速度和加速度分布概率.

片段重组的基本思想是局部最优法，要求每拼接一个新的片段，该新片段的D值相较其他未拼接片段最小，且前后拼接处的速度、加速度变化符合实际发动机工作原理，这样的曲线才比较符合实际情况，即片段与片段应该是基本连续的，不能出现大的跳跃；最后要求片段到片段是可以转移的，即相应转移概率要大于0.以此制订出片段重组的4个步骤：

（1）计算70个行驶片段的速度-加速度概率分布，计算出D值，选择D值最小的片段作为初始片段；

（2）选择下一个片段的起始速度与前个一片段的末速度差距必须保持在可达范围之内，即试验数据中允许的加速度范围之内；

（3）前一片段与后一片段所属的状态转移概率要大于0，能够转移才能进行片段重组；

（4）每选择一次片段，对于新合成的工况，都要使D值相对最小.

按上述标准得到重组后的简化工况，和简化前的工况进行对比，得到的结果如图2.11所示.

图2.11 简化工况前后的对比

为了做对比，将原实际工况约两周的时间段压缩为一小时左右，和简化后工况同长.可以看出，简化后的工况相比简化前的工况更加平滑，也非常接近于简化前的工况.用简化后的工况代表汽车在公路上的行驶状态，从而为研究如何降低重型汽车耗油量提供帮助.

0 Z: ?/ c' `3 v, H0 U0 f


% D( T! [8 t) f4 C
$ \- X8 ?8 s* x4 k