2.4　马尔可夫链法建模

马尔可夫链模型在自然科学、工程技术、社会科学、经济研究等领域都有广泛的应用。它可以解决一些随机系统的随机转移和确定系统的状态转移问题。马尔可夫链预测法是应用概率论中马尔可夫链（Markov Chain）的理论和方法来研究分析时间序列的变化规律，并由此预测其未来变化趋势的一种动态预测技术。这种技术已在市场预测分析、市场管理决策、卫生事业管理、卫生经济研究、生物遗传等方面得到了广泛的应用。本节主要介绍马尔可夫链的基本原理以及运用原理去解决一些简单问题的基本方法，最后通过实例演示马尔可夫链建模的基本思想及建模思路。

2.4.1　马尔可夫链的基本知识

（1）随机过程的概念

定义2.1　设集合{Xt：t∈T}是一族随机变量，T是一个实数集合，如果对于任意t∈T，Xt是一个随机变量，则称{Xt：t∈T}是一个随机过程。

t为参数，一般情况下可以认为是时间，称t的取值集合T为参数集合；随机变量Xt的取得每一个可能值，称为随机过程的一个状态。其全体可能值构成的集合，称为随机过程的状态空间，用E表示；当参数集合T为非负整数集时，随机过程又称为随机序列，随机序列可用{Xn:n=1，2，3，…}表示。当T为时间时，该随机序列就是时间序列。马尔可夫链，也称为马氏链，就是一种特殊的随机时间序列。

（2）马尔可夫链。

定义2.2　设{Xn:n=1，2，3，…}是一个随机序列，状态空间E为有限或可列集。若对于任意正整数m，n。如果i∈E，j∈E，ik∈E（k=1，2，…，n-1）满足P（Xn+m=j｜Xn=i，Xn-1=in-1，…，X1=i1）=P（Xn+m=j｜Xn=i）成立，则称随机序列{Xn:n=1，2，3，…}为一个马尔可夫链，简称为马氏链。

从该定义可知，随机变量Xn=i，是指第n步时，随机变量Xn处于状态i。条件概率P（Xn+m=j｜Xn=i）是指第n步时，随机变量Xn处于状态i的条件下，第n+m步随机变量Xn+m处于状态j的条件概率。

两个条件概率相等，说明第n+m步时的随机变量Xn+m所处的状态j，只与第n步时的随机变量Xn所处的状态i有关，而与第n步前面的所有n-1步的随机变量们所处的状态无关，将此称为随机变量序列的无后效性。具有无后效性的随机变量序列称为马尔可夫链，即为马氏链。

（3）时齐的马尔可夫链

定义2.3　如果马尔可夫链{Xn:n=1，2，3，…}中的条件概率P（Xn+m=j｜Xn=i）与n无关，则称此马尔可夫链{Xn:n=1，2，3，…}是时齐的马尔可夫链。

时齐的马尔可夫链{Xn:n=1，2，3，…}中的条件概率P（Xn+m=j｜Xn=i）与n无关，只与m有关，即无论从第几步的状态i出发，再经过m步，到达状态j的概率都相等，所以可将该条件概率记为P（Xn+m=j｜Xn=i）=pij（m），即系统由状态i出发，经过m步（m个时间段）的转移到达系统状态j的概率是pij（m）。

概率pij（m）一般称为转移概率。由时齐性的定义可知，系统由状态i转移到状态j的转移概率pij（m），只依赖于时间间隔m（步数）的长短，与起始的时刻无关。

（4）转移概率矩阵

定义2.4　对于一个马尔可夫链{Xn:n=1，2，3，…}，称以m步转移概率pij（m）为元素的矩阵为转移概率矩阵，记为P（m）。

m=1时的转移概率矩阵称为一步转移概率矩阵，记为P=P（1）=（pij），其中pij=pij（1），简称为转移矩阵，记为

马尔可夫链{Xn:n=1，2，3，…}的转移概率矩阵P（m）=（pij（m））的性质如下。

性质2.1　对于一切i∈E，j∈E，有0≤pij（m）≤1成立。

性质2.2　对于一切i∈E，有成立。

性质2.3　对于一切i∈E，j∈E，使得成立。

转移概率矩阵P（m）既可以是有限矩阵，也可以是无限矩阵，本节主要考虑状态集E为有限集的马尔可夫链，此时对应的转移概率矩阵P和P（m）都是有限矩阵，并且有P（m+1）=P（m）P，P（m）=Pm成立。

由此可见，对具有N个状态的马尔可夫链而言，描述它的概率性质，最重要的是它在n时刻处于状态i下一时刻转移到状态j的一步转移概率：P（Xn+1=j｜Xn=i）=pij。

【例2.3】　设某抗病毒药销售情况分为“畅销”和“滞销”两种，以“1”代表“畅销”，“2”代表“滞销”。以Xn表示第n个季度的销售状态，则Xn可以取值1或2。若未来的抗病毒药销售状态，只与现在的市场状态有关，而与以前的市场状态无关，则抗病毒药的市场状态{Xn，n≥1}就构成一个马尔可夫链。设一步转移概率矩阵为

则其二步转移矩阵为

同理可以求出其k步转移概率矩阵P（k）=Pk。

（5）状态概率

系统的状态用随机变量Xn表示，ai（n）=P（Xn=i），i=1，2，…，N称为状态概率。

m步时系统所处的状态概率向量为a（m）=（a1（m），a2（m），…，aN（m）），利用马尔可夫链的一步转移概率矩阵P，可以计算出关系式a（m+1）=a（m）P成立，且有m步时系统所处的状态概率向量a（m）=（a1（m），a2（m），…，aN（m））=（a1（0），a2（0），…，aN（0））Pm，记a（0）为初始状态概率向量，上式又可以记为a（m）=a（0）Pm。

（6）马尔可夫链的正则链

定义2.5　设{Xn:n=1，2，3，…，N}为马尔可夫链，如果存在一个正整数k，使得马尔可夫链从任意状态i出发，经过k步转移都以大于零的概率到达状态j，则称此时的马尔可夫链{Xn:n=1，2，3，…，N}为正则链。

定理2.1　若马尔可夫链的转移矩阵为P，则它是正则链的充分必要条件是，存在正整数k，使得Pk>0（矩阵的每一个元素大于零）。

（7）马尔可夫链的遍历性

定义2.6　如果马尔可夫链{Xn:n=1，2，3，…，N}为正则链，则该马尔可夫链{Xn:n=1，2，3，…，N}具有遍历性（或称是平稳的）。

即存在一个向量（有限维向量）W=（p1，p2，…，pN）使得，与初始状态概率向量a（0）无关，W称为稳态概率或极限状态分布。

且有，即向量W=（p1，p2，…，pN）是它的稳定解。

若当转移步数n→+∞时，系统所处的状态记为（a1（+∞），a2（+∞），…，ak（+∞）），则有（a1（+∞），a2（+∞），…，ak（+∞））=W。同时由上面的结果又可以推出，。其中P（m）=（pij（m））是m步转移概率矩阵，P=P（1）=（pij）是1步转移概率矩阵。

（8）状态转移概率的估算

在马尔可夫链建模方法中，系统状态的转移概率的估算非常重要。估算的方法通常有两种。一是主观概率法。它是根据人们长期积累的经验以及对预测事件的了解，对事件发生的可能性大小的一种主观估计，这种方法一般是在缺乏历史统计资料或资料不全的情况下使用。二是统计估算法。主观概率法是一种需要根据以往积累的经验，提前给定概率分布，然后再利用马氏链进行分析的方法。而在已知历史统计资料的情况下，可以使用统计估算法。统计估算法估计转移概率的方法如下：

假定系统有N种状态S1，S2，…，SN，根据系统的状态转移的历史记录，得到表2.7的统计表格，以表示系统从状态i转移到状态j的转移概率估计值，则由表2.7的数据计算估计值的公式如下：

　　 （2.62） 　　

【例2.4】　设某系统有3种状态S1，S2和S3，系统状态的转移情况见表2.8。试求系统的状态转移概率矩阵。

【解】　由式（2.62），得

故系统的转移概率矩阵为。

（9）马尔可夫链的应用步骤

第一步：根据实际问题设定随机变量Xn，建立随机变量序列，并判断随机变量序列是否为时齐的马尔可夫链。

第二步：建立第n（步）个时段与第n+1（步）个时段间的系统状态之间的联系，并用转移概率pij表示出来，构造出转移概率矩阵P。在系统具有N个状态的情况下，对系统在时刻k的状态进行预测。

第三步：利用转移概率矩阵P=（pij），寻求是否存在正整数k，使得P（k）=Pk>0，判断系统是否为正则链，并求其平稳分布（即极限状态分布）。

【例2.5】　市场占有率预测　预测A、B、C三个厂家生产的某种抗病毒药在未来的市场占有情况。

【解】　首先进行市场调查。主要调查以下两件事。

①目前的市场占有情况.如在购买该药的总共1000家对象（购买力相当的医院、药店等）中，买A、B、C三药厂的各有400家、300家、300家，那么A、B、C三药厂目前的市场占有份额分别为：40％、30％、30％。则S（0）=（0.4，0.3，0.3）为目前市场的占有分布或称初始概率分布。

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

其次，根据上述马尔可夫链的应用步骤，计算转移概率矩阵，建立数学模型，并利用k步的转移概率矩阵进行预测。

假定在未来的时期内，顾客相同间隔时间的流动情况不因时期的不同而发生变化，以1、2、3分别表示顾客买A、B、C三厂家的药这三个状态，以季度为模型的步长（即转移一步所需的时间），那么根据表2.9，我们可以得模型的转移概率矩阵：

矩阵中的第1行（0.4，0.3，0.3）表示目前是A厂的顾客下季度有40％仍买A厂的药，转为买B厂和C厂的各有30％.同样，第2行、第3行分别表示目前是B厂和C厂的顾客下季度的流向。

由P我们可以计算任意的k步转移矩阵，如三步转移矩阵：

从这个矩阵的各行可知三个季度以后各厂家顾客的流动情况.如从第二行（0.504，0.252，0.244）知，B厂的顾客三个季度后有50.4％转向买A厂的药，25.2％仍买B厂的，24.4％转向买C厂的药。

设表示预测对象k季度以后的市场占有率，初始分布则为，市场占有率的预测模型为

S（k）=S（0）Pk=S（k-1）P

现在，由第一步，计算有S（0）=（0.4，0.3，0.3），由此，可预测任意时期A、B、C三厂家的市场占有率。例如，三个季度以后的预测值为：

大致上，A厂占有一半的市场，B厂、C厂各占四分之一。

最后，判别马尔可夫链的正则性，对最终市场的占有率进行预测。

由于一步转移概率矩阵P>0，故马尔可夫链为正则链，稳定分布存在。当市场出现平衡状态时，令稳定分布为S=（p1，p2，p3），由方程SP=S，即得

再加上条件p1+p2+p3=1，经整理得

该方程组是三个变量四个方程的方程组，在前三个方程中只有二个是独立的，任意删去一个，从剩下的三个方程中，可求出唯一解：

p1=0.5，p2=0.25，p3=0.25

即为A、B、C三家的最终市场占有率。

2.4.2　有利润的马尔可夫链

在马尔可夫链模型中，随着时间的推移，系统的状态可能发生转移，这种转移常常会引起某种经济指标的变化。这种随着系统的状态转移，赋予一定利润的马尔可夫链，称为有利润的马尔可夫链。对于一般的具有转移矩阵

的马尔可夫链，当系统由i转移到j时，赋予利润rij（i，j=1，2，…，N），则称

为系统的利润矩阵，rij>0称为盈利，rij<0称为亏本，rij=0称为不亏不盈。

随着时间的变化，系统的状态不断地转移，从而可得到一系列利润，由于状态的转移是随机的，因而一系列的利润是随机变量，其概率关系由马氏链的转移概率决定。

一般地定义k步转移利润随机变量的分布为：

则系统处于状态i经过k步转移后所得的期望利润的递推计算式为：

当k=1时，规定边界条件。

称一步转移的期望利润为即时的期望利润，并记

【例2.6】　期望利润预测　企业追逐市场占有率的真正目的是使利润增加，因此竞争各方无论是为了夺回市场份额，还是为了保住或者提高市场份额，在制订对策时都必须对期望利润进行预测。

预测主要分以下两步进行。①市场统计调查。首先调查销路的变化情况，即查清由畅销到滞销或由滞销到畅销，连续畅销或连续滞销的可能性是多少。其次统计出由于销路的变化，获得的利润和亏损情况。②建立数学模型，列出预测公式进行预测。

通过市场调查，得到如下的销路转移表（表2.10）和利润变化表（表2.11）。由此，我们来建立数学模型。

销路转移表说明连续畅销的可能性为50％，由畅销转入滞销的可能性也是50％，由滞销到畅销为40％，连续滞销的可能性为60％。利润表说明的是连续畅销获利900万元，由畅销到滞销或由滞销到畅销均获利300万元，连续滞销则亏损700万元。从而得到销售状态的转移矩阵P和利润矩阵R如下：

P和R便构成一个有利润的马尔可夫链。由前面所述的基本原理及公式得即时期利润的预测公式：

k步以后的期望利润：将调查数据代入则可预测各时期的期望利润值。如：

q1=9×0.5+3×0.5=6

q2=3×0.4-7×0.6=-3

由此可知，当本季度处于畅销时，在下一季度可以期望获得利润600万元；当本季度处于滞销时，下一季度将期望亏损300万元。

　　 同样算得： 　　

由此可预测本季度处于畅销时，两个季度后可期望获利750万元，三个季度后可期望获利855万元；当本季度处于滞销时，两个季度后将亏损240万元，三个季度后亏损144万元。

2.4.3　案例分析

（1）马尔可夫模型在流行病监测中的应用

马尔可夫模型是用于描述时间和状态都是离散的随机过程的数学模型。应用其理论和方法，可以对疾病发病情况随时间序列的变化规律进行分析和研究，预测疾病的发展变化趋势，为预防和控制疾病提供依据。表2.12统计了某市1980～1995年肾综合征出血热（HFRS）的发病率分别为（单位：1/10万）：2.95、6.28、10.28、7.01、7.36、13.78、33.93、35.87、33.40、28.38、30.50、33.79、39.70、30.39、39.70、33.59，表2.13给出了状态取值和初始概率分布。试建立模型对疾病的未来发展趋势进行预测。

首先根据资料将发病率划分为四个状态，统计各数据的状态归属及各状态出现的频率（初始概率），得表2.12和表2.13。

由表2.12可得各状态的转移频率即状态转移概率的估计值，从而得模型的一步转移概率矩阵：

可认为HFRS下一年的发病率只与当年发病率有关，而与过去的发病率无关，且任意时期的一步转移概率矩阵不变，从而满足无后效性和平稳性的假设，因而可用初始分布为（4/16，2/16，1/16，9/16），转移概率矩阵为P的马氏链模型来预测HFRS发病率未来的情况。

计算多步转移矩阵：

计算极限或解方程（p1，p2，p3，p4）=（p1，p2，p3，p4）P，，得模型的极限概率分布（稳态分布）：（0，0，1/9，8/9）。

分析预测：由于1995年处于状态4，比较P的第4行的四个数字知，p44=0.875最大，所以预测1996年仍处于状态4，即发病率大于30/10万。同样，从二、三、四步转移矩阵知，依然是状态4转入状态4的概率最大，所以预测1996～1999年该市的HFRS发病率将持续在大于30/10万（高发区）水平，这提醒我们应该对此高度重视，采取相应对策。

如果转移概率矩阵始终不变，从极限分布看，最终HFRS发病率将保持在高发区水平，当然，这应该是不会符合实际情况的，因为随着各方面因素的改变，转移概率矩阵一般也会发生变化。所以马尔可夫链模型主要适用于短期预测。

在用马尔可夫链模型进行预测的过程中，无后效性和平稳性是最基本的要求，而模型是否合理有效，状态的划分和转移概率矩阵的估算是关键，不同的状态划分可能会得到不同的结果，通常我们根据有关预测对象的专业知识和数据的多少及范围来确定系统状态。

在卫生管理事业中，用马尔可夫链模型还可预测医疗器械、药品的市场占有率，药品的期望利润收益等。

（2）最佳进货策略

一个小型的水族馆专营各种规格的水族箱，每个周末，店老板都要清点存货，确定下一周是否进货。老板的进货策略是：如果本周某种规格的水族箱存货全部售出的话，下周初就再进货3个，否则便不再进货。这样的策略可能会造成部分时间顾客买不到货，造成一定的潜在利润损失。表2.14给出了该店过去两年的需求情况，根据这组统计数据，确定该店的缺货情况。

顾客的购买是随机行为，通过数据分析发现，该水族馆平均每周需求是1个水族箱，需求近似服从参数为1的泊松分布。

表2.14的数据的直方图见图2.16。

模型假设如下。

（ⅰ）第n周水族馆的需求Qn服从参数为1的泊松分布。

（ⅱ）第n周水族馆的存货用随机变量Xn表示，且X0=3。

模型分析及求解如下。

根据模型假设得到：

于是

P（Xn+1=1｜Xn=1）=P（Qn=0）=0.3679，P（Xn+1=2｜Xn=1）=0

P（Xn+1=3｜Xn=1）=1-P（Qn=1）=1-0.3679=0.6321

P（Xn+1=1｜Xn=2）=P（Qn=1）=0.3679

P（Xn+1=2｜Xn=2）=P（Qn=0）=0.3679

P（Xn+1=3｜Xn=2）=P（Qn=2）+P（Qn≥3）=0.2642

P（Xn+1=1｜Xn=3）=P（Qn=2）=0.1840

P（Xn+1=2｜Xn=3）=P（Qn=1）=0.3679

P（Xn+1=3｜Xn=3）=P（Qn=0）+P（Qn≥3）=0.4481

如果记pij=（Xn+1=j｜Xn=i），则有

下面进一步计算缺货概率P（Qn>Xn）。一般来说这个概率依赖于n，为了得到关于缺货的更一般的信息，我们需要对Xn的信息再作一些分析。

{Xn}是一个遍历的马尔可夫链，其一定存在唯一的渐进稳定的单位概率向量π，它可以通过求解稳定状态方程计算出来。令π=πP得到

π=（0.2847，0.2631，0.4521）

因此，对于充分大的n，近似地有

P（Xn=1）=0.2847，P（Xn=2）=0.2631，P（Xn=3）=0.4521，

上式表明，每年约有10.49％的时间（约5周时间）水族馆是缺货的。

（3）定岗定编问题

社会系统中，常常因为职务、地位等的不同，划分出许多的等级，各等级的人数比例称为等级结构，定岗定编问题即是保持一个稳定合理的等级机构，这类问题在许多单位都可以看到它的缩影，大学教师的各种职称、工程技术人员的各种级别等。那么等级结构是怎样随时间变化的呢？

等级结构的变化依赖系统内部的等级随时间的转移（即通常所说的职务升降），以及系统内外部的交流（即通常所说的调入、调出、退休、死亡等）。通过数学语言将等级结构随时间变化关系恰当地表示出来，就构成这个问题的数学模型。

①模型假设及符号约定如下。

ⅰ.将一个系统由低向高分成m个等级，每隔s年进行一次正常的等级调整。

ⅱ.ni（k）表示第k次调整时第i个等级的人数，记n（k）=［n1（k），n2（k），…，nm（k）］，不妨称之为等级结构为系统第k年的总人数。

ⅲ.记ai（k）=ni（k）/N（k），a（k）=［a1（k），a2（k）…，am（k）］，a（k）称为等级结构向量。

ⅳ.记，表示每次从等级i升到等级j的人数占等级i中人数的比例。

ⅴ.w=［w1，w2，…，wm］，wi表示每次从等级i中退出人数的比例。

ⅵ.r=［r1，r2，…，rm］，ri表示每次调入等级i的人数占总调入人数的比例。

ⅶ.记R（k）为第k次调入总人数。

ⅷ.记W（k）为第k次退出总人数。

ⅸ.从k到k+1年总人数的增长量记为M（k）。

一般地，P0，w，r分别称为内部转移矩阵、退出向量、调入向量。为简便起见，不妨假设其均与时间无关。

②模型建立及求解。根据假设可知，

同时又可以得到，wi≥0，ri≥0，且，。

则第k+1次的系统总人数满足方程

N（k+1）=N（k）+R（k）-W（k）

每个等级人数的转移方程为

n（k+1）=n（k）P0+R（k）r　　（2.63）

从k到k+1年总人数的增长量M（k）有

R（k）=W（k）+M（k）=n（k）wT+M（k）　　（2.64）

将式（2.64）代入式（2.63）得到

n（k+1）=n（k）（P0+wTr）+M（k）r　　（2.65）

记P=P0+wTr，则P也是随机矩阵，式（2.65）可以表示为

n（k+1）=n（k）P+M（k）r

通常称该式为等级分布基本方程。

假如系统的总人数每年以固定的比例增长，即M（k）=βN（k），则

a（k+1）=（1+β）-1［a（k）P+βr］

特别地，如果每年进出系统的人数大致相等，即系统总人数N（k）保持不变。那么，M（k）=0，β=0，上述方程可以简化为a（k+1）=a（k）P，即为一个马尔可夫链。

对于给出的等级分布基本方程，下面考虑如下问题：给定初始等级结构a（0），如何确定调入比例r*，使等级变化尽快达到或接近给定的理想等级结构a*。需要指出的是，如果等级结构满足，则称等级结构为稳定的。

系统是否有稳定的等级结构是有条件的，如果存在，则r必须满足

，，且有

　　 （2.66） 　　

保证式（2.66）成立的充分必要条件是存在非负向量a（每个分量非负），使

a≥aP0

如果矩阵Q=E-P0可逆，由式（2.66）得到

　　 （2.67） 　　

令e=（1，1，…，1），由于的各分量之和为1，即。利用式（2.67）得

　　 （2.68） 　　

再将式（2.68）代入式（2.67）得到

在处理实际问题时，通常会比较两个等级的接近程度，以便确定当前等级的状态。为此，我们引入等级距离的概念。定义两个等级a（1），a（2）之间的距离D（a（1），a（2））如下：

其中λi为加权因子，由对各等级的关注程度确定。一个满意的等级分布应该满足如下优化问题：

　　 （2.69） 　　

由于

如果记，则与Y-r成正比，式（2.69）等价于

对于上面的优化问题，在某种程度上，它只是条件极值问题，可以用拉格朗日乘子法求解。

这个模型不仅可以用来描述社会中的等级结构，也可以研究不同部门人员的流动，甚至是不同产业中的劳动力构成，应用范围非常广泛。本模型只是讨论的最简单、最基本的情况，在增加条件和控制变量的情况下，可以给出更灵活、更复杂的模型。