如何在R中构建二阶及更高阶马尔可夫链的可复现示例？

一、数据准备 #

首先定义数值向量并转换为状态序列，最终整理为数据框：

# 定义数值向量
data <- c(2.3, 1.7, 1.7, 1.4, 1.7, 2.0, 2.2, 1.8, 1.7, 1.6, 2.0, 1.5, 1.1,
          1.4, 1.8, 2.0, 1.5, 1.2, 1.0, 1.2, 1.5, 1.6, 1.1,
          1.5, 2.0, 2.1, 2.1, 2.0, 1.7, 1.7, 1.7, 1.3, 0.8, -0.1, 0.0, -0.1, -0.2, 0.0, 0.1, 0.2, 0.2, 0.0, 0.2, 0.5, 0.7, 1.4,
          1.0, 0.9, 1.1, 1.0, 1.0, 0.8, 1.1, 1.5, 1.6, 1.7, 2.1, 2.5, 2.7, 2.4, 2.2, 1.9, 1.6, 1.7, 1.9, 2.2, 2.0, 2.2, 2.1,
          2.1, 2.2, 2.4, 2.5, 2.8, 2.9, 2.9, 2.7, 2.3, 2.5, 2.2, 1.9, 1.6, 1.5, 1.9, 2.0, 1.8, 1.6, 1.8, 1.7, 1.7, 1.8, 2.1,
          2.3, 2.5, 2.3, 1.5, 0.3, 0.1, 0.6, 1.0, 1.3, 1.4, 1.2, 1.2, 1.4, 1.4, 1.7, 2.6, 4.2, 5.0, 5.4, 5.4, 5.3, 5.4, 6.2,
          6.8, 7.0, 7.5, 7.9, 8.5, 8.3, 8.6, 9.1, 8.5, 8.3, 8.2, 7.7, 7.1, 6.5, 6.4, 6.0, 5.0, 4.9, 4.0, 3.0, 3.2, 3.7, 3.7)

# 转换为状态序列：Increase（上升）、Decrease（下降）、Stagnant（平稳）
direction <- c(0, diff(data))
direction_states <- ifelse(direction > 0, "Increase", ifelse(direction < 0, "Decrease", "Stagnant"))

# 整理为数据框并移除初始无状态差异的行
df_track <- as.data.frame(cbind(data, direction_states))
df_track <- df_track[-1, ]

二、一阶马尔可夫链实现 #

使用markovchain包构建一阶模型，代码及输出如下：

library(markovchain)

# 生成状态转移对
direction_states_diff <- df_track$direction_states[2:nrow(df_track)]
prev_states <- df_track$direction_states[1:length(direction_states_diff)]

# 计算转移矩阵
transition_table <- table(prev_states, direction_states_diff)
transition_matrix <- transition_table / rowSums(transition_table)
transition_matrix[is.na(transition_matrix)] <- 0

# 创建马尔可夫链对象
markov_chain <- new("markovchain", 
                    states = c("Decrease", "Increase", "Stagnant"), 
                    transitionMatrix = as.matrix(transition_matrix), 
                    name = "一阶状态马尔可夫链")

print(markov_chain)

输出结果：

Decrease  Increase   Stagnant
Decrease 0.5614035 0.3333333 0.10526316
Increase 0.2537313 0.6567164 0.08955224
Stagnant 0.5833333 0.3333333 0.08333333

模拟后续状态：

set.seed(123) 
simulated_states <- rmarkovchain(n = 10, object = markov_chain, t0 = tail(df_track$direction_states, 1))
print(simulated_states)

三、高阶马尔可夫链的问题 #

直接尝试构建二阶马尔可夫链时，由于markovchain包的transitionMatrix仅支持二维矩阵，若传入三维数组会触发错误：

Error in validObject(.Object) : 
  invalid class “markovchain” object: invalid object for slot "transitionMatrix" in class "markovchain": got class "array", should be or extend class "matrix"

四、解决方案：复合状态转换法 #

高阶马尔可夫链可以通过将连续k个状态合并为一个复合状态，转化为一阶马尔可夫链处理。以二阶为例：

1. 生成二阶复合状态序列 #

将相邻两个原始状态拼接为复合状态，比如"Decrease-Increase"表示前一步是Decrease、当前是Increase：

# 生成二阶复合状态：前i-1和第i个状态的组合
k <- 2 # 阶数
compound_states <- sapply(1:(length(direction_states)-k+1), function(i) {
  paste(direction_states[i:(i+k-1)], collapse = "-")
})

2. 计算复合状态的转移矩阵 #

复合状态的转移是从S_t-S_{t+1}到S_{t+1}-S_{t+2}，本质是一阶转移：

# 生成复合状态的转移对
prev_compound <- compound_states[1:(length(compound_states)-1)]
next_compound <- compound_states[2:length(compound_states)]

# 计算转移矩阵
compound_transition_table <- table(prev_compound, next_compound)
compound_transition_matrix <- compound_transition_table / rowSums(compound_transition_table)
compound_transition_matrix[is.na(compound_transition_matrix)] <- 0

3. 创建二阶马尔可夫链（复合状态版） #

# 获取所有唯一的复合状态
unique_compound_states <- unique(compound_states)

# 构建复合状态马尔可夫链
compound_markov <- new("markovchain",
                       states = unique_compound_states,
                       transitionMatrix = as.matrix(compound_transition_matrix),
                       name = "二阶复合状态马尔可夫链")

print(compound_markov)

4. 模拟并还原原始状态序列 #

模拟复合状态后，拆分出每个时间步的原始状态：

set.seed(123)
# 以最后一个复合状态为初始状态，模拟10个复合状态
simulated_compound <- rmarkovchain(n = 10, object = compound_markov, 
                                   t0 = tail(compound_states, 1))

# 还原原始状态：初始状态取最后k-1个原始状态，后续从复合状态中提取第二个元素
original_last_k <- tail(direction_states, k-1)
simulated_original <- c(original_last_k, sapply(simulated_compound, function(x) strsplit(x, "-")[[1]][2]))

print(simulated_original)

扩展到更高阶 #

若要构建k阶马尔可夫链，仅需调整k值，重复上述步骤：

• 将连续k个原始状态拼接为复合状态
• 计算复合状态的一阶转移矩阵
• 模拟后拆分复合状态得到原始状态序列

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