基于两种总体基准分布的分层抽样实现问题咨询

核心思路 #

你需要的是同时控制两个变量边际分布的抽样方法，这类问题可以通过「交叉分层+迭代比例拟合（IPF）」解决：

• 先将gender和education的所有交叉组合作为基础分层单元
• 用IPF算法计算每个交叉层的目标样本量，确保最终样本的gender和education比例同时匹配总体基准
• 针对抽样框架中层样本量不足的情况，先全选该层样本，再把剩余待抽取的样本量重新分配到有剩余容量的层，迭代调整直到尽可能满足边际要求

R代码实现 #

以下是完整的实现代码，包含IPF样本量分配、层不足处理和抽样过程：

library(dplyr)

# 参数设置
N <- 1000  # 抽样框架大小
n <- 300   # 目标样本量

# 生成抽样框架（模拟数据）
set.seed(123) # 固定种子保证可复现
framework <- data.frame(
  id        = seq(1:N),
  gender    = sample(c("M","F"), N, replace = TRUE, prob = c(0.3, 0.7)),
  education = sample(c("1. Low", "2. Mid", "3. High"), N, replace = TRUE, prob = c(0.2, 0.3, 0.5))
)

# 总体基准分布
pop_gender <- data.frame(gender = c("M", "F"), prop = c(0.5, 0.5))
pop_education <- data.frame(education = c("1. Low", "2. Mid", "3. High"), prop = c(0.4, 0.3, 0.3))

# --------------------------
# 步骤1：构建交叉分层并统计各层框架数量
# --------------------------
cross_strata <- framework %>%
  count(gender, education, name = "stratum_size") %>%
  mutate(target = n) # 初始目标样本量先设为总样本量（后续IPF调整）

# --------------------------
# 步骤2：迭代比例拟合（IPF）计算目标样本量
# --------------------------
# 转换基准为目标频数
gender_target <- pop_gender %>% mutate(target_n = prop * n)
edu_target <- pop_education %>% mutate(target_n = prop * n)

# 迭代调整直到边际误差足够小
max_iter <- 100
tolerance <- 0.01
iter <- 1
converged <- FALSE

while(iter <= max_iter & !converged) {
  # 按gender调整
  cross_strata <- cross_strata %>%
    left_join(gender_target, by = "gender") %>%
    group_by(gender) %>%
    mutate(target = target * (target_n / sum(target))) %>%
    ungroup() %>%
    select(-target_n)
  
  # 按education调整
  cross_strata <- cross_strata %>%
    left_join(edu_target, by = "education") %>%
    group_by(education) %>%
    mutate(target = target * (target_n / sum(target))) %>%
    ungroup() %>%
    select(-target_n)
  
  # 检查收敛：边际和与目标的差异是否小于 tolerance*n
  gender_diff <- abs(cross_strata %>% group_by(gender) %>% summarize(sum_target = sum(target)) %>% left_join(gender_target, by = "gender") %>% mutate(diff = abs(sum_target - target_n)) %>% pull(diff))
  edu_diff <- abs(cross_strata %>% group_by(education) %>% summarize(sum_target = sum(target)) %>% left_join(edu_target, by = "education") %>% mutate(diff = abs(sum_target - target_n)) %>% pull(diff))
  
  if(all(gender_diff < tolerance*n) & all(edu_diff < tolerance*n)) {
    converged <- TRUE
  }
  iter <- iter + 1
}

# 转换为整数样本量（四舍五入后调整总样本量为n）
cross_strata <- cross_strata %>%
  mutate(target_n = round(target)) %>%
  mutate(target_n = ifelse(row_number() == nrow(.), n - sum(target_n[-nrow(.)]), target_n))

# --------------------------
# 步骤3：处理层样本量不足的情况
# --------------------------
cross_strata <- cross_strata %>%
  mutate(actual_n = pmin(stratum_size, target_n), # 实际抽取量：不超过层内框架数量
         remaining = target_n - actual_n) # 记录未满足的样本量

# 计算剩余需要抽取的总样本量
total_remaining <- sum(cross_strata$remaining)

# 如果有剩余，将其分配到有剩余容量的层（按当前层的剩余容量比例分配）
if(total_remaining > 0) {
  available_strata <- cross_strata %>%
    filter(stratum_size > actual_n) %>%
    mutate(available = stratum_size - actual_n)
  
  if(nrow(available_strata) > 0) {
    available_strata <- available_strata %>%
      mutate(additional = round(total_remaining * (available / sum(available)))) %>%
      mutate(additional = ifelse(row_number() == nrow(.), total_remaining - sum(additional[-nrow(.)]), additional))
    
    # 更新实际抽取量
    cross_strata <- cross_strata %>%
      left_join(available_strata %>% select(gender, education, additional), by = c("gender", "education")) %>%
      mutate(additional = replace_na(additional, 0)) %>%
      mutate(actual_n = actual_n + additional) %>%
      select(-remaining, -additional)
  } else {
    warning("部分层样本量不足且无其他可用层，无法达到目标样本量n")
    n <- sum(cross_strata$actual_n)
    message(paste("实际样本量调整为：", n))
  }
}

# --------------------------
# 步骤4：按层抽取样本
# --------------------------
selected_ids <- c()

for(i in 1:nrow(cross_strata)) {
  stratum <- cross_strata[i,]
  # 从该层中抽取actual_n个样本
  stratum_ids <- framework %>%
    filter(gender == stratum$gender, education == stratum$education) %>%
    pull(id) %>%
    sample(size = stratum$actual_n, replace = FALSE)
  
  selected_ids <- c(selected_ids, stratum_ids)
}

# 提取最终样本
final_sample <- framework %>%
  filter(id %in% selected_ids)

# --------------------------
# 验证结果
# --------------------------
cat("样本gender分布：\n")
print(prop.table(table(final_sample$gender)))
cat("\n样本education分布：\n")
print(prop.table(table(final_sample$education)))
cat("\n各交叉层实际抽取量：\n")
print(final_sample %>% count(gender, education))

关键细节说明 #

• 迭代比例拟合（IPF）：通过反复调整交叉层的目标样本量，让gender和education的边际频数逐步逼近总体基准，无需假设联合分布或条件分布一致
• 层不足处理：当交叉层的框架数量小于目标样本量时，先全选该层，再将剩余样本量按比例分配到有剩余容量的层，尽可能接近目标边际分布
• 整数调整：IPF得到的是连续值，需转换为整数样本量，最后调整总样本量确保等于目标n

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