装备升级问题的Markov解法
理论预备
一、一般随机过程
1、 定义：
依赖于一个变动参数t的一族随机变量{X(t), t ∈ T}。其中，变动参数t的所有可取值
的集合T为参数空间。X(t)的值所构成的集合S成为随机过程的状态空间。
例如，从时间t = 0开始记录某电话总机的呼叫次数，设t = 0时没有呼叫，至时刻t的
呼叫次数记作Nt，则随机变量族{Nt, t ≥ 0}是随机过程。
2、 马氏过程：
如果已知在时间t系统处于状态X的条件下，在时刻τ(τ > 𝑡𝑡)系统所处状态与时刻t以
前系统所处的状态无关，此过程便为马尔可夫过程（随机过程的一个子类）。
例如，在布朗运动中，已知时刻t下的运动状态条件下，微粒在t后的运动情况和微
粒在t以前的情况无关。若X(t)表示微粒在时刻t的位置，则X(t)是马尔可夫过程。
二、马尔可夫链
1、 定义：
设{Xn, n = 1,2, … }是一个随机变量序列，用“Xn = i”表示时刻n系统处于状态i这一事
件，称pij (n) = p(Xn+1 = j ∣ Xn = i)为事件“Xn = i”出现的条件下，事件“Xn+1 = j”出
现的概率，又称它为系统的一步转移概率。若对任意的非负整数i1、i2、i3 … in−1、i、j以及一切n ≥ 0，有 p(Xn+1 = j ∣ Xn = i, Xk = ik, k = 1,2, … n − 1) = p(Xn+1 = j ∣ Xn = i) = pij (n)
则称{Xn}是一个马尔可夫链。一步转移概率有以下的性质： pij ≥ 0，(i, j = 1,2, … , n)
�pij = 1 ，nj=1 (i, j = 1,2, … , n)
把各个状态之间的一步转移概率排成矩阵，成为状态矩阵
P = �p11 ⋯ p1n⋮ ⋱ ⋮pn1 ⋯ pnn�
每个状态i对应状态矩阵P的第i行。
三、k步转移概率与k步转移矩阵 1、 k步转移概率
系统从状态i恰好经过k步转移到状态j的概率。记作pij(k) = p(Xk+1 = j ∣ X1 = i)。 2、 k步转移矩阵 P(k) = �pij(k)�n×n
显然，P(k)为概率矩阵，即有： pij(k) ≥ 0，(i, j = 1,2, … , n)
�pij(k) = 1 ，nj=1 (i, j = 1,2, … , n)
3、 切普曼∙柯尔莫哥洛夫方程 pij(n) = �pik(m)pkj(n−m)k P(n) = P(m)P(n−m)
应用以上方程可以推出： P(n) = P(n−1)P = Pn
问题及解决
开门见山，直接提出问题：
已知：某装备基础等级为 1，最高等级为N。当装备在等级i的状态时，进行一次升级行
为之后，到达等级j的状态的概率为pij，即一步转移概率矩阵是 P = �p11 ⋯ p1N⋮ ⋱ ⋮pN1 ⋯ pNN�
问题：此装备从等级x到等级y(1 ≤ x < 𝑦𝑦 ≤ 𝑁𝑁)，平均需要进行多少次升级（升级次数
的数学期望）？
分析思路：
首先需要计算进行k次升级到达等级y的概率。注意这个概率并不是k步转移矩阵中的那
个pxy(k)。
理由：我们要求的这个事件的停止条件是“只要升到等级y就停止”；而k步转移中的停
止条件是“只要升到次数到k就停止”，而不管k步之前是否有无到达过等级y。但是，这个
概率却可以藉由k步转移方法计算得出。
“只要升到等级y就停止，共升级了k次”的言外之意是“前k − 1次都没有到达过等级y”，
也就是说，前k − 1次的各个一步转移概率矩阵中，第y列元素必须都等于 0。
记： Ix = (0 … 1 … 0)，第 x个元素为 1
Iy =
⎝
⎜
⎛
0
⋮1
⋮0⎠⎟
⎞
，第 y个元素为 1
Ey = �1 ⋯ 0⋮ 0 ⋮0 ⋯ 1�，单位矩阵对角线上第 y个元素为 0
那么，事件“从等级x到等级y，只要升到等级y就停止，共升级了k次”的概率记作px→y(k) ，就
有： px→y(k) = Ix(PEy)k−1PIy
于是，升级次数tx→y的数学期望（即平均次数）： Rx→y = E�tx→y� = � k ∙ px→y(k) =∞k=1 Ix �� k ∙ (PEy)k−1∞k=1 �PIy
具体计算：
由于PEy是一般矩阵，所以可以使用一般的方法。将PEy进行 Jordan分解如下
PEy = AJA−1
那么
Rx→y = � k ∙ Ix(AJA−1)k−1PIy∞k=1 = IxA�� k ∙ Jk−1∞k=1 �A−1PIy
由于J的特殊结构，大多数情况下
� k ∙ Jk−1∞k=1
是可以手算出结果的，即便不能也可以由电脑计算。
如果I − PEy可逆且PEy的所有特征根的模都< 1，那么有以下计算Rx→y的更简单方法些：
记
S(n) = �� k ∙ (PEy)k−1nk=1 �
两边都乘PEy
S(n)PEy = �� k ∙ (PEy)knk=1 �
上面两式相减，得
S(n)(I − PEy) = ��(PEy)k−1nk=1 � − n(PEy)n
两边都乘I − PEy
S(n)(I − PEy)(I − PEy) = ��(PEy)k−1nk=1 � (I − PEy) − n(PEy)n(I − PEy)
即 S(n)(I − PEy)(I − PEy) = I + n(PEy)n+1 − (n + 1)(PEy)n
由于前面假设I − PEy可逆，所以 S(n) = [I + n(PEy)n+1 − (n + 1)(PEy)n][(I − PEy)−1]2
由PEy的所有特征根的模都< 1，得(PEy)n收敛，所以 n(PEy)n+1 − (n + 1)(PEy)n → 0
所以 S(∞) = [(I− PEy)−1]2
从而 Rx→y = Ix[(I− PEy)−1]2PIy
举例模拟
例子：承接上面的问题，给出一步转移矩阵如下
P =
⎝
⎜
⎛
0
0 0
0 0 0
0 0 0 1 ⎠⎟
⎞
计算由 1级到 5级，平均需要的升级次数。
PEy =
⎝
⎜
⎛
0
0 0
0 0 0
0 0 0 0 0⎠⎟
⎞
计算可得PEy最大特征根的模= < 1，且
I − PEy =
⎝
⎜
⎛
−
−
−
− − − 0
0 − − 0
0 0 − 0
0 0 0
− 1 ⎠⎟
⎞
其行列式|I − PEy| = ≠ 0，即I − PEy可逆。
所以针对这个问题，有四种解决方法（因计算量庞大，故借助软件，这里使用的是
Matlab）：
1、基础的循环方法；
2、非 Jordan分解的 Markov方法；
3、用 Jordan分解的 Markov方法；
4、逆矩阵方法。
Matlab程序代码：
clc;clear all;
format long g;
P=[ 0 0 0
0 0
0
0 0 0 0 1];
disp('======================================================');
x=input('输入基础等级: ');
y=input('输入目标等级: ');
disp('======================================================');
Ey=eye(5,5);Ey(y,y)=0;
Ix=zeros(1,5);Ix(x)=1;
Iy=zeros(5,1);Iy(y)=1;
%一、循环方法=====================
tic;
N=10000;
total_sum=0;
%===概率判定矩阵，用以判断升级结果指向
T=zeros(5,6);
for j=2:6
i=1;
while i<=j-1
T(:,j)=T(:,j)+P(:,i);
i=i+1;
end
end
for ii=1:N
n=x;
m=y;
sum=0;
while n<m
a=rand(1);
for i=1:5
if a>=T(n,i)&&a<T(n,i+1)
n=i;
sum=sum+1;
break;
end
end
end
total_sum=total_sum+sum;
end
fprintf('一、基础的循环方法结果：%g次；用时：%g秒。\n',total_sum/N,toc);
disp('======================================================');
% 二、非Jordan分解的Markov方法====
tic;
R=0;
for n=1:1000
R=R+n*Ix*(P*Ey)^(n-1)*P*Iy;
end
fprintf('二、非Jordan分解的Markov方法结果：%g次；用时：%g秒。\n',R,toc);
disp('======================================================');
%三、进行Jordan分解的Markov方法====
tic;
R=0;
[A,J]=jordan(P*Ey);
for n=1:100
R=R+n*Ix*A*J^(n-1)*inv(A)*P*Iy;
end
fprintf('三、用Jordan分解的Markov方法结果：%g次；用时：%g秒。\n',R,toc);
disp('======================================================');
%四、逆矩阵解法====================
tic;
R=0;
R=Ix*(inv(eye(5)-P*Ey))^2*P*Iy;
fprintf('四、逆矩阵方法结果：%g次；用时：%g秒。\n',R,toc);
disp('======================================================');
运行结果如下：
==================================================================
输入基础等级: 1
输入目标等级: 5
==================================================================
一、基础的循环方法结果：次；用时：秒。
==================================================================
二、非 Jordan分解的Markov方法结果：次；用时：秒。
==================================================================
三、用 Jordan分解的Markov方法结果：次；用时：秒。
==================================================================
四、逆矩阵方法结果：次；用时：秒。
==================================================================
需要说明的是：前三种方法是普适的，当然运算量也很大；第四种方法却有其局限性，
即要求矩阵I − PEy可逆，否则将导致错误，例如：
计算从 1级到 3级的平均升级次数的时候，由于
I − PEy =
⎝
⎜
⎛
−
−
− 0
− − − 0
0 0 1 0 0
0 0 − 0
0 0 0 − 0 ⎠⎟
⎞
其行列式|I − PEy| = 0，即I − PEy不可逆，故而不能使用第四种方法。如果使用，结果
如下：
==================================================================
输入基础等级: 1
输入目标等级: 3
==================================================================
一、基础的循环方法结果：次；用时：秒。
==================================================================
二、非 Jordan分解的Markov方法结果：次；用时：秒。
==================================================================
三、用 Jordan分解的Markov方法结果：次；用时：秒。
==================================================================
Warning: Matrix is singular to working precision.
> In Up_grade at 76
四、逆矩阵方法结果：NaN次；用时：秒。
==================================================================
从上面的结果看出，前三种方法都得出了正确结果，而在调用第四种方法的时候给出
错误提示“Warning: Matrix is singular to working precision.”意为“矩阵奇异，无法正确计算”。
Python代码：
# -*- coding: utf-8 -*-
"""
装备升级问题的处理方法：
1、基本的循环方法；
2、Markov链方法（没有使用 Jordan分解，原因是 Python里没有找到对应的函数，如
果现编写 Jordan分解函数，耗时耗力）；
3、Markov链特殊情况：逆矩阵方法。
"""
import time
import random
import numpy as np
from scipy import linalg as la
#确定并输出一步转移概率矩阵 P如下：
P=([[,,,,],
[,,,,],
[,,,,],
[,,,,],
[,,,,]])
print "=================================================="
print P
print "=================================================="
#从键盘获取基础等级 x和目标等级 y：
x=int(raw_input('Input basic grade: '))
y=int(raw_input('Input target grade: '))
print "=================================================="
Ey=(5)
Ey[y-1,y-1]=0
Ix=(5)
Ix[x-1]=1
Iy=(5).reshape(5,1)
Iy[y-1]=1
begin=()
#一、循环方法=====================
#计算并输出累计矩阵 T：
T=([[]*6]*5)
for i in range(5):
T[:,i+1]=T[:,i]+P[:,i]
print T
print "=================================================="
#模拟过程，总模拟次数为 N:
N=
total_sum=0
for i in range(int(N)):
n=x-1
m=y-1
s=0
while n<m:
a=()
for i in range(5):
if (a>=T[n,i])and(a<T[n,i+1]):
n=i
s=s+1
break
total_sum=total_sum+s
print "The result of method 1:",total_sum/N,"times."
end1=()
print 'Total time of method 1:',end1-begin,'seconds.'
print "=================================================="
#二、非 Jordan分解的Markov方法====
R=0
for i in range(1,1000):
R=R+i*Ix*(P*Ey)**(i-1)*P*Iy
print "The result of method 2:",R[0,0],"times."
end2=()
print 'Total time of method 2:',end2-end1,'seconds.'
print "=================================================="
#三、逆矩阵解法====================
I=(5)
if (I-P*Ey)!=0:
R=Ix*((I-P*Ey).I)**2*P*Iy
print "The result of method 3:",R[0,0],"times."
end3=()
print 'Total time of method 3:',end3-end2,'seconds.'
else:
print 'Matrix I-Q is singular, progress ends.'
print "=================================================="
结果：
初始等级 1，目标等级 5：
==================================================
[[ 0. 0. 0. ]
[ 0. 0. ]
[ 0. ]
[ ]
[ 0. 0. 0. 0. 1. ]]
==================================================
Input basic grade: 1
Input target grade: 5
==================================================
[[ 0. 1. 1. 1. 1. ]
[ 0. 1. 1. 1. ]
[ 0. 1. 1. ]
[ 0. 1. ]
[ 0. 0. 0. 0. 0. 1. ]]
==================================================
The result of method 1: times.
Total time of method 1: seconds.
==================================================
The result of method 2: times.
Total time of method 2: seconds.
==================================================
The result of method 3: times.
Total time of method 3: seconds.
==================================================
初始等级 1，目标等级 3：
==================================================
[[ 0. 0. 0. ]
[ 0. 0. ]
[ 0. ]
[ ]
[ 0. 0. 0. 0. 1. ]]
==================================================
Input basic grade: 1
Input target grade: 3
==================================================
[[ 0. 1. 1. 1. 1. ]
[ 0. 1. 1. 1. ]
[ 0. 1. 1. ]
[ 0. 1. ]
[ 0. 0. 0. 0. 0. 1. ]]
==================================================
The result of method 1: times.
Total time of method 1: seconds.
==================================================
The result of method 2: times.
Total time of method 2: seconds.
==================================================
Matrix I-Q is singular, progress ends.
==================================================
价值衡量
仅仅计算升级次数的数学期望是不够的，理由是，在每个等级进行一次升级操作（目
标等级自然比当前等级高）的成本不是一成不变的。如果在等级i进行一次升级操作的成本
记作ci，那么成本向量就是 c = (c1 c2 … cN )
在从等级x到等级y的过程中，如果计算出在各个等级上进行升级次数的数学期望 r = (r1 r2 … rN )
那么，从等级x到等级y的总成本期望就是
Cx→y = � ciriNi=1 = c ∙ r
很显然
Rx→y = � ciNi=1
剩下的问题就是计算r = (r1 r2 … rN )。
之前计算的事件“从等级x到等级y，只要升到等级y就停止，共升级了k次”的概率为px→y(k) ，
而且： px→y(k) = Ix(PEy)k−1PIy
那么由归一化特点，必然有
� px→y(k) = Ix ��(PEy)k−1P∞k=1 � Iy∞k=1 = 1
即在上述事件序列中，装备处于等级y的积累概率为 1，考虑到概率的意义，也即此过
程中到达等级y的次数期望。那么对于
�(PEy)k−1P∞k=1
其第y列自然全是 1。第x行中第y个元素之前的各个数值就是装备从等级x到等级y过程
中对应等级出现的次数期望。
另外，由于计数发生在升级操作之后，也即对升级结果状态计数，为了得到升级前状
态计数，只要将初始状态计数结果加 1，结束状态计数结果减 1即可。
所以
r = Ix ���(PEy)k−1P∞k=1 � + I�Ey
具体模拟计算向量r（借用之前的程序，只选 Matlab中的基本循环和非 Jordan方法）：
clc;clear all;
format short;
P=[ 0 0 0
0 0
0
0 0 0 0 1];
disp('======================================================');
x=input('输入基础等级: ');
y=input('输入目标等级: ');
disp('======================================================');
Ix=zeros(1,5);Ix(x)=1;
Ey=eye(5,5);Ey(y,y)=0;
r=[0 0 0 0 0];
tic;
N=100000;
total_sum=0;
T=zeros(5,6);
for j=2:6
i=1;
while i<=j-1
T(:,j)=T(:,j)+P(:,i);
i=i+1;
end
end
for ii=1:N
n=x;
m=y;
su=0;
while n<m
a=rand(1);
r(n)=r(n)+1;
for i=1:5
if a>=T(n,i)&&a<T(n,i+1)
n=i;
su=su+1;
break;
end
end
end
total_sum=total_sum+su;
end
disp('循环方法结果：')
r=r/N
disp('======================================================');
% 二、非 Jordan分解的Markov方法====
R=0;
for n=1:1000
R=R+(P*Ey)^(n-1)*P;
end
disp('理论方法结果：')
r=Ix*(R+eye(5))*Ey
disp('======================================================');
模拟结果：
======================================================
输入基础等级: 1
输入目标等级: 5
======================================================
循环方法结果：
r =
0
======================================================
理论方法结果：
r=
0
======================================================
当然如果(PEy)k−1可逆，也可以使用逆矩阵的方法，此处不赘述。
得到r之后，用
Cx→y = � ciriNi=1 = c ∙ r
计算最终成本即可。
QQ：707509279
