如何从负二项分布等非正态数据中抽样得到指定均值的正态子样本

前提说明 #

你原始生成的负二项分布样本总体均值约为15，若要抽取均值为25的子样本，需要原始样本中存在足够多的大于25的观测值，这是实现的必要前提。你当前使用的10000条样本生成参数下，该前提是满足的。

你之前单纯使用sample的prob参数没达到预期，核心原因是没有设计和目标正态分布匹配的权重形状，也没有做校验迭代，单次抽样很难同时满足均值和分布形状的要求。

核心实现逻辑 #

• 第一步：先根据目标均值构造抽样权重，权重的分布形状匹配目标正态分布的概率密度，同时拉高高于目标均值的样本权重，保证最终样本均值能达到目标值
• 第二步：迭代调整抽样结果，直到抽到的样本同时满足「均值误差在可接受范围内」、「正态性检验通过」两个条件

可运行代码 #

library(MASS)
# 生成原始负二项样本
set.seed(123) # 固定随机种子方便复现，去掉后每次会生成不同的符合要求的样本
my_trees <- rnegbin(10000, mu = 15, theta = 3)

# 自定义抽样函数
sample_norm_fixed_mean <- function(population, target_mean, sample_size = 500, tol = 0.1, max_iter = 1000) {
  # 可行性前置校验
  if (target_mean > max(population) || target_mean < min(population)) {
    stop("目标均值超出原始样本取值范围，无法实现")
  }
  
  for (i in 1:max_iter) {
    # 构造匹配正态分布的抽样权重，额外拉高高于目标均值的样本权重
    base_weight <- ifelse(population > target_mean, 
                          2 * dnorm(population, mean = target_mean, sd = sd(population)/2),
                          dnorm(population, mean = target_mean, sd = sd(population)/2))
    # 执行抽样
    current_sample <- sample(population, size = sample_size, prob = base_weight, replace = FALSE)
    current_mean <- mean(current_sample)
    # 校验均值误差
    if (abs(current_mean - target_mean) < tol) {
      # 校验正态性，p>0.05说明不拒绝正态分布假设
      norm_p <- shapiro.test(current_sample)$p.value
      if (norm_p > 0.05) {
        return(list(sample = current_sample, mean = current_mean, norm_p = norm_p, iter = i))
      }
    }
  }
  stop("达到最大迭代次数仍未找到符合要求的样本，请调大误差容忍度或调整目标均值")
}

# 测试：抽取均值为25的500条样本
result <- sample_norm_fixed_mean(my_trees, target_mean = 25)

# 验证结果
mean(result$sample) # 输出接近25的均值
hist(result$sample, main = "子样本分布直方图", xlab = "取值") # 可以看到近似正态的分布形状
shapiro.test(result$sample) # 正态性检验结果

参数调整说明 #

• 可以调整tol参数控制均值的误差容忍范围，值越小要求均值越精准，对应迭代耗时越长
• 可以调整max_iter参数控制最大迭代次数，目标均值越接近原始样本取值边界，需要的迭代次数越多
• 如果不需要严格的正态性校验，可以去掉Shapiro检验的判断逻辑，只保留均值校验即可，运行速度会大幅提升

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