考虑行列嵌套结构的0-1矩阵Bootstrap重采样置信区间估计

问题背景 #

你需要对存在行列嵌套相关性的0-1矩阵进行Bootstrap重采样，要求满足两个核心约束：

• 每列固定抽取10次
• 每行在两列中的总抽取次数不超过2次
  最终通过重采样结果估计总体中1的比例置信区间。

实现思路 #

放弃随机打乱列顺序的步骤，改为按列依次约束抽样，确保严格满足行列限制：

1. 对第一列进行有放回重采样，严格限制每行抽取次数≤2，累计完成10次抽取
2. 基于第一列的抽取结果，计算第二列每行剩余可抽取次数（2 - 第一列该行抽取次数），再对第二列进行有放回重采样，满足剩余次数约束，累计完成10次抽取
3. 对每次重采样的结果，模拟对应列的1出现概率，计算总比例

R代码实现 #

library(mosaic)

# 1. 定义示例数据
c1 <- c(1,1,1,1,1,0,0,0,0,0)  # 列1的1占比50%
c2 <- c(1,0,0,0,0,0,0,0,0,0)  # 列2的1占比10%
d <- cbind.data.frame(c1, c2)

# 2. 定义带约束的重采样函数：按指定概率抽取，限制每行最大次数，满足总抽取量
resample_constrained <- function(probs_col, max_per_row, total_draws) {
  row_counts <- rep(0, length(probs_col))
  remaining <- total_draws
  
  while (remaining > 0) {
    # 仅对未达抽取上限的行计算抽样概率
    valid_probs <- ifelse(row_counts < max_per_row, probs_col, 0)
    valid_probs <- valid_probs / sum(valid_probs)  # 归一化概率
    
    # 抽取一行并更新计数
    selected_row <- sample(seq_along(probs_col), size = 1, prob = valid_probs)
    row_counts[selected_row] <- row_counts[selected_row] + 1
    remaining <- remaining - 1
  }
  return(row_counts)
}

# 3. 执行Bootstrap重采样
n_reps <- 10000  # 重采样次数
per_col_draws <- 10  # 每列抽取次数
p1 <- mean(c1)
p2 <- mean(c2)

bt_constrained <- numeric(n_reps)

for (i in 1:n_reps) {
  # 第一列抽样：每行最多2次，共抽10次
  counts_col1 <- resample_constrained(probs_col = c1, max_per_row = 2, total_draws = per_col_draws)
  # 模拟第一列的1出现次数：每行counts_col1次试验，每次成功概率p1
  success_col1 <- sum(rbinom(nrow(d), counts_col1, p1))
  
  # 第二列抽样：每行最多(2 - 第一列抽取次数)次，共抽10次
  max_col2 <- 2 - counts_col1
  counts_col2 <- resample_constrained(probs_col = c2, max_per_row = max_col2, total_draws = per_col_draws)
  # 模拟第二列的1出现次数
  success_col2 <- sum(rbinom(nrow(d), counts_col2, p2))
  
  # 计算本次重采样的1的比例
  bt_constrained[i] <- (success_col1 + success_col2) / (per_col_draws * 2)
}

# 查看结果
mean(bt_constrained)  # 均值应接近0.3
sd(bt_constrained)    # 标准差会比按列无约束抽样更小，符合约束下的分布集中特性

优化说明 #

• 上述函数通过循环逐步调整抽样，确保严格满足约束条件，适合小样本场景
• 如果矩阵规模较大（行数/列数多），可以改用多分类分布初始抽样+超限调整的方式，减少循环次数提升效率：
  1. 用rmultinom生成初始抽取计数
  2. 将超过上限的次数重新分配给未达上限的行
  3. 重复调整直到所有行满足约束

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