基于分组的重采样：R语言硬币翻转序列统计验证

我用R开展分析，现有一组学生重复抛掷带相关性硬币的数据集（已整理为final数据集），需完成以下步骤：

• 可放回随机抽取100个student_id
• 每个抽取的id对应学生的所有数据全部取出，重复抽取则数据重复出现
• 统计每个抽取id对应的硬币翻转序列（HH、HT、TH、TT）的次数，禁止统计不同抽取id间的序列转换
• 重复上述步骤多次

我编写了尝试代码，但运行耗时较长，且不确定逻辑正确性：担心同一学生多次被抽取时，前一次数据的末尾序列会与后一次数据的开头序列错误关联，影响结果准确性。

原代码 #

数据生成代码 #

library(tidyverse)

set.seed(123)
ids <- 1:100
student_id <- sort(sample(ids, 100000, replace = TRUE))
coin_result <- character(1000)
coin_result[1] <- sample(c("H", "T"), 1)

for (i in 2:length(coin_result)) {
  if (student_id[i] != student_id[i-1]) {
    coin_result[i] <- sample(c("H", "T"), 1)
  } else if (coin_result[i-1] == "H") {
    coin_result[i] <- sample(c("H", "T"), 1, prob = c(0.6, 0.4))
  } else {
    coin_result[i] <- sample(c("H", "T"), 1, prob = c(0.4, 0.6))
  }
}

#tidy up
my_data <- data.frame(student_id, coin_result)
my_data <- my_data[order(my_data$student_id),]

final <- my_data %>%
    group_by(student_id) %>%
    mutate(flip_number = row_number())

尝试实现代码 #

set.seed(123) 
ids <- 1:100 

library(dplyr)

results <- list()

for (j in 1:100) {
  selected_ids <- sample(ids, 100, replace = TRUE)
  
  resampled_data <- data.frame()
  
  for (i in 1:length(selected_ids)) {
    current_id <- selected_ids[i]
    current_data <- final %>% filter(student_id == current_id)
    current_data$resampled_id <- i
    resampled_data <- rbind(resampled_data, current_data)
  }
  
  current_result <- resampled_data %>%
    group_by(resampled_id) %>%
    summarize(Sequence = str_c(coin_result, lead(coin_result)), .groups = 'drop') %>%
    filter(!is.na(Sequence)) %>%
    count(Sequence)
  
  results[[j]] <- current_result
}

结果可视化代码（可选） #

final_result <- data.frame(iteration = 1:100, HH = numeric(100), HT = numeric(100), TH = numeric(100), TT = numeric(100))

for (i in 1:100) {
  current_result <- results[[i]]
  total_count <- sum(current_result$n)
  final_result$HH[i] <- current_result$n[current_result$Sequence == "HH"] / total_count
  final_result$HT[i] <- current_result$n[current_result$Sequence == "HT"] / total_count
  final_result$TH[i] <- current_result$n[current_result$Sequence == "TH"] / total_count
  final_result$TT[i] <- current_result$n[current_result$Sequence == "TT"] / total_count
}

library(ggplot2)

final_result_long <- final_result %>%
  pivot_longer(cols = c(HH, HT, TH, TT), names_to = "Sequence", values_to = "Probability")

ggplot(final_result_long, aes(x = iteration, y = Probability, color = Sequence)) +
  geom_line()

核心逻辑验证 #

你的原代码中，通过resampled_id分组计算序列的逻辑是正确的：每个抽取的学生数据被分配独立的resampled_id，lead()只会在同一resampled_id内计算相邻翻转的序列，不会跨id产生错误关联。但代码效率低的核心原因是循环中反复用rbind拼接数据，这在R中是非常耗时的操作。

优化后的代码 #

1. 预计算每个学生的序列统计（仅需执行一次）

提前计算每个学生的四种序列次数，后续抽样直接复用统计结果，避免每次迭代都处理原始数据：

library(tidyverse)

# 预计算每个学生的序列统计
student_seq_counts <- final %>%
  group_by(student_id) %>%
  mutate(sequence = str_c(coin_result, lead(coin_result))) %>%
  filter(!is.na(sequence)) %>%
  count(sequence) %>%
  pivot_wider(names_from = sequence, values_from = n, values_fill = 0) %>%
  ungroup()

2. 高效抽样与统计（替代嵌套循环）

用向量化操作替代循环拼接，大幅提升运行速度：

set.seed(123)
n_iterations <- 100
ids <- 1:100

results <- map_dfr(1:n_iterations, function(j) {
  # 可放回抽取100个id
  selected_ids <- sample(ids, 100, replace = TRUE)
  
  # 匹配预计算的统计数据，重复抽取则重复对应行
  sampled_stats <- student_seq_counts %>%
    slice(match(selected_ids, student_id))
  
  # 计算本次迭代的序列比例
  total_seq <- sum(sampled_stats %>% select(HH, HT, TH, TT))
  seq_probs <- sampled_stats %>%
    summarize(
      HH = sum(HH)/total_seq,
      HT = sum(HT)/total_seq,
      TH = sum(TH)/total_seq,
      TT = sum(TT)/total_seq
    ) %>%
    mutate(iteration = j)
  
  return(seq_probs)
})

3. 简化可视化代码

results_long <- results %>%
  pivot_longer(cols = c(HH, HT, TH, TT), names_to = "Sequence", values_to = "Probability")

ggplot(results_long, aes(x = iteration, y = Probability, color = Sequence)) +
  geom_line()

优化说明 #

1. 预计算复用：仅需一次计算所有学生的序列统计，后续抽样直接调用，避免重复处理原始数据
2. 向量化操作：用slice(match(...))直接获取抽样后的统计行，比循环rbind效率提升数倍
3. 逻辑严谨性：每个学生的序列统计独立计算，抽样后直接累加，完全避免跨抽取id的序列关联问题

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