如何在R中自动计算多状态马尔可夫链的累积百分比

如果你的代码目前需要手动为每个状态（如A、B、C）编写cumsum语句，可以通过以下几种方式实现自动化，无需针对单个状态硬编码：

方法1：基础R的指示矩阵法 #

通过将状态序列转换为指示矩阵（每列对应一个状态，值为1表示当前时间步处于该状态），再对每列批量计算累积百分比：

# 模拟马尔可夫链状态序列
set.seed(123)
state_seq <- sample(c("A", "B", "C"), size = 1000, replace = TRUE)

# 生成指示矩阵（自动识别所有唯一状态）
indicator_mat <- model.matrix(~ 0 + state_seq)
colnames(indicator_mat) <- unique(state_seq)

# 批量计算各状态的累积百分比
cumulative_pct <- t(apply(indicator_mat, 2, function(col) {
  cumsum(col) / seq_along(col)
}))

# 查看前5步结果
cumulative_pct[, 1:5]

方法2：tidyverse宽表处理法 #

利用tidyverse的向量化操作，自动生成所有状态的累积百分比列：

library(tidyverse)

set.seed(123)
state_data <- tibble(
  time_step = 1:1000,
  state = sample(c("A", "B", "C"), 1000, replace = TRUE)
)

# 自动生成所有状态的累积百分比
cumulative_data <- state_data %>%
  # 将每个状态转为二进制指示列
  mutate(across(state, ~map_dfc(unique(state), ~as.integer(.x == .y)) %>% 
                  set_names(unique(state)))) %>%
  # 对每个状态列计算累积百分比
  mutate(across(all_of(unique(state)), ~cumsum(.x) / time_step)) %>%
  select(time_step, everything())

# 查看前5行
head(cumulative_data, 5)

方法3：tidyverse分组累积法 #

通过分组统计，自动处理所有状态的累积计算：

library(tidyverse)

set.seed(123)
state_data <- tibble(
  time_step = 1:1000,
  state = sample(c("A", "B", "C"), 1000, replace = TRUE)
)

cumulative_data <- state_data %>%
  # 按时间步和状态计数
  count(time_step, state) %>%
  # 补全所有时间步-状态组合（缺失的计数填0）
  complete(time_step, state, fill = list(n = 0)) %>%
  # 按状态分组计算累积计数和百分比
  group_by(state) %>%
  mutate(
    cumulative_count = cumsum(n),
    cumulative_pct = cumulative_count / time_step
  ) %>%
  ungroup() %>%
  # 转换为宽表格式（每列对应一个状态的累积百分比）
  pivot_wider(names_from = state, values_from = cumulative_pct)

# 查看前5行
head(cumulative_data, 5)

以上三种方法的核心优势是：自动识别所有唯一状态，无论状态数量是3个还是更多，都无需手动修改代码添加新的cumsum语句，完全适配动态变化的状态集合。

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