Markov链在很多领域中被用来建立数学模型,该模型常用来描述用同一种方法进行多次实验或测量,且每次实验结果属于几个指定的可能结果之一,每次实验的结果仅仅依赖于最近的上一次实验结果。
例:实验室动物每天可以吃三种食物中的任一种。结果表明,如果第一次测试中,这个动物选择一种食物,则下一次测试中它选同样的食物的概率是50%,下一次测试它选择其他两种实物的概率都是25%。
(1)求该实验的随机矩阵。
(2)如果这个动物第一次测试选了1号食物,在第二次测试中选择2号食物的概率、第三次选择食物2和3的概率是多少?
解:
(1)随机矩阵PPP:
(2)在MATLAB中求P2=P×PP^2=P\times PP2=P×P:
P =
0.5000 0.2500 0.2500
0.2500 0.5000 0.2500
0.2500 0.2500 0.5000
>> P^2
ans =
0.3750 0.3125 0.3125
0.3125 0.3750 0.3125
0.3125 0.3125 0.3750
即P2P^2P2为:
由矩阵PPP可求出该动物第二次测试选择各种食物的概率表:v1=Pv0=[0.50000.25000.2500]v_1=Pv_0=\begin{bmatrix}0.5000\\0.2500\\0.2500\end{bmatrix}v1=Pv0=⎣⎡0.50000.25000.2500⎦⎤,所以第二次选择2号食物概率为v1第二行=0.2500v_{1第二行}=0.2500v1第二行=0.2500。
由矩阵P2P^2P2可以求出该动物第三次测试选择各种食物的概率,设第一次选了第一种食物的初始概率向量v0=[100]v_0=\begin{bmatrix}1\\0\\0\end{bmatrix}v0=⎣⎡100⎦⎤,则第三次的概率向量为:v2=P2v0=[0.37500.31250.3125]v_2=P^2v_0=\begin{bmatrix}0.3750\\0.3125\\0.3125\end{bmatrix}v2=P2v0=⎣⎡0.37500.31250.3125⎦⎤,所以第三次选择食物2和3概率分别为v2第二行=0.3125v_{2第二行}=0.3125v2第二行=0.3125和v2第三行=0.3125v_{2第三行}=0.3125v2第三行=0.3125。
概率向量,随机矩阵(Stochastic Matrix),马尔科夫链
一个具有非负元素且各元素相加为1的向量为概率向量。由概率向量作为列向量组合而成的方阵为随机矩阵(也叫马尔科夫矩阵、转移矩阵等)。
例如,上例中,矩阵PPP为随机矩阵,PPP的每一列为概率向量。
马尔科夫链是一个向量序列x0,x1,x2…x_0,x_1,x_2\dotsx0,x1,x2…与一个随机矩阵PPP的关系,满足:
x1=Px0, x2=Px1, x3=Px2, …x_1=Px_0,\space x_2=Px_1,\space x_3=Px_2,\space \dotsx1=Px0, x2=Px1, x3=Px2, …
此关系可以用一阶差分方程来描述:
xk+1=Pxk, k=0,1,2,…x_{k+1}=Px_k, \space k=0,1,2,\dotsxk+1=Pxk, k=0,1,2,…
当向量在RnR^nRn中的一个马尔科夫链描述一个系统或实验的序列时,xkx_kxk中的元素分别列出了系统在nnn个可能的状态中的概率,或者实验结果可能是nnn个可能的结果之一的概率。所以xkx_kxk也称为状态向量。
稳态向量(Steady-State Vector)
如果PPP是一个随机矩阵,如果一个状态向量qqq满足Pq=qPq=qPq=q,则此向量qqq为相对于PPP的稳态向量(也可叫平衡向量)。
每个随机矩阵都有一个稳态向量。
例:求上例的稳态向量xxx。
根据稳态向量的定义,满足Px=xPx=xPx=x,即(P−I)x=0(P-I)x=0(P−I)x=0
用上面的随机矩阵PPP构建齐次方程:(P−I)x=0(P-I)x=0(P−I)x=0
再解方程:
>> I=eye(3) //生成3阶单位矩阵I
I =
1 0 0
0 1 0
0 0 1
>> rref(P-I) //将矩阵(P-I)化简为简化阶梯阵
ans =
1 0 -1
0 1 -1
0 0 0
对齐次方程组,可以直接对系数矩阵进行消元称为标准阶梯阵,效果和对增广矩阵[P−IO][P-I\quad O][P−IO]消元是一样的,下面是对增广矩阵的消元法以解线性方程组:
>> O = [0;0;0]
O =
0
0
0
>> A = [P-I O]
A =
-0.5000 0.2500 0.2500 0
0.2500 -0.5000 0.2500 0
0.2500 0.2500 -0.5000 0
>> rref(A)
ans =
1 0 -1 0
0 1 -1 0
0 0 0 0
可以发现,消元出来的效果,和直接用系数矩阵来消元是一样的。
通解为:
x=[x1x2x3]=[x3x3x3]=x3[111]x=\begin{bmatrix}x_1\\x_2\\x_3\end{bmatrix}=\begin{bmatrix}x_3\\x_3\\x_3\end{bmatrix}=x_3\begin{bmatrix}1\\1\\1\end{bmatrix}x=⎣⎡x1x2x3⎦⎤=⎣⎡x3x3x3⎦⎤=x3⎣⎡111⎦⎤
令自由变量x3=1x_3=1x3=1,则此特殊解为:w=[x1x2x3]=[111]w=\begin{bmatrix}x_1\\x_2\\x_3\end{bmatrix}=\begin{bmatrix}1\\1\\1\end{bmatrix}w=⎣⎡x1x2x3⎦⎤=⎣⎡111⎦⎤;
也可令x3=2x_3=2x3=2,则此特殊解为:u=[x1x2x3]=[222]u=\begin{bmatrix}x_1\\x_2\\x_3\end{bmatrix}=\begin{bmatrix}2\\2\\2\end{bmatrix}u=⎣⎡x1x2x3⎦⎤=⎣⎡222⎦⎤。
则求P×[111]P\times \begin{bmatrix}1\\1\\1\end{bmatrix}P×⎣⎡111⎦⎤和P×[222]P\times \begin{bmatrix}2\\2\\2\end{bmatrix}P×⎣⎡222⎦⎤:
>> w = [1;1;1]
w =
1
1
1
>> product = P*w
product =
1
1
1
>> u = [2;2;2]
u =
2
2
2
>> product2 = P*u
product2 =
2
2
2
>> q = w/sum(w) //注意,这里的q向量是概率向量
q =
0.3333
0.3333
0.3333
>> product3 = P*q
product3 =
0.3333
0.3333
0.3333
可见,P×w=wP\times w=wP×w=w,P×u=uP\times u=uP×u=u和P×q=qP\times q=qP×q=q。
但是根据定义,向量qqq是Px=xPx=xPx=x,即(P−I)x=0(P-I)x=0(P−I)x=0的全体解的结合的概率向量(求法:使用任意一个解,如上面的www,除以此解www每个元素之和即得到概率向量qqq,可验算qqq的各元素之和为1,且满足P×q=qP\times q=qP×q=q),所以向量qqq随机矩阵PPP的一个稳态向量。上面的www和uuu因为各元素和不为1,所以不能叫概率向量,所以就不能是稳态向量。
正则随机矩阵(Regular Stochastic Matrix):
如果一个随机矩阵的某次幂PkP^kPk(注意,是有某次幂就行)的每个元素仅包含严格正的元素,则此随机矩阵可称为正则随机矩阵。
对于上面的例子的随机矩阵PPP和P2P^2P2,每项元素都是正数,所以可以说PPP是正则随机矩阵。
例:判断下列矩阵是否为正则随机矩阵:
P1=[0.70.20.600.20.40.30.60]P_1=\begin{bmatrix}0.7&0.2&0.6\\0&0.2&0.4\\0.3&0.6&0\end{bmatrix}P1=⎣⎡0.700.30.20.20.60.60.40⎦⎤
P2=[0.70.20.600.200.30.60.4]P_2=\begin{bmatrix}0.7&0.2&0.6\\0&0.2&0\\0.3&0.6&0.4\end{bmatrix}P2=⎣⎡0.700.30.20.20.60.600.4⎦⎤
P3=[010100001]P_3=\begin{bmatrix}0&1&0\\1&0&0\\0&0&1\end{bmatrix}P3=⎣⎡010100001⎦⎤
解:
MATLAB计算
P1 =
0.7000 0.2000 0.6000
0 0.2000 0.4000
0.3000 0.6000 0
>> (P1)^2
ans =
0.6700 0.5400 0.5000
0.1200 0.2800 0.0800
0.2100 0.1800 0.4200
可见,找到了一个P12P_1^2P12是全部严格正的(没有负数和零元素),且P1P_1P1每列和为1,所以P1P_1P1是正则随机矩阵。
P2 =
0.7000 0.2000 0.6000
0 0.2000 0
0.3000 0.6000 0.4000
>> (P2)^2
ans =
0.6700 0.5400 0.6600
0 0.0400 0
0.3300 0.4200 0.3400
>> (P2)^3
ans =
0.6670 0.6380 0.6660
0 0.0080 0
0.3330 0.3540 0.3340
>> (P2)^4
ans =
0.6667 0.6606 0.6666
0 0.0016 0
0.3333 0.3378 0.3334
>> (P2)^5
ans =
0.6667 0.6654 0.6667
0 0.0003 0
0.3333 0.3343 0.3333
...
可见,P2P_2P2矩阵无论多少次方,都含有0,虽然每列和都为1(所以是随机矩阵),所以不是正则随机矩阵的。
P3 =
0 1 0
1 0 0
0 0 1
>> (P3)^2
ans =
1 0 0
0 1 0
0 0 1
>> (P3)^3
ans =
0 1 0
1 0 0
0 0 1
>> (P3)^4
ans =
1 0 0
0 1 0
0 0 1
>> (P3)^5
ans =
0 1 0
1 0 0
0 0 1
>> (P3)^6
ans =
1 0 0
0 1 0
0 0 1
矩阵P3P_3P3也因为存在0,所以不是正则随机矩阵。
对此三个矩阵求稳态向量:
>> P1O=[P1-I O]
P1O =
-0.3000 0.2000 0.6000 0
0 -0.8000 0.4000 0
0.3000 0.6000 -1.0000 0
>> rref(P1O)
ans =
1.0000 0 -2.3333 0
0 1.0000 -0.5000 0
0 0 0 0
>> P2O=[P2-I O]
P2O =
-0.3000 0.2000 0.6000 0
0 -0.8000 0 0
0.3000 0.6000 -0.6000 0
>> rref(P2O)
ans =
1.0000 0 -2.0000 0
0 1.0000 0 0
0 0 0 0
>> P3O=[P3-I O]
P3O =
-1 1 0 0
1 -1 0 0
0 0 0 0
>> rref(P3O)
ans =
1 -1 0 0
0 0 0 0
0 0 0 0
对于P1P_1P1,通解:x=[x1x2x3]=x3[2.30.51]x=\begin{bmatrix}x_1\\x_2\\x_3\end{bmatrix}=x_3\begin{bmatrix}2.3\\0.5\\1\end{bmatrix}x=⎣⎡x1x2x3⎦⎤=x3⎣⎡2.30.51⎦⎤,令x3=1x_3=1x3=1,q=x/sum(x)=[0.60530.13160.2632]q=x/sum(x)=\begin{bmatrix}0.6053\\0.1316\\0.2632\end{bmatrix}q=x/sum(x)=⎣⎡0.60530.13160.2632⎦⎤。
对于P2P_2P2,通解:x=[x1x2x3]=[000]x=\begin{bmatrix}x_1\\x_2\\x_3\end{bmatrix}=\begin{bmatrix}0\\0\\0\end{bmatrix}x=⎣⎡x1x2x3⎦⎤=⎣⎡000⎦⎤,q=x/sum(x)q=x/sum(x)q=x/sum(x)不存在。
对于P3P_3P3,通解:x=[x1x2x3]=x2[110]+x3[001]x=\begin{bmatrix}x_1\\x_2\\x_3\end{bmatrix}=x_2\begin{bmatrix}1\\1\\0\end{bmatrix}+x_3\begin{bmatrix}0\\0\\1\end{bmatrix}x=⎣⎡x1x2x3⎦⎤=x2⎣⎡110⎦⎤+x3⎣⎡001⎦⎤,令x2=1x_2=1x2=1,x3=1x_3=1x3=1,q=x/sum(x)=[0.33330.33330.3333]q=x/sum(x)=\begin{bmatrix}0.3333\\0.3333\\0.3333\end{bmatrix}q=x/sum(x)=⎣⎡0.33330.33330.3333⎦⎤。
可见,对于仅有平凡解的方程组,稳态向量是不存在的;对于每个随机矩阵,稳态向量若存在,则唯一。
关于马尔科夫链的收敛:
定理:如果PPP是一个n×nn\times nn×n的正则随机矩阵,则PPP具有唯一的稳态向量qqq。即,如果x0x_0x0是初始状态,且xk+1=Pxk, k=0,1,2,…x_{k+1}=Px_k,\space k=0,1,2,\dotsxk+1=Pxk, k=0,1,2,…,则当kkk趋近于无穷大的时候,马尔科夫链{xk}\{x_k\}{xk}收敛于qqq(即xkx_kxk中的元素无线接近qqq中对应的元素)。
例:
把天气状态分为好、中和差三种。如果今天天气好,明天天气有60%可能性是好的,30%可能性是中的,10%可能性是差的;如果今天天气中等,明天40%可能是好的,30%可能是中;如果今天天气差,明天40%是好的,50%是中的。
(1)求随机矩阵;
(2)若今天是好天气概率50%,中的概率50%,求明天是差天气的概率;
(3)如果周一的天气预报40%是中等天气,60%为差天气,求星期三是好天气的概率;
(4)求对应的稳态向量qqq。
解:
(1)随机矩阵WWW为:
(2)今天的概率向量为:v0=[0.500.500]v_0= \begin{bmatrix}0.50\\0.50\\0\end{bmatrix}v0=⎣⎡0.500.500⎦⎤
明天的天气的概率向量:
v1=Wv=[0.50.30.2]v_1=Wv=\begin{bmatrix}0.5\\0.3\\0.2\end{bmatrix}v1=Wv=⎣⎡0.50.30.2⎦⎤
MATLAB计算:
W =
0.6000 0.4000 0.4000
0.3000 0.3000 0.5000
0.1000 0.3000 0.1000
v =
0.5000
0.5000
0
>> Wv = W*v
Wv =
0.5000
0.3000
0.2000
所以,明天天气是差的概率是v1第三行=0.2v_{1第三行}=0.2v1第三行=0.2
(3)周一的天气概率向量:vmon=[00.40.6]v_{mon}=\begin{bmatrix}0\\0.4\\0.6\end{bmatrix}vmon=⎣⎡00.40.6⎦⎤
MATLAB计算:
>> vmon = [0;.4;.6]
vmon =
0
0.4000
0.6000
>> vtue = W*vmon
vtue =
0.4000
0.4200
0.1800
>> vwedn = W*vtue
vwedn =
0.4800
0.3360
0.1840
周二的天气概率:vtue=Wvmon=[0.40.420.18]v_{tue}=Wv_{mon}=\begin{bmatrix}0.4\\0.42\\0.18\end{bmatrix}vtue=Wvmon=⎣⎡0.40.420.18⎦⎤
周三的天气概率:vwedn=Wvtue=[0.480.3360.184]v_{wedn}=Wv_{tue}=\begin{bmatrix}0.48\\0.336\\0.184\end{bmatrix}vwedn=Wvtue=⎣⎡0.480.3360.184⎦⎤
所以,周三是好天气的概率为vwedn第一行=48%v_{wedn第一行}=48\%vwedn第一行=48%。
(4)求对应的稳态向量qqq:
稳态向量满足方程Wv=vWv=vWv=v,即(W−I)v=0(W-I)v=0(W−I)v=0,把增广矩阵[W−IO][W-I\quad O][W−IO]化简为简化阶梯阵:
>> W
W =
0.6000 0.4000 0.4000
0.3000 0.3000 0.5000
0.1000 0.3000 0.1000
>> I
I =
1 0 0
0 1 0
0 0 1
>> WO = [W-I O]
WO =
-0.4000 0.4000 0.4000 0
0.3000 -0.7000 0.5000 0
0.1000 0.3000 -0.9000 0
>> rref(WO)
ans =
1 0 -3 0
0 1 -2 0
0 0 0 0
通解为:x=[x1x2x3]=[3x32x3x3]=x3[321]x=\begin{bmatrix}x_1\\x_2\\x_3\end{bmatrix}=\begin{bmatrix}3x_3\\2x_3\\x_3\end{bmatrix}=x_3\begin{bmatrix}3\\2\\1\end{bmatrix}x=⎣⎡x1x2x3⎦⎤=⎣⎡3x32x3x3⎦⎤=x3⎣⎡321⎦⎤
任选一个特殊解,令x3=1x_3=1x3=1,则x特=[321]x_特=\begin{bmatrix}3\\2\\1\end{bmatrix}x特=⎣⎡321⎦⎤
所以,对应随机矩阵WWW的稳态向量q=13+2+1[321]=[121316]q=\frac{1}{3+2+1}\begin{bmatrix}3\\2\\1\end{bmatrix}=\begin{bmatrix}\frac{1}{2}\\\frac{1}{3}\\\frac{1}{6}\end{bmatrix}q=3+2+11⎣⎡321⎦⎤=⎣⎡213161⎦⎤
由此稳态向量可看出,经过长期的迭代,在某一天中好天气的概率为50%,中的天气33.33%,差天气16.67%。