R的极客理想系列文章，涵盖了R的思想，使用，工具，创新等的一系列要点，以我个人的学习和体验去诠释R的强大。
R语言作为统计学一门语言，一直在小众领域闪耀着光芒。直到大数据的爆发，R语言变成了一门炙手可热的数据分析的利器。随着越来越多的工程背景的人的加入，R语言的社区在迅速扩大成长。现在已不仅仅是统计领域，教育，银行，电商，互联网….都在使用R语言。
要成为有理想的极客，我们不能停留在语法上，要掌握牢固的数学，概率，统计知识，同时还要有创新精神，把R语言发挥到各个领域。让我们一起动起来吧，开始R的极客理想。
关于作者：
- 张丹，分析师/程序员/Quant: R,Java,Nodejs
- blog: http://blog.fens.me
- email: bsspirit@gmail.com
转载请注明出处：
http://blog.fens.me/r-markov/
前言
马尔科夫决策过程是强化学习的数学框架基础，通过状态、动作、奖励和状态转移概率的形式化描述，为智能体提供了与环境交互的理论模型。在强化学习中，MDP为价值函数、策略优化和贝尔曼方程等关键概念提供了理论基础，是理解Q学习、策略梯度等算法的重要前提，帮助智能体学习在不确定环境中做出最优决策。
新的方向，开始强化学习的学习。
目录
- 马尔可夫(Markov)决策过程介绍
- 用R语言实现马尔可夫链
1. 马尔可夫(Markov)决策过程介绍
马尔可夫决策过程（Markov decision process, MDP），是一种关于事件发生的概率预测方法。它是根据事件的目前状况来预测其将来各个时刻（或时期)变动状况的一种预测方法。马尔可夫决策过程包含一组交互对象，即智能体（agent）和环境，并定义了5个模型要素：状态（state）、动作（action）、策略（policy）、奖励（reward）和回报（return），其中策略是状态到动作的映射，回报是奖励随时间步的折现或积累。
马尔可夫决策过程于随机过程领域，随机过程与概率论中静态的随机现象不同，随机过程的研究对象是随时间演变的随机现象（例如天气随时间的变化、城市交通随时间的变化）。在随机过程中，随机现象在某时刻的取值是一个向量随机变量，用表示，所有可能的状态组成状态集合。随机现象便是状态的变化过程。在某时刻的状态St通常取决于t时刻之前的状态。
状态转移
马尔可夫过程举例，一个具有 6 个状态的马尔可夫过程的例子。
其中每个绿色圆圈（s1-s6）表示一个状态，每个状态都有一定概率（包括概率为 0）转移到其他状态。其中s6 被称为终止状态，它不再转出到其他状态。状态之间的虚线箭头表示状态的转移，箭头旁的数字表示该状态转移发生的概率。从每个状态出发转移到其他状态的概率总和为 1。例如，有s1有 90%概率到自己保持不变，s1有 10%概率转移到s2。而在s2又有 50%概率回到s1，s2有 50%概率转移到s3。
把上面的状态转移过程，用矩阵表示，即状态转移矩阵。其中，矩阵第i行j列的值，则代表从状态Si转移到Sj的概率。
给定一个马尔可夫过程，我们就可以从某个状态出发，根据它的状态转移矩阵生成一个状态序列（episode），这个步骤也被叫做采样（sampling）。例如，从S1出发，可以生成序列S1->S2->S3->S6，或S1->S1->S2->S3->S4->S5->S3->S6等。生成这些序列的概率和状态转移矩阵有关。
马尔可夫的奖励过程
在马尔可夫过程的基础上加入奖励函数r 和 折扣因子γ，就可以得到马尔可夫奖励过程（Markov reward process）。一个马尔可夫奖励过程由(S,P,r,γ)构成，各个组成元素的含义如下所示。
- S是有限状态的集合。
- P是状态转移矩阵。
- r是奖励函数，某个状态s的奖励r(s) 指转移到该状态时可以获得奖励的期望。
- γ是折扣因子（discount factor），γ的取值范围为[0,1)。引入折扣因子的理由为远期利益具有一定不确定性，有时我们更希望能够尽快获得一些奖励，所以我们需要对远期利益打一些折扣。接近 1 的γ更关注长期的累计奖励，接近 0 的γ更考虑短期奖励。
回报的计算：在一个马尔可夫奖励过程中，从第t时刻状态St开始，直到终止状态时，所有奖励的衰减之和称为回报Gt（Return）。
计算公式：
其中，Rt表示在时刻t获得的奖励。对马尔可夫过程添加奖励函数，构建成一个马尔可夫奖励过程。例如，进入状态s2可以得到-2奖励，表明我们不希望进入s2。进入s4可以获得最高的奖励10，但是进入s6之后奖励为0，并且此时序列也终止了。
手动计算G1的奖励，选取s1为起始状态，设置γ=0.5，采样到一条状态序列为s1->s2->s3->s6，就可以计算s1的回报G1=-2.5，得到 G1 = -1 + 0.5*(-2) + 0.5^2*(-2) = -2.5。
状态价值计算
在马尔可夫奖励过程中，一个状态的期望回报（即从这个状态出发的未来累积奖励的期望）被称为这个状态的价值（value）。所有状态的价值就组成了价值函数（value function），价值函数的输入为某个状态，输出为这个状态的价值。
计算公式：V(s)为价值，Gt为某个状态的回报。
在上式的最后一个等号中，一方面，即时奖励的期望正是奖励函数的输出，另一方面，等式中剩余部分可以根据从状态出发的转移概率得到。获得马尔可夫奖励过程中非常有名的贝尔曼方程（Bellman equation）。
求解得到各个状态的价值V(s)，
我们现在用贝尔曼方程来进行简单的验证。例如，对于状态s4来说，当γ=0.5时，则
可以发现左右两边的值是相等的，说明我们求解得到的价值函数是满足状态为时的贝尔曼方程。
参考文章：https://hrl.boyuai.com/chapter/1/马尔可夫决策过程/
2. 用R语言实现马尔可夫链
我们用R语言模拟上面的整个马尔可夫决策的计算过程，先使用手搓的方法，后面再掉包实现。
2.1 马尔可夫过程
首先，定义状态集合s，共6种状态 s1,s2,s3,s4,s5,s6。
# 6种状态
> s <- c("s1", "s2", "s3", "s4", "s5", "s6")
> s
[1] "s1" "s2" "s3" "s4" "s5" "s6"
定义转移概率矩阵p
> p <- matrix(c(
+ 0.9,0.1, 0, 0, 0, 0,
+ 0.5, 0,0.5, 0, 0, 0,
+ 0, 0, 0,0.6, 0,0.4,
+ 0, 0, 0, 0,0.3,0.7,
+ 0,0.2,0.3,0.5, 0, 0,
+ 0, 0, 0, 0, 0, 1
+ ), nrow = 6, byrow = TRUE, dimnames = list(s, s))
查看矩阵结果
> p
s1 s2 s3 s4 s5 s6
s1 0.9 0.1 0.0 0.0 0.0 0.0
s2 0.5 0.0 0.5 0.0 0.0 0.0
s3 0.0 0.0 0.0 0.6 0.0 0.4
s4 0.0 0.0 0.0 0.0 0.3 0.7
s5 0.0 0.2 0.3 0.5 0.0 0.0
s6 0.0 0.0 0.0 0.0 0.0 1.0
画图转移概率矩阵
> library(diagram)
> plotmat(p,lwd = 1, box.lwd = 1, cex.txt = 0.8,
+ box.size = 0.02, box.type = "circle", shadow.size = 0.001,
+ box.col = "light blue", arr.length=0.1, arr.width=0.1,
+ self.cex = 1.5,self.shifty = -0.01,self.shiftx = 0.04,
+ main = "Markov Chain")
模拟随机的状态转移过程。
# 模拟马尔科夫链的函数
> simulate_markov_chain <- function(transition_matrix, initial_state, n_steps, stop=NULL) {
+ states <- rownames(transition_matrix)
+ chain <- character(n_steps + 1)
+ chain[1] <- initial_state
+
+ for (i in 1:n_steps) {
+ current_state <- chain[i]
+ next_state <- sample(states, size = 1, prob = transition_matrix[current_state, ])
+ chain[i + 1] <- next_state
+ }
+
+ # 停止状态
+ if(!is.null(stop)){
+ chain<-chain[1:which(chain==stop)[1]]
+ }
+
+ return(chain)
+ }
模拟马尔科夫链的结果。
# 设置初始状态和步数
> initial_state <- "s1" # 初始状态
> n_steps <- 300 # 执行步数
> stop_state<-'s6' # 结束状态
> markov_chain <- simulate_markov_chain(p, initial_state, n_steps,stop_state)
> markov_chain
[1] "s1" "s1" "s1" "s1" "s1" "s2" "s3" "s6"
# 再次随机模拟结果
> markov_chain <- simulate_markov_chain(p, initial_state, n_steps)
> markov_chain
[1] "s1" "s1" "s1" "s2" "s3" "s4" "s5" "s4" "s6"
可视化展示马尔科夫链过程
> library(ggplot2)
> library(dplyr)
> # 准备数据用于绘图
> chain_df <- data.frame(
+ step = 1:length(markov_chain),
+ state = factor(markov_chain, levels = s)
+ )
> # 绘制马尔科夫链路径
> ggplot(chain_df, aes(x = step, y = state, group = 1)) +
+ geom_line(color = "gray") +
+ geom_point(aes(color = state), size = 3) +
+ theme_minimal()
模拟多条马尔科夫链观察收敛性
# 5条路径
> n_simulations <- 5
> n_steps <- 50
> simulations <- list()
> for (i in 1:n_simulations) {
+ simulations[[i]] <- simulate_markov_chain(p, initial_state, n_steps)
+ }
# 准备数据用于绘图
> sim_df <- data.frame(
+ step = rep(0:n_steps, n_simulations),
+ state = unlist(simulations),
+ simulation = rep(1:n_simulations, each = n_steps + 1)
+ )
# 绘制多条马尔科夫链
> ggplot(sim_df, aes(x = step, y = state, group = simulation, color = factor(simulation))) +
+ geom_line(alpha = 0.6) +
+ geom_point(size = 1.5) +
+ labs(title = "多条马尔科夫链模拟",
+ x = "时间步", y = "状态",
+ color = "模拟编号") +
+ theme_minimal()
我们模拟了5条路径的马尔科夫链
计算马尔科夫链的稳态分布
# 计算稳态分布
> calculate_steady_state <- function(transition_matrix, tolerance = 1e-6) {
+ n_states <- nrow(transition_matrix)
+ # 初始猜测
+ pi <- rep(1/n_states, n_states)
+
+ while (TRUE) {
+ pi_new <- pi %*% transition_matrix
+ if (max(abs(pi_new - pi)) < tolerance) break
+ pi <- pi_new
+ }
+
+ return(as.vector(pi_new))
+ }
> steady_state <- calculate_steady_state(p)
> names(steady_state) <- s
# 稳态分布:
> steady_state
s1 s2 s3 s4 s5 s6
1.663447e-05 1.794239e-06 1.016132e-06 7.652676e-07 2.406675e-07 9.999795e-01
# 可视化稳态分布
> barplot(steady_state, col = c("gold", "gray", "blue"),
+ main = "马尔科夫链稳态分布",
+ ylab = "概率", xlab = "状态",
+ ylim = c(0, 0.6))
2.2 马尔可夫奖励过程
奖励过程：为每个状态定义了即时奖励，模拟了带奖励的马尔科夫链，并计算了总奖励和折现奖励。
编写带奖励的马尔科夫链函数，用于奖励计算。
# 模拟带奖励的马尔科夫链
> simulate_markov_reward <- function(transition_matrix, rewards, initial_state, n_steps, gamma = 0.9, stop=NULL) {
+ states <- rownames(transition_matrix)
+ chain <- character(n_steps + 1)
+ reward_sequence <- numeric(n_steps + 1)
+ discounted_rewards <- numeric(n_steps + 1)
+
+ chain[1] <- initial_state
+ reward_sequence[1] <- rewards[initial_state]
+ discounted_rewards[1] <- (gamma^0) * rewards[initial_state]
+
+ for (i in 1:n_steps) {
+ current_state <- chain[i]
+ next_state <- sample(states, size = 1, prob = transition_matrix[current_state, ])
+ chain[i + 1] <- next_state
+ reward_sequence[i + 1] <- rewards[next_state]
+ discounted_rewards[i + 1] <- (gamma^i) * rewards[next_state]
+ }
+
+ return(list(
+ chain = chain,
+ reward = unlist(reward_sequence),
+ discounted_reward = unlist(discounted_rewards),
+ total_reward = sum(unlist(reward_sequence)),
+ total_discounted_reward = sum(unlist(discounted_rewards))
+ ))
+ }
设置模拟参数
> initial_state <- "s1" # 初始状态
> n_steps <- 100 # 步长
> gamma <- 0.9 # 折扣因子
# 奖励规则
> rewards<-data.frame("s1"=-1,"s2"=-2,"s3"=-2,"s4"=10,"s5"=1,"s6"=0)
运行带奖励的马尔科夫链函数，
> result <- simulate_markov_reward(p, rewards, initial_state, n_steps, gamma)
# 总奖励
> result$total_reward
[1] -23
# 总折现奖励
> result$total_discounted_reward
[1] -9.48585606660266
查看完整和运行结果
> result
$chain # 马尔科夫链过程
[1] "s1" "s1" "s1" "s1" "s1" "s1" "s1" "s1" "s1" "s1" "s1"
[12] "s1" "s1" "s1" "s2" "s1" "s1" "s1" "s1" "s1" "s1" "s1"
[23] "s1" "s2" "s1" "s1" "s1" "s2" "s3" "s4" "s6" "s6" "s6"
[34] "s6" "s6" "s6" "s6" "s6" "s6" "s6" "s6" "s6" "s6" "s6"
[45] "s6" "s6" "s6" "s6" "s6" "s6" "s6" "s6" "s6" "s6" "s6"
[56] "s6" "s6" "s6" "s6" "s6" "s6" "s6" "s6" "s6" "s6" "s6"
[67] "s6" "s6" "s6" "s6" "s6" "s6" "s6" "s6" "s6" "s6" "s6"
[78] "s6" "s6" "s6" "s6" "s6" "s6" "s6" "s6" "s6" "s6" "s6"
[89] "s6" "s6" "s6" "s6" "s6" "s6" "s6" "s6" "s6" "s6" "s6"
[100] "s6" "s6"
$reward # 各个步骤的奖励值
[1] -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -2 -1 -1 -1
[19] -1 -1 -1 -1 -1 -2 -1 -1 -1 -2 -2 10 0 0 0 0 0 0
[37] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[55] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[73] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[91] 0 0 0 0 0 0 0 0 0 0 0
$discounted_reward # 各个步骤的折现奖励值
[1] -1.00000000 -0.90000000 -0.81000000 -0.72900000
[5] -0.65610000 -0.59049000 -0.53144100 -0.47829690
[9] -0.43046721 -0.38742049 -0.34867844 -0.31381060
[13] -0.28242954 -0.25418658 -0.45753585 -0.20589113
[17] -0.18530202 -0.16677182 -0.15009464 -0.13508517
[21] -0.12157665 -0.10941899 -0.09847709 -0.17725876
[25] -0.07976644 -0.07178980 -0.06461082 -0.11629947
[29] -0.10466953 0.47101287 0.00000000 0.00000000
[33] 0.00000000 0.00000000 0.00000000 0.00000000
[37] 0.00000000 0.00000000 0.00000000 0.00000000
[41] 0.00000000 0.00000000 0.00000000 0.00000000
[45] 0.00000000 0.00000000 0.00000000 0.00000000
[49] 0.00000000 0.00000000 0.00000000 0.00000000
[53] 0.00000000 0.00000000 0.00000000 0.00000000
[57] 0.00000000 0.00000000 0.00000000 0.00000000
[61] 0.00000000 0.00000000 0.00000000 0.00000000
[65] 0.00000000 0.00000000 0.00000000 0.00000000
[69] 0.00000000 0.00000000 0.00000000 0.00000000
[73] 0.00000000 0.00000000 0.00000000 0.00000000
[77] 0.00000000 0.00000000 0.00000000 0.00000000
[81] 0.00000000 0.00000000 0.00000000 0.00000000
[85] 0.00000000 0.00000000 0.00000000 0.00000000
[89] 0.00000000 0.00000000 0.00000000 0.00000000
[93] 0.00000000 0.00000000 0.00000000 0.00000000
[97] 0.00000000 0.00000000 0.00000000 0.00000000
[101] 0.00000000
$total_reward # 总奖励
[1] -23
$total_discounted_reward # 总折现奖励
[1] -9.485856
可视化马尔科夫链状态序列，以上面带奖励的马尔科夫链生成的数据集result为基础。
> library(ggplot2)
> library(gridExtra)
> # 准备数据
> sim_data <- data.frame(
+ step = 0:n_steps,
+ state = factor(result$chain, levels = s),
+ reward = result$reward,
+ discounted_reward = result$discounted_reward
+ )
> # 绘制马尔科夫链状态序列状态
> p1 <- ggplot(sim_data, aes(x = step, y = state, color = state)) +
+ geom_point(size = 3) +
+ geom_line(aes(group = 1), color = "gray") +
+ labs(title = "马尔科夫链状态序列", x = "时间步", y = "状态")
# 绘制马尔科夫链状态序列奖励
> p2 <- ggplot(sim_data, aes(x = step, y = reward)) +
+ geom_col(aes(fill = state)) +
+ scale_fill_manual(values = c("Sunny" = "gold", "Cloudy" = "gray", "Rainy" = "blue")) +
+ labs(title = "即时奖励", x = "时间步", y = "奖励")
# 绘制马尔科夫链状态序列折现奖励
> p3 <- ggplot(sim_data, aes(x = step, y = discounted_reward)) +
+ geom_col(aes(fill = state)) +
+ scale_fill_manual(values = c("Sunny" = "gold", "Cloudy" = "gray", "Rainy" = "blue")) +
+ labs(title = paste0("折现奖励 (γ=", gamma, ")"), x = "时间步", y = "折现奖励")
>
> # 组合图形
> grid.arrange(p1, p2, p3, ncol = 1)
价值函数计算，通过解析解法，直接求解Bellman方程得到精确解。
# 计算状态价值函数
> calculate_state_values <- function(transition_matrix, rewards, gamma = 0.9, epsilon = 1e-6) {
+ n_states <- nrow(transition_matrix)
+ I <- diag(n_states)
+
+ # 构建线性方程组: V = R + γPV → (I - γP)V = R
+ A <- I - gamma * transition_matrix
+ b <- rewards
+
+ # 解线性方程组
+ V <- solve(A, b)
+
+ return(V)
+ }
查看参数
> p
s1 s2 s3 s4 s5 s6
s1 0.9 0.1 0.0 0.0 0.0 0.0
s2 0.5 0.0 0.5 0.0 0.0 0.0
s3 0.0 0.0 0.0 0.6 0.0 0.4
s4 0.0 0.0 0.0 0.0 0.3 0.7
s5 0.0 0.2 0.3 0.5 0.0 0.0
s6 0.0 0.0 0.0 0.0 0.0 1.0
> rewards
s1 s2 s3 s4 s5 s6
1 -1 -2 -2 10 1 0
> gamma<-0.5
> gamma
[1] 0.5
计算状态价值
> state_values <- calculate_state_values(p, rewards, gamma)
> state_values
s1 s2 s3 s4 s5 s6
-2.019502 -2.214518 1.161428 10.538093 3.587286 0.000000
本文初步探索了，马尔科夫决策过程中的基本概念，包括奖励过程、价值函数和策略优化。通过调整折扣因子γ、奖励值或转移矩阵来观察这些变化如何影响最优策略和价值函数。
本文的代码已上传到github:https://github.com/bsspirit/ml/blob/master/mavkov.r
本文有部分代码是通过Deepseek生成，再次惊叹Deepseek的强大。
转载请注明出处：
http://blog.fens.me/r-markov/
