如何用R语言markovchain包计算马尔可夫链各状态周期？

核心问题 #

1. 初始错误原因：markovchain 对象要求 states 参数必须是字符类型，而非整数，因此首次运行时抛出类型不匹配错误。
2. 修正后警告与错误结果：
   • 警告提示矩阵不可约（即马尔可夫链包含多个通信类），period() 函数对非不可约链仅返回0，无法给出每个状态的周期。
   • 代码中存在笔误：创建的对象是 mcP，但调用 period(mc) 时引用了未定义的 mc。

正确解决思路 #

马尔可夫链的周期是基于通信类计算的：同一通信类内的所有状态周期相同。因此需要：

1. 识别链中的所有通信类；
2. 对每个通信类单独计算周期；
3. 将每个状态映射到其所属通信类的周期。

修正代码 #

library(markovchain)

# 定义转移矩阵
P <- matrix(c(0, 0.5, 0.5, 0, 0, 0,
              1, 0, 0, 0, 0, 0,
              1, 0, 0, 0, 0, 0,
              0, 0, 0, 0.25, 0.5, 0.25,
              0, 0, 0, 0.5, 0.25, 0.25,
              0, 0, 0, 0, 0, 1), nrow = 6, ncol = 6, byrow = TRUE)

# 创建markovchain对象（状态转为字符型）
states <- as.character(1:6)
mcP <- new("markovchain", states = states, transitionMatrix = P, name = "X")

# 获取所有通信类
comm_classes <- communicatingClasses(mcP)

# 初始化周期向量
periods <- numeric(length(states))
names(periods) <- states

# 逐个计算每个通信类的周期
for (cls in comm_classes) {
  if (length(cls) == 1) {
    # 单状态类（如状态6）：自转移存在，周期为1
    periods[cls] <- 1
  } else {
    # 提取通信类对应的子转移矩阵
    sub_indices <- match(cls, states)
    sub_P <- P[sub_indices, sub_indices]
    # 创建子马尔可夫链对象
    sub_mc <- new("markovchain", states = cls, transitionMatrix = sub_P)
    # 计算该不可约子链的周期
    cls_period <- period(sub_mc)
    # 为类内所有状态赋值周期
    periods[cls] <- cls_period
  }
}

# 输出结果
cat("各状态周期:\n")
print(periods)

运行结果 #

各状态周期:
1 2 3 4 5 6 
2 2 2 1 1 1 

结果说明 #

• 通信类 {1,2,3}：状态1只能在偶数步返回自身，周期为2，因此类内所有状态周期均为2；
• 通信类 {4,5}：存在自转移（如状态4有0.25概率留在自身），周期为1；
• 通信类 {6}：自转移概率为1，周期为1。

内容的提问来源于stack exchange，提问作者Sumiya Akter Nisher