R语言构造截尾均值置信区间的代码错误排查

问题背景 #

需要测试截尾均值（trimmed means）的覆盖概率，采用Wilcox著作给出的置信区间公式构建区间：其中s_w为Winsorised方差，γ为截尾比例系数，代码中记为alpha。初始编写的代码输出的置信区间始终包含0，导致覆盖概率计算结果恒为1，区间构造逻辑存在错误。

初始问题代码 #

sample_var <- function(data, alpha){
  n <- length(data)
  data <- sort(data)
  data_t <- data[(floor(n*alpha)+1):(n-floor(alpha*n))]
  m <- length(data_t)
  t_mean <- mean(data_t)
 sigma <-  (1/(1-2*alpha)^2)* ((1/n) *sum((data_t-t_mean)^2)+ alpha*(data_t[1]-t_mean)^2 +
                       alpha*(data_t[m]-t_mean)^2)
 sigma
}
sample_var <- Vectorize(sample_var, vectorize.args = "alpha")

    conf_int <- function(data,alpha){
      a <- floor(alpha * n)
      n <- length(data)
      df <- n-2*a-1
      data_t <- data[a:(n-a)]
      t_mean <- mean(data_t)
      t_quantile <- qt(p = alpha, df = df)
      sw <- sample_var(data = data, alpha = alpha)
      ul <- t_mean + t_quantile * sw / ((1-2*alpha)*sqrt(n))
      ll <- t_mean - t_quantile * sw / ((1-2*alpha)*sqrt(n))
     c(ll, ul)
    }

错误点定位 #

• 变量定义顺序错误：conf_int函数中先执行a <- floor(alpha * n)，之后才定义n <- length(data)，运行时会直接找不到对象n，逻辑顺序完全颠倒。
• 截尾子集索引错误：构造截尾数据data_t时使用索引a:(n-a)，和方差计算函数中的截尾索引规则不一致，正确截尾索引应为(floor(n*alpha)+1):(n-floor(alpha*n))，原索引会多保留两端待截除的观测值，当a=0时还会调用不存在的0索引。
• t分位数取值错误：构造双侧置信区间时，qt(p = alpha, df = df)取的是下侧alpha分位数，数值为负，直接导致上下限计算逻辑颠倒、区间宽度异常偏大，正确取值应为置信度对应的双侧分位数1-(1-置信水平)/2。
• Winsorised方差计算逻辑错误：原计算仅对截尾后的数据计算离均差平方和，正确计算逻辑需要先将两端被截掉的观测值替换为相邻的截尾临界值，构造完整的Winsorised样本后再计算方差，原公式漏算了被替换观测对应的离均差平方和部分。

wilcox.test方案的问题 #

尝试使用wilcox.test构造区间的代码同样存在错误，且该函数本身不适用于截尾均值区间计算：

set_seed(1)
data <- rnorm(100)
wilcox_test <- function(data, alpha){
  n <- length(alpha)
  a <- floor(alpha*n)+1
  b <- n-floor(alpha)
  wilcox.test(data, subset = data[a:b], conf.int = TRUE)
}

对应错误点：

• 函数拼写错误：R中设置随机种子的正确函数为set.seed()，不存在set_seed()写法。
• 变量定义错误：wilcox_test函数中n <- length(alpha)逻辑完全错误，n应为输入数据的长度即length(data)。
• 索引计算错误：截尾上界b的计算式写为n-floor(alpha)，漏乘样本量n，正确应为n - floor(alpha*n)。
• 函数适用性错误：wilcox.test默认计算的是样本伪中位数/中位数的置信区间，不能直接用于计算截尾均值的置信区间，这也是调整alpha参数时区间结果无变化的核心原因。

修正后可运行代码 #

# 计算Winsorised方差与截尾均值标准误
sample_var <- function(data, alpha){
  n <- length(data)
  data_sorted <- sort(data)
  a <- floor(n * alpha)
  # 构造Winsorised样本：两端截尾值替换为相邻临界值
  data_win <- data_sorted
  if(a > 0){
    data_win[1:a] <- data_sorted[a+1]
    data_win[(n-a+1):n] <- data_sorted[n-a]
  }
  # 计算截尾均值
  t_mean <- mean(data_sorted[(a+1):(n-a)])
  # 计算Winsorised样本方差（分母为n）
  s_w <- sum((data_win - mean(data_win))^2)/n
  # 计算截尾均值标准误
  se <- sqrt(s_w) / ((1 - 2*alpha) * sqrt(n))
  return(list(trimmed_mean = t_mean, s_w = s_w, se = se))
}
sample_var <- Vectorize(sample_var, vectorize.args = "alpha")

# 构造截尾均值t置信区间
conf_int <- function(data, alpha, conf_level = 0.95){
  n <- length(data)
  a <- floor(n * alpha)
  # t分布自由度
  df <- n - 2*a - 1
  var_res <- sample_var(data, alpha)
  t_mean <- var_res$trimmed_mean
  se <- var_res$se
  # 双侧t临界值
  t_crit <- qt(1 - (1-conf_level)/2, df = df)
  ll <- t_mean - t_crit * se
  ul <- t_mean + t_crit * se
  return(c(lower = ll, upper = ul))
}

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