R语言求解转移矩阵极限/唯一平稳分布并绘制可视化结果

你给出的3阶转移矩阵对应不可约非周期马尔可夫链，存在唯一平稳分布，且该平稳分布等同于极限分布，以下方法均为确定性计算，不会产生随机结果。

1. 定义转移矩阵 #

首先运行给定的矩阵定义代码：

P <- matrix(c(0.2,0.3,0.5,0.1,0.8,0.1,0.4,0.2,0.4), nrow=3, ncol=3, byrow=TRUE)

2. 求解平稳分布的两种确定性方法 #

平稳分布π满足两个核心条件：π %*% P = π、sum(π) = 1，可以通过以下两种方式求解：

方法一：线性方程组求解 #

将平稳条件转换为线性方程组直接求解，精度最高：

# 构造系数矩阵：转置P减单位矩阵，对应π(P-I)=0的约束
coef_mat <- t(P) - diag(nrow(P))
# 替换最后一行为全1，对应各状态概率和为1的约束
coef_mat[nrow(coef_mat), ] <- rep(1, ncol(coef_mat))
# 构造右端向量，前n-1项为0，最后一项为1
rhs <- rep(0, nrow(P))
rhs[length(rhs)] <- 1
# 求解线性方程组得到平稳分布
pi_stationary <- solve(coef_mat, rhs)
# 打印结果
print(pi_stationary)

方法二：特征值法求解 #

平稳分布是转移矩阵转置的特征值为1对应的左特征向量，归一化后即可得到：

# 计算P转置的特征值与特征向量
eigen_result <- eigen(t(P))
# 提取特征值最接近1的特征向量（避免数值计算误差的影响）
target_vec <- eigen_result$vectors[, which.min(abs(eigen_result$values - 1))]
# 对特征向量做归一化，使各分量和为1
pi_stationary_eigen <- Re(target_vec / sum(target_vec))
# 打印结果
print(pi_stationary_eigen)

两种方法输出的结果完全一致，均为恒定的精确值，不存在随机性。

3. 结果可视化 #

可以使用基础绘图包绘制平稳分布的柱状图，直观展示各状态的稳态概率：

# 绘制柱状图
bar_pos <- barplot(pi_stationary,
                   names.arg = paste0("状态", 1:3),
                   col = c("lightblue", "lightgreen", "lightcoral"),
                   main = "转移矩阵对应唯一平稳分布",
                   ylab = "稳态概率",
                   ylim = c(0, max(pi_stationary) * 1.1))
# 添加数值标签
text(x = bar_pos,
     y = pi_stationary + 0.02,
     labels = round(pi_stationary, 3))

注意事项 #

• 以上所有方法均为确定性数值计算，全程没有随机模拟步骤，输出结果完全恒定
• 仅当转移矩阵对应的马尔可夫链不可约非周期时，平稳分布等同于极限分布，本次你提供的转移矩阵满足该条件

内容的提问来源于stack exchange，提问作者yoyo