Markov链，稳态向量及其求法，概率向量，（正则）随机矩阵及其判别方法，马尔科夫链收敛定理

  Markov链在很多领域中被用来建立数学模型，该模型常用来描述用同一种方法进行多次实验或测量，且每次实验结果属于几个指定的可能结果之一，每次实验的结果仅仅依赖于最近的上一次实验结果。

例：实验室动物每天可以吃三种食物中的任一种。结果表明，如果第一次测试中，这个动物选择一种食物，则下一次测试中它选同样的食物的概率是50%，下一次测试它选择其他两种实物的概率都是25%。
（1）求该实验的随机矩阵。
（2）如果这个动物第一次测试选了1号食物，在第二次测试中选择2号食物的概率、第三次选择食物2和3的概率是多少？

解：

（1）随机矩阵$P$：

（2）在MATLAB中求$P^2=P\times 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

即$P^2$为：

  由矩阵$P$可求出该动物第二次测试选择各种食物的概率表：$v_1=P{v}_0=\left[\begin{array}{c}0.5000\\ 0.2500\\ 0.2500\end{array}\right]$，所以第二次选择2号食物概率为$v_{1第二行}=0.2500$。

  由矩阵$P^2$可以求出该动物第三次测试选择各种食物的概率，设第一次选了第一种食物的初始概率向量$v_0=\left[\begin{array}{c}1\\ 0\\ 0\end{array}\right]$，则第三次的概率向量为：$v_2=P^2{v}_0=\left[\begin{array}{c}0.3750\\ 0.3125\\ 0.3125\end{array}\right]$，所以第三次选择食物2和3概率分别为$v_{2第二行}=0.3125$和$v_{2第三行}=0.3125$。

概率向量，随机矩阵（Stochastic Matrix），马尔科夫链

  一个具有非负元素且各元素相加为1的向量为概率向量。由概率向量作为列向量组合而成的方阵为随机矩阵（也叫马尔科夫矩阵、转移矩阵等）。
  例如，上例中，矩阵$P$为随机矩阵，$P$的每一列为概率向量。
  马尔科夫链是一个向量序列$x_0,x_1,x_2\dots$与一个随机矩阵$P$的关系，满足：
$$x_1=P{x}_0,x_2=P{x}_1,x_3=P{x}_2,\dots$$
此关系可以用一阶差分方程来描述：
$$x_{k+1}=P{x}_k,k=0,1,2,\dots$$
  当向量在$R^n$中的一个马尔科夫链描述一个系统或实验的序列时，$x_k$中的元素分别列出了系统在$n$个可能的状态中的概率，或者实验结果可能是$n$个可能的结果之一的概率。所以$x_k$也称为状态向量。

稳态向量（Steady-State Vector）

  如果$P$是一个随机矩阵，如果一个状态向量$q$满足$Pq=q$，则此向量$q$为相对于$P$的稳态向量（也可叫平衡向量）。
  每个随机矩阵都有一个稳态向量。

例：求上例的稳态向量$x$。
根据稳态向量的定义，满足$Px=x$，即$(P-I)x=0$
用上面的随机矩阵$P$构建齐次方程：$(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]$消元是一样的，下面是对增广矩阵的消元法以解线性方程组：

>> 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=\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\end{array}\right]=\left[\begin{array}{c}{x}_3\\ x_3\\ x_3\end{array}\right]=x_3\left[\begin{array}{c}1\\ 1\\ 1\end{array}\right]$
令自由变量$x_3=1$，则此特殊解为：$w=\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\end{array}\right]=\left[\begin{array}{c}1\\ 1\\ 1\end{array}\right]$；
也可令$x_3=2$，则此特殊解为：$u=\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\end{array}\right]=\left[\begin{array}{c}2\\ 2\\ 2\end{array}\right]$。

则求$P\times \left[\begin{array}{c}1\\ 1\\ 1\end{array}\right]$和$P\times \left[\begin{array}{c}2\\ 2\\ 2\end{array}\right]$：

>> 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\times w=w$，$P\times u=u$和$P\times q=q$。
  但是根据定义，向量$q$是$Px=x$，即$(P-I)x=0$的全体解的结合的概率向量（求法：使用任意一个解，如上面的$w$，除以此解$w$每个元素之和即得到概率向量$q$，可验算$q$的各元素之和为1，且满足$P\times q=q$），所以向量$q$随机矩阵$P$的一个稳态向量。上面的$w$和$u$因为各元素和不为1，所以不能叫概率向量，所以就不能是稳态向量。

正则随机矩阵（Regular Stochastic Matrix）：

  如果一个随机矩阵的某次幂$P^k$（注意，是有某次幂就行）的每个元素仅包含严格正的元素，则此随机矩阵可称为正则随机矩阵。
  对于上面的例子的随机矩阵$P$和$P^2$，每项元素都是正数，所以可以说$P$是正则随机矩阵。

例：判断下列矩阵是否为正则随机矩阵：
$P_1=\left[\begin{array}{ccc}0.7& 0.2& 0.6\\ 0& 0.2& 0.4\\ 0.3& 0.6& 0\end{array}\right]$

$P_2=\left[\begin{array}{ccc}0.7& 0.2& 0.6\\ 0& 0.2& 0\\ 0.3& 0.6& 0.4\end{array}\right]$

$P_3=\left[\begin{array}{ccc}0& 1& 0\\ 1& 0& 0\\ 0& 0& 1\end{array}\right]$

解：
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

可见，找到了一个$P_1^2$是全部严格正的（没有负数和零元素），且$P_1$每列和为1，所以$P_1$是正则随机矩阵。

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
...

可见，$P_2$矩阵无论多少次方，都含有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

矩阵$P_3$也因为存在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

对于$P_1$，通解：$x=\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\end{array}\right]=x_3\left[\begin{array}{c}2.3\\ 0.5\\ 1\end{array}\right]$，令$x_3=1$，$q=x/sum(x)=\left[\begin{array}{c}0.6053\\ 0.1316\\ 0.2632\end{array}\right]$。

对于$P_2$，通解：$x=\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ 0\end{array}\right]$，$q=x/sum(x)$不存在。

对于$P_3$，通解：$x=\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\end{array}\right]=x_2\left[\begin{array}{c}1\\ 1\\ 0\end{array}\right]+x_3\left[\begin{array}{c}0\\ 0\\ 1\end{array}\right]$，令$x_2=1$，$x_3=1$，$q=x/sum(x)=\left[\begin{array}{c}0.3333\\ 0.3333\\ 0.3333\end{array}\right]$。

可见，对于仅有平凡解的方程组，稳态向量是不存在的；对于每个随机矩阵，稳态向量若存在，则唯一。

关于马尔科夫链的收敛：

  定理：如果$P$是一个$n\times n$的正则随机矩阵，则$P$具有唯一的稳态向量$q$。即，如果$x_0$是初始状态，且$x_{k+1}=P{x}_k,k=0,1,2,\dots$，则当$k$趋近于无穷大的时候，马尔科夫链{xk}\{x_k\}{xk​}收敛于$q$（即$x_k$中的元素无线接近$q$中对应的元素）。

例：
把天气状态分为好、中和差三种。如果今天天气好，明天天气有60%可能性是好的，30%可能性是中的，10%可能性是差的；如果今天天气中等，明天40%可能是好的，30%可能是中；如果今天天气差，明天40%是好的，50%是中的。
（1）求随机矩阵；
（2）若今天是好天气概率50%，中的概率50%，求明天是差天气的概率；
（3）如果周一的天气预报40%是中等天气，60%为差天气，求星期三是好天气的概率；
（4）求对应的稳态向量$q$。
解：
（1）随机矩阵$W$为：

（2）今天的概率向量为：$v_0=\left[\begin{array}{c}0.50\\ 0.50\\ 0\end{array}\right]$
明天的天气的概率向量：
$v_1=Wv=\left[\begin{array}{c}0.5\\ 0.3\\ 0.2\end{array}\right]$

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

所以，明天天气是差的概率是$v_{1第三行}=0.2$

（3）周一的天气概率向量：$v_{mon}=\left[\begin{array}{c}0\\ 0.4\\ 0.6\end{array}\right]$

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

周二的天气概率：$v_{tue}=W{v}_{mon}=\left[\begin{array}{c}0.4\\ 0.42\\ 0.18\end{array}\right]$
周三的天气概率：$v_{wedn}=W{v}_{tue}=\left[\begin{array}{c}0.48\\ 0.336\\ 0.184\end{array}\right]$

所以，周三是好天气的概率为$v_{wedn第一行}=48\%$。

（4）求对应的稳态向量$q$：
稳态向量满足方程$Wv=v$，即$(W-I)v=0$，把增广矩阵$[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=\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\end{array}\right]=\left[\begin{array}{c}3x_3\\ 2x_3\\ x_3\end{array}\right]=x_3\left[\begin{array}{c}3\\ 2\\ 1\end{array}\right]$

任选一个特殊解，令$x_3=1$，则$x_{特}=\left[\begin{array}{c}3\\ 2\\ 1\end{array}\right]$

所以，对应随机矩阵$W$的稳态向量$q=\frac{1}{3+2+1}\left[\begin{array}{c}3\\ 2\\ 1\end{array}\right]=\left[\begin{array}{c}\frac{1}{2}\\ \frac{1}{3}\\ \frac{1}{6}\end{array}\right]$

由此稳态向量可看出，经过长期的迭代，在某一天中好天气的概率为50%，中的天气33.33%，差天气16.67%。