R实现指数分布污染样本置信区间平均幅度计算及绘图

1. 两类样本的置信区间构造 #

采用大样本正态近似方法构造期望值倒数（即参数λ）的95%近似置信区间，逻辑如下：

• 对任意长度为n的样本，先计算样本均值$\bar{x}$，指数分布下样本均值渐近服从正态分布，通过delta方法推导可得λ的估计量$1/\bar{x}$同样渐近正态
• 95%置信水平对应的标准正态分位数为z = qnorm(0.975) ≈ 1.96，置信区间形式为：

  $\frac{1}{\bar{x}} \pm z \cdot \frac{1}{\bar{x}\sqrt{n}}$

• 无污染样本直接通过rexp(n, rate=2.5)生成，服从Exp(λ=2.5)分布
• 污染样本生成规则：先生成原始无污染样本，随机抽取占比20%的观测位置，将对应值替换为Exp(λ_c=0.02)分布的生成值即可
• 单个置信区间的幅度（长度）为上下界差值，即2*z/(bar_x * sqrt(n))

2. 分样本量计算平均区间幅度 #

固定随机种子为740保证结果可复现，n的取值范围为100到2500、步长100，每个n下生成m=1050个两类样本，分别计算所有样本的置信区间幅度后取均值，对应得到MA(n)和MAc(n)，实现代码如下：

library(ggplot2)
# 固定全局参数
set.seed(740)
m <- 1050
lambda_raw <- 2.5
lambda_contam <- 0.02
contam_ratio <- 0.2
z <- qnorm(0.975)
n_list <- seq(100, 2500, 100)

# 初始化结果存储表
res <- data.frame(
  n = integer(0),
  ma = numeric(0),
  sample_type = character(0)
)

# 遍历所有样本量计算
for (n in n_list) {
  contam_num <- round(contam_ratio * n)
  # 计算无污染样本的平均区间幅度
  width_clean <- replicate(m, {
    x <- rexp(n, lambda_raw)
    2 * z / (mean(x) * sqrt(n))
  })
  # 计算污染样本的平均区间幅度
  width_contam <- replicate(m, {
    x <- rexp(n, lambda_raw)
    replace_idx <- sample(1:n, contam_num, replace = F)
    x[replace_idx] <- rexp(contam_num, lambda_contam)
    2 * z / (mean(x) * sqrt(n))
  })
  # 结果存入表
  res <- rbind(res,
               data.frame(n = n, ma = mean(width_clean), sample_type = "无污染MA(n)"),
               data.frame(n = n, ma = mean(width_contam), sample_type = "污染MAc(n)")
  )
}

3. 对比可视化绘制 #

以n为横轴、平均区间幅度为纵轴，按样本类型分组绘制折线+散点图即可实现对比，代码如下：

ggplot(res, aes(x = n, y = ma, color = sample_type)) +
  geom_line(linewidth = 1) +
  geom_point(size = 1.8) +
  labs(
    x = "样本量n",
    y = "95%置信区间平均幅度",
    color = "样本类型",
    title = "污染/无污染样本置信区间平均幅度随样本量变化对比"
  ) +
  theme_bw()

运行代码后可直接得到计算结果与对比图：随着样本量n增大，两类样本的平均区间幅度都会单调下降；由于混入的异常值来自均值为50的指数分布（远高于原分布均值0.4），会拉高样本均值，导致污染样本用普通正态近似计算得到的区间幅度明显低于无污染样本，存在区间宽度被低估的问题。

内容的提问来源于stack exchange，提问作者CC-SAM