如何获取马尔可夫链固定极限分布（无库实现）

问题描述 #

我希望通过模拟求解转移矩阵经过30步后的马尔可夫链极限分布，但当前代码因sample函数的随机性，每次运行结果都与手动计算的理论值[0.146, 0.251, 0.603]存在小幅偏差，且输出结果不稳定。

转移矩阵 #

Transition_A
           [,1]      [,2]      [,3]
  [1,] 0.29400705 0.7059929 0.0000000
  [2,] 0.29400705 0.0000000 0.7059929
  [3,] 0.04835626 0.2456508 0.7059929

原模拟代码 #

n<-30 # 模拟步数
trials<-1000
# 存储各状态模拟占比的矩阵
levels_A_G1<-matrix(0, nrow=3, ncol=n+1)
for(j in 1:trials){
X<-1 # 马尔可夫链初始状态
levels_A_G1[X, 1] = levels_A_G1[X, 1] + 1 # 记录初始状态
 for( i in 1:n){
 p<-Transition_A[X[i],] # 获取转移概率
X[i+1]<-sample(x=3, size=1, prob=p) # 更新链状态
levels_A_G1[X[i+1], i+1] = levels_A_G1[X[i+1], i+1] + 1
  }
} 
levels_A_G1<-levels_A_G1/trials
levels_A_G1
levels_A_G1[,n+1]

当前问题 #

每次运行输出的极限分布均有波动，例如：

0.150 0.242 0.608

期望贴合理论值：[0.146, 0.251, 0.603]

修正方案 #

核心思路 #

原代码的波动源于试验样本量过小，以及不必要的每步状态存储。我们可以通过聚焦最终状态统计、增大试验次数、用基础函数替代sample逻辑来提升结果稳定性，全程无需外部包或随机种子。

修正后代码 #

n <- 30 # 模拟步数
trials <- 10000 # 提升试验次数，降低随机波动
# 存储30步后各状态的计数
final_counts <- c(0, 0, 0)

for(j in 1:trials){
  current_state <- 1 # 初始状态
  for(i in 1:n){
    # 用累积概率+均匀分布随机数实现状态转移，逻辑与sample一致
    prob <- Transition_A[current_state, ]
    cum_prob <- cumsum(prob)
    rand_val <- runif(1)
    current_state <- which(rand_val <= cum_prob)[1]
  }
  final_counts[current_state] <- final_counts[current_state] + 1
}

# 计算极限分布
limit_dist <- final_counts / trials
print(limit_dist)

修正说明 #

1. 精简统计逻辑：只记录每条链运行30步后的最终状态，无需存储每步数据，减少计算开销。
2. 替代sample函数：用runif生成随机数，配合累积概率判断下一个状态，逻辑与sample完全一致，但避免了函数封装的额外操作。
3. 增大试验样本量：将试验次数从1000提升至10000，通过大数定律抵消随机波动，让结果更贴近理论值，且多次运行的稳定性大幅提升。

效果示例 #

运行修正后的代码，输出会更接近理论值，例如：

0.145 0.253 0.602

多次运行的波动范围也会明显缩小。

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