R语言硬币翻转相关概率Bootstrap抽样代码优化与正确性验证

你的Bootstrap抽样步骤逻辑是合理的：

• 有放回抽样学生：符合Bootstrap核心思想，通过重抽样模拟样本分布
• 随机截取学生翻转序列：保留了序列的前后依赖特性，不会破坏翻转间的关联关系
• 计算条件概率与置信区间：通过多次重抽样得到条件概率的分布，可用于判断是否存在记忆效应（无记忆时各条件概率应为0.5）

1. 生成模拟测试数据 #

先构造符合场景的模拟数据，用于验证代码逻辑：

# 生成单个学生的带记忆硬币翻转序列
generate_student_seq <- function(n_flips) {
  seq <- character(n_flips)
  seq[1] <- sample(c("H", "T"), 1, prob = c(0.5, 0.5))
  for (i in 2:n_flips) {
    prev <- seq[i-1]
    seq[i] <- sample(c("H", "T"), 1, prob = ifelse(prev == "H", c(0.6, 0.4), c(0.4, 0.6)))
  }
  seq
}

# 生成100名学生的数据，每人翻转次数在20-50间随机
set.seed(123)
student_data <- lapply(sample(20:50, 100, replace = TRUE), generate_student_seq)

2. 高效Bootstrap核心实现 #

优化点说明： #

• 预分配结果矩阵，减少内存动态分配开销
• 用data.table加速分组统计，替代低效循环
• 过滤无效短序列，避免计算错误

library(data.table)

bootstrap_transition_probs <- function(data, k = 1000) {
  n_students <- length(data)
  # 预分配结果矩阵，存储每次Bootstrap的4个条件概率
  results <- matrix(NA, nrow = k, ncol = 4, 
                    dimnames = list(NULL, c("P(H|H)", "P(T|H)", "P(H|T)", "P(T|T)")))
  
  for (b in 1:k) {
    # 1. 有放回抽样学生
    sampled_idx <- sample(1:n_students, n_students, replace = TRUE)
    sampled_students <- data[sampled_idx]
    
    # 2. 随机截取序列（长度至少为2才能计算条件概率）
    truncated_seqs <- lapply(sampled_students, function(seq) {
      len <- length(seq)
      if (len < 2) return(NULL)
      start <- sample(1:(len-1), 1)
      end <- sample(start:len, 1)
      seq[start:end]
    })
    valid_seqs <- Filter(function(x) length(x) >= 2, truncated_seqs)
    
    # 3. 提取所有相邻翻转对
    transition_pairs <- do.call(rbind, lapply(valid_seqs, function(seq) {
      data.table(prev = seq[-length(seq)], curr = seq[-1])
    }))
    
    # 4. 计算条件概率
    counts <- transition_pairs[, .N, by = .(prev, curr)]
    total_prev <- counts[, .(total = sum(N)), by = prev]
    probs <- counts[total_prev, on = "prev", prob = N / total]
    
    # 填充结果
    results[b, ] <- c(
      probs[prev == "H" & curr == "H", prob],
      probs[prev == "H" & curr == "T", prob],
      probs[prev == "T" & curr == "H", prob],
      probs[prev == "T" & curr == "T", prob]
    )
  }
  
  # 计算95%置信区间
  ci <- apply(results, 2, function(x) quantile(x, c(0.025, 0.975), na.rm = TRUE))
  list(bootstrap_results = results, confidence_intervals = ci)
}

# 运行Bootstrap（k=1000为例）
set.seed(456)
bootstrap_output <- bootstrap_transition_probs(student_data, k = 1000)

# 查看置信区间
print(bootstrap_output$confidence_intervals)

3. 大规模并行优化 #

若k值较大（如10000），可通过并行处理进一步提速：

library(parallel)
n_cores <- detectCores() - 1

# 单轮Bootstrap处理函数
single_bootstrap <- function(data, n_students) {
  sampled_idx <- sample(1:n_students, n_students, replace = TRUE)
  sampled_students <- data[sampled_idx]
  
  truncated_seqs <- lapply(sampled_students, function(seq) {
    len <- length(seq)
    if (len < 2) return(NULL)
    start <- sample(1:(len-1), 1)
    end <- sample(start:len, 1)
    seq[start:end]
  })
  valid_seqs <- Filter(function(x) length(x) >= 2, truncated_seqs)
  
  transition_pairs <- do.call(rbind, lapply(valid_seqs, function(seq) {
    data.table(prev = seq[-length(seq)], curr = seq[-1])
  }))
  
  counts <- transition_pairs[, .N, by = .(prev, curr)]
  total_prev <- counts[, .(total = sum(N)), by = prev]
  probs <- counts[total_prev, on = "prev", prob = N / total]
  
  c(
    probs[prev == "H" & curr == "H", prob],
    probs[prev == "H" & curr == "T", prob],
    probs[prev == "T" & curr == "H", prob],
    probs[prev == "T" & curr == "T", prob]
  )
}

# 并行运行
set.seed(789)
cl <- makeCluster(n_cores)
clusterExport(cl, c("student_data", "single_bootstrap"))
clusterEvalQ(cl, library(data.table))

bootstrap_results_parallel <- parSapply(cl, 1:1000, function(x) {
  single_bootstrap(student_data, length(student_data))
})
stopCluster(cl)

# 转换结果并计算置信区间
bootstrap_results_parallel <- t(bootstrap_results_parallel)
colnames(bootstrap_results_parallel) <- c("P(H|H)", "P(T|H)", "P(H|T)", "P(T|T)")
ci_parallel <- apply(bootstrap_results_parallel, 2, function(x) quantile(x, c(0.025, 0.975), na.rm = TRUE))
print(ci_parallel)

得到置信区间后：

• 若P(H|H)的95%置信区间不包含0.5，说明前一次为正面时，下一次正面的概率显著偏离无记忆状态
• 若P(T|T)的95%置信区间不包含0.5，同理说明前一次为反面时，下一次反面的概率显著偏离无记忆状态
• 只要其中一个置信区间不包含0.5，即可支持“前一次翻转结果影响下一次”的假设

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