决策论——马尔科夫决策精解
马尔可夫过程（Markov process）由俄国数学家A.A.马尔可夫于1907年提出，是一类重要的随机过程，广泛应用于自然科学、社会科学、工程及机器学习等领域。其核心特性是“无后效性”，即未来的状态仅依赖于当前的状态，而与过去的状态无关。这种“记忆无关性”使得马尔可夫过程在研究复杂系统时具有极大的简化作用。
以动物群体的变化为例，假设我们研究一个森林中动物数量的变化。如果当前时刻森林中有若干只动物，那么未来动物数量的变化（如增加、减少或保持不变）仅取决于当前动物的数量，而与过去的变化路径无关。类似的动态过程在生物学中很常见，如物种繁衍、灭绝的模拟等。在现实世界中，很多现象都可以被建模为马尔可夫过程。例如，布朗运动、流行病传播、交通流量等。这些现象的未来状态均只依赖于当前状态，而非之前的演变过程。在机器学习中，典型的例子是马尔可夫决策过程（Markov Decision Process, MDP），它被用于描述序列决策问题。MDP常用于强化学习（Reinforcement Learning, RL）中，用来建模智能体（agent）与环境的交互，如自动驾驶、游戏AI等。
一、马尔科夫决策理论
马尔科夫过程中最核心的几个概念：过去，现在，将来。马尔可夫模型（Markov Model）是一种统计模型，广泛应用在语音识别，词性自动标注，音字转换，概率文法、序列分类等各个自然语言处理等应用领域。经过长期发展，尤其是在语音识别中的成功应用，使它成为一种通用的统计工具。到目前为止，它一直被认为是实现快速精确的语音识别系统的最成功的方法之一。
1.1马尔科夫链
状态空间:马尔可夫链是随机变量\(X_1,X_2,X_3,…,X_n\)所组成的一个数列，每一个变量\(X_i\)都有几种不同的可能取值，他们所有可能取值的集合，被称为“状态空间”，而\(X_n\)的值则是在时间\(n\)的状态，状态空间记为\(S\)。
转移概率(Transition Probability）：\(p_{ij}=p(x_{t+1}=j|x_t=i)\)表示在时刻\(t\)处于状态\(i\)，在下一时刻\(t+1\)处于状态\(j\)概率；\(n\)是系统所有可能的状态个数，称\(P=[p_{ij}]_{n∗n}\)为一步转移概率矩阵，对于任意\(i \in S\),都有 \(\sum_{j=1}^{n}p_{ij} = 1\)，亦即具有每一行元素和等于1的特点，这样的向量称为概率向量。
例1：假定某大学有1万学生，每人每月用1支牙膏，并且只使用“中华”牙膏与“云白”牙膏两者之一。根据本月（12月）调查，有3000人使用云白牙膏，7000人使用中华牙膏。又据调查，使用云白牙膏的3000人中，有60％的人下月将继续使用云白牙膏，40％的人将改用中华牙膏； 使用中华牙膏的7000人中， 有70％的人下月将继续使用中华牙膏，30％的人将改用云白牙膏。据此可得到如下表所示。
| 云白牙膏 | 中华牙膏 | |
|---|---|---|
| 云白牙膏 | 60％ | 40％ |
| 中华牙膏 | 30％ | 70％ |
上表中的4个概率就称为状态转移概率，而这四个转移概率组成的矩阵为
可以看出， 转移概率矩阵的一个特点是其各行元素之和为１，其经济意义是：现在使用某种牙膏的人中，将来使用各种品牌牙膏的人数百分比之和为1。
马尔可夫链可用条件概率模型来描述。设随机过程\(X_t(t\in T)\)的时间集合\(T = \{1,2,3,…\}\)，状态空间\(S = \{1,2,3,…,n\}\)，即\(X_t (t\in T)\)是时间离散、状态离散的随机过程。若对任意的整数\(t\in T\)，满足
则称\(X_t(t\in T)\)为马尔可夫链，简称马氏链。上式称为过程的马尔可夫性或无后效性。
转移概率矩阵：矩阵各元素都是非负的，并且各行元素之和等于1，各元素用概率表示，在一定条件下是互相转移的，故称为转移概率矩阵。如用于市场决策时，矩阵中的元素是市场或顾客的保留、获得或失去的概率。\(P(k)\)表示\(k\)步转移概率矩阵，转移概率只与出发状态、转移步数、到达状态相关。设\(k\)步转移矩阵为$P(k)=[p_{ij}(k)] $，其中
这样一步转移概率：\(p_{ij}(1) \leftrightarrow p_{ij}\)，可证明\(k\)步转移概率矩阵为一步转移概率矩阵的\(k\)次幂。
1.2 平稳分布
平稳分布：若存在概率向量\(x = (x_1,x_2,…,x_n)\)，使得概率矩阵\(P\)满足：\(xP = x\)，则称\(x\)为\(P\)的固定概率向量（特征向量）。特别地，若\(x\)为一状态概率向量，\(P\)为状态转移概率矩阵，则称\(x\)为对应马尔可夫链的一个平稳分布。
显然，若对任意的\(i,j\in S:\lim_{m\rightarrow +\infty}p_{ij}(m) = \pi_{j}\)，则\(\pi = (\pi_1,\pi_2,…,\pi_N)\)为稳态分布。设存在\(\pi = (\pi_1,\pi_2,…,\pi_n)\)稳态分布，则下式恒成立：\(P(k) = P(k-1)P\)，令\(k\rightarrow \infty\)就得就得\(\pi = \pi P\)
平稳分布性质：
- 若随机过程某时刻的状态概率向量为平稳分布，则称过程处于平衡状态。 一旦过程处于平衡状态，将永远处于平衡状态。
- 对于有限状态的马尔可夫链，平稳分布必定存在。特别地，当状态转移矩阵为正规概率矩阵时，平稳分布唯一。
例2：假设在某一固定选区美国国会选举的投票结果可用向量表示为
假设我们用这种类型的向量每两年记录一次美国国会选举的结果，同时每次选举的结果仅依赖于前一次选举的结果，则刻画每两年选举的向量构成的序列是一个马尔可夫链。作为一个随机矩阵的例子，取
标志为 “ D ” 的第一行中的数值刻画在一次选举中为民主党投票的人在下一次选举中将如何投票的百分比。这里我们已经假设 70%的人在下一次选举中再一次投 “ D ” 的票，20%的人将投 “ R ” 的票，10%的人将投 “ L ” 的票。 对$ P $的其他两行有类似的解释，下图给出这个矩阵的图示：
例3：社会学家们发现决定一个人的收入阶层的最重要的因素就是其父母的收入阶层。如果一个人的收入属于下层类别，那么他的孩子属于下层收入的概率是0.65，属于中层收入的概率是0.28，属于上层收入的概率 是 0.07。事实上，从父代到子代，收入阶层的变化的转移概率如下：
使用矩阵的表示方式，一步转移概率矩阵记为
假设当前这一代人处在下层、中层、上层的人的比例是概率分布向量\(π_0= [ π_0(1)，π_0(2)，π_0(3)]\)，那么他们的子女的分布比例将是\(π_1＝π_0P\)，他们的孙子代的分布比例将是\(π_2= π_1P＝π_0P^2, ….\)，第\(n\)代子孙的收入分布比例将是\(π_n= π_｛n-1｝P＝π_0P^n, ….\)。假设初始概率分布为\(π_0 =[0.21,0.68, 0.11]\)，则我们可以计算前 \(n\) 代人的分布状况如下：
| 第\(n\)代人 | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 |
|---|---|---|---|---|---|---|---|---|---|---|
| 下层 | 0.210 | 0.252 | 0.270 | 0.278 | 0.282 | 0.285 | 0.286 | 0.286 | 0.289 | 0.286 |
| 中层 | 0.680 | 0.554 | 0.512 | 0.497 | 0.490 | 0.489 | 0.489 | 0.489 | 0.488 | 0.489 |
| 上层 | 0.110 | 0.194 | 0.218 | 0.225 | 0.226 | 0.225 | 0.225 | 0.225 | 0.225 | 0.225 |
我们发现从第7代人开始，这个分布就稳定不变了，这个是偶然的吗？我们换一个初始概率分布\(π_0 =[0.75, 0.15, 0.10]\) 试试看，继续计算前\(n\) 代人的分布状况如下
| 第\(n\)代人 | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 |
|---|---|---|---|---|---|---|---|---|---|---|
| 下层 | 0.750 | 0.522 | 0.407 | 0.349 | 0.318 | 0.303 | 0.295 | 0.291 | 0.289 | 0.286 |
| 中层 | 0.150 | 0.347 | 0.426 | 0.459 | 0.475 | 0.482 | 0.485 | 0.487 | 0.488 | 0.489 |
| 上层 | 0.100 | 0.132 | 0.167 | 0.192 | 0.207 | 0.215 | 0.220 | 0.222 | 0.225 | 0.225 |
我们发现，到第9代人的时候，分布又收敛了。最为奇特的是，两次给定不同的初始概率分布，最终都收敛到概率分布\(π =[0.286,0.489,0.225]\)，也就是说收敛的行为和初始概率分布\(π_0\)无关。这说明这个收敛行为主要是由概率转移矩阵\(P\)决定的。我们计算一下\(P^n\):
我们发现，当\(n\)足够大的时候，这个\(P^n\)矩阵的每一行都是稳定地收敛到\(π = [0.286, 0.489, 0.225]\) 这个概率分布。
例4：某地区有甲、乙、丙三家药厂生产板蓝根，有1600个用户，假定在研究期间无新用户加入也无老用户退出，只有用户的转移。已知 8月份有480 户是甲厂的顾客；320 户是乙厂的顾客；800户是丙厂的顾客。9 月份，甲厂的顾客有 48 户转乙厂，96户转丙厂；乙厂的顾客有32户转甲厂，64户转丙厂；丙厂有的顾客有 64户转甲厂，32户转乙厂。假设顾客保持相同的流转，请预测
（1）这三家药厂在10月和11月的顾客人数，
（2）稳态时市场的占有率。
| 从-到 | 甲 | 乙 | 丙 | 合计 |
|---|---|---|---|---|
| 甲 | 336 | 48 | 96 | 480 |
| 乙 | 32 | 224 | 64 | 320 |
| 并丙 | 224 | 32 | 704 | 800 |
| 合计 | 432 | 304 | 864 | 1600 |
9月份的状态向量为（432，304，864），由
可预测，10月份，甲、乙、丙三家的顾客数分别为（402，291，908）。
同理，
可预测，11月份，甲、乙、丙三厂的顾客数分别为383、280、937。
求得稳态分布：\(x^*= (0.2187,0.1563,0.6250)\)
所以三家药厂在均衡时的市场占有率分别是：甲22%，乙16%，丙62%。
二、马尔可夫决策案例
2.1 案例1
考虑某地区农业收成变化的3个状态，即“丰收”、“平收”和“歉收”。记E1为“丰收”状态，E2为“平收”状态，E3为“歉收”状态。下表给出了该地区1991—2030年期间农业收成的状态变化情况。试计算该地区农业收成变化的状态转移概率矩阵。
| 年份 | 1991 | 1992 | 1993 | 1994 | 1995 | 1996 | 1997 | 1998 | 1999 | 2000 |
|---|---|---|---|---|---|---|---|---|---|---|
| 序号 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 |
| 状态 | E1 | E1 | E2 | E3 | E2 | E1 | E3 | E2 | E1 | E2 |
| 年份 | 2001 | 2002 | 2003 | 2004 | 2005 | 2006 | 2007 | 2008 | 2009 | 2010 |
| 序号 | 11 | 12 | 13 | 14 | 15 | 16 | 17 | 18 | 19 | 20 |
| 状态 | E3 | E1 | E2 | E3 | E1 | E2 | E1 | E3 | E3 | E1 |
| 年份 | 2011 | 2012 | 2013 | 2014 | 2015 | 2016 | 2017 | 2018 | 2019 | 2020 |
| 序号 | 21 | 22 | 23 | 24 | 25 | 26 | 27 | 28 | 29 | 30 |
| 状态 | E3 | E3 | E2 | E1 | E1 | E3 | E2 | E2 | E1 | E2 |
| 年份 | 2021 | 2022 | 2023 | 2024 | 2025 | 2026 | 2027 | 2028 | 2029 | 2030 |
| 序号 | 31 | 32 | 33 | 34 | 35 | 36 | 37 | 38 | 39 | 40 |
| 状态 | E1 | E3 | E2 | E1 | E1 | E2 | E2 | E3 | E1 | E2 |
#python作马尔科夫状态转移图不尽如意！
import rpy2.robjects as ro
import rpy2.robjects.packages as rpackages
from rpy2.robjects.lib import grdevices
import matplotlib.pyplot as plt
import numpy as np
# 加载必要的 R 包
base = rpackages.importr('base')
utils = rpackages.importr('utils')
# 加载必要的 R 包
utils.chooseCRANmirror(ind=1) # 选择镜像
rpackages.importr('markovchain')
rpackages.importr('tidyverse')
rpackages.importr('expm')
rpackages.importr('diagram')
# 定义 R 代码
r_code = """
# 加载所需的库
library(markovchain)
library(tidyverse)
library(expm)
library(diagram)
# 定义提供的数据
tb456 <- data.frame(
年份 = c(1991:2030),
状态 = c("E1", "E1", "E2", "E3", "E2", "E1", "E3", "E2", "E1", "E2",
"E3", "E1", "E2", "E3", "E1", "E2", "E1", "E3", "E3", "E1",
"E3", "E3", "E2", "E1", "E1", "E3", "E2", "E2", "E1", "E2",
"E1", "E3", "E2", "E1", "E1", "E2", "E2", "E3", "E1", "E2")
)
# 创建前一年的状态（滞后）
tb456 <- tb456 %>%
mutate(state1 = lag(状态))
# 删除第一行（因为第一行没有滞后值）
tb456 <- tb456[-1, ]
# 创建状态转移的交叉表
tss <- table(tb456$state1, tb456$状态)
# 返回状态转移概率矩阵，保留小数点后两位
tmA <- round(prop.table(tss, 1), 2)
# 初始状态向量，假设系统最初从状态 E2 开始
initial <- matrix(c(0, 1, 0), nrow = 1, byrow = TRUE) # E1, E2, E3 对应的初始概率
# 预测未来 n 年的函数
predict_years <- function(n) {
initial %*% (tmA %^% n)
}
# 存储20年预测结果的矩阵
mats <- sapply(1:20, predict_years) %>%
t() %>%
data.frame()
# 将结果保留小数点后两位
mats <- round(mats, 2)
print(mats)
# 绘制马尔可夫链的可视化图
png(filename="markov_chain.png")
plotmat(tmA, pos = c(1,2), lwd = 1, box.lwd = 2, cex.txt = 0.8,
box.size = 0.1, box.type = "circle", box.prop = 0.5,
box.col = "light blue", arr.length = .1, arr.width = .1,
self.cex = .6, self.shifty = -.01, self.shiftx = .15, main = "Markov Chain")
dev.off()
"""
# 执行 R 代码
ro.r(r_code)
# 使用 matplotlib 显示保存的图像
img = plt.imread("markov_chain.png")
plt.imshow(img)
plt.axis('off') # 关闭坐标轴
plt.show()
2.2 案例2
在外地旅游时，为了出行方便，许多游客会选择租车出行。某公司在北京和天津两地开展汽车租赁业务，共投入7000辆车。消费者可以从任意一个城市租车，也可以任意还车到各个城市的不同网点。需要解决的问题是如何分配最初的车辆布局可以让车辆利用率达到最大化。例如：从天津借出的车辆0.7还车在天津，0.3还车在北京。从北京借出的车辆0.6还车在北京，0.4还车在天津。可用一步转移矩阵来表示：
矩阵中的每个元素表示从某地借车还车到某地的意愿，比如第一天从某地借车第二天到某地还车的意愿，可以综合考虑两地之间的距离，出行目的等因素或直接进行调查分析，得出上面一步转移概率。具体地：\(a,b,c,d\)为四个城市，其中从\(a\)城市租借的车辆在\(a\)城市归还的比例为p[1][1]（第一个状态转第一个状态的概率）、\(b\)城市为p[1][2]、\(c\)城市为p[1][3]、\(d\)城市为p[1][4]；同理在\(b\)和\(c\)城租借的车辆在各个城市归还的比例分别表示为p[2][1]、p[2][2]、p[2][3]、p[2][4]和p[3][1]、p[3][2]、p[3][3]、p[3][4]。根据租客还车的偏好，可以统计每个租车人的换车比例，例如第一天归还的比例占0.7，第二天归还比例占0.3，等等。
利用矩阵运算可以很好处理这样的问题，计算该转移矩阵的平稳分布，可知租车最后稳定状态的概率，就可推断各个城市网点所需要的配置车辆的多少。根据平稳分布的定义，
解上面方程即得\(\pi=(0.447,0.263,0.158,0.132)\)，后面根据公司初始运行的车辆总数进行相应的分配即可。从实际问题出发，可知车辆数量一定的情况下，车辆在各个城市网点的分布最终会达到平衡，以此为基础分配最初的车辆分布情况，可提高车辆的利用率直至提高运行效率。
import matplotlib.pyplot as plt
import numpy as np
#计算稳态概率
p = np.array([[0.7,0.2,0.0,0.1],[0.2,0.5,0.2,0.1],[0.1,0.1,0.5,0.3],[0.5,0.2,0.2,0.1]])
r_1 = np.array([[0.8,0.6,0.4,0.2],[0.3,0.7,0.1,0.4],[0.4,0.3,0.8,0.2],[0.2,0.5,0.4,0.8]])
r_2 = np.array([[0.2,0.4,0.6,0.8],[0.7,0.3,0.9,0.6],[0.6,0.7,0.2,0.8],[0.8,0.5,0.6,0.2]])
amount_per_day = {}
amount = np.array([1000,300,600,800])
amount_per_day[0] = amount
amount_per_day[1] = np.dot(np.multiply(p,r_1).transpose(),amount.transpose())#第一天单独处理
for i in range(2,21):
amount_per_day[i] = np.dot(np.multiply(p,r_1).transpose(),amount_per_day[i - 1].transpose())+\
np.dot(np.multiply(p,r_2).transpose(),amount_per_day[i - 2].transpose())
amount_per_day[i] = np.array(list(map(int, amount_per_day[i][:])) ) #每次循环都取整
print(amount_per_day[20])
#各相应网点初始运行时的车辆向量为[1000,300,600,800]，运行平稳后的车辆向量为[833 490 294 245]。
#结果可视化
a = []
b = []
c = []
d = []
for j in range(21):
a.append(amount_per_day[j][0])
b.append(amount_per_day[j][1])
c.append(amount_per_day[j][2])
d.append(amount_per_day[j][3])
wang_numpy = np.arange(0,21,1)
plt.plot(wang_numpy,a,'r--')
plt.plot(wang_numpy,b,'b--')
plt.plot(wang_numpy,c,'g--')
plt.plot(wang_numpy,d,'m--')
plt.xlabel('iteraitons or days')
plt.ylabel('number of cars')
plt.xticks(np.arange(0, 21, 1))
plt.grid(True)
plt.show()
从上图中可看出，大约运行6期后车辆就处于一种稳定状态，也就是最后便于管理的状态。
总结
马尔可夫链的一个重要应用成果就是平稳分布，它是描述一个系统的长期状态。平稳分布可以用来描述一个系统的长期状态，也可以用于描述一个系统的稳定性，可以用来预测系统的未来状态。马尔科夫过程具有在给定当前知识或信息的情况下，过去（即当期以前的历史状态）对于将来（即当期以后的未来状态）是无关的。这样其概率转移矩阵就表示马尔可夫链状态的转移概率，其平稳分布反映了其变化稳定态势，从而可以利用这个平稳分布进行管理决策。平稳分布可应用于各种领域，如数据挖掘、信息检索、机器研究、自然语言处理、网络安全等。例如，在数据挖掘中，可以用马尔可夫链来分析数据，推断数据之间的联系；在机器研究中，可以用马尔可夫链来预测系统的未来状态；在自然语言处理中，可以用马尔可夫链来分析语言结构；在网络安全中，可以用马尔可夫链来模拟网络攻击，从而更好地保护网络安全。
