为何qqnorm生成的理论分位数与小数据集手动计算结果存在差异？

问题描述 #

我用R的qqnorm()函数从整数数据集生成理论分位数，原本以为不管数据集大小n是多少，这个函数算出的分位数和手动用公式qnorm((1:n - 0.5)/n)得到的结果都一致，但实际发现只有n=1和n>10时两者才相同。

复现代码 #

# Define the function to compare quantiles
compare_quantiles <- function(n) {
  data <- 1:n
  qq_result <- qqnorm(data, plot = FALSE)
  theoretical_quantiles_qq_norm <- qq_result$x
  theoretical_quantiles_by_hand <- qnorm((1:length(data) - 0.5) / length(data))
  
  # Check if the quantiles are identical
  identical_results <- all.equal(theoretical_quantiles_qq_norm, theoretical_quantiles_by_hand)
  
  return(identical_results)
}
# Check for n = 1 to n = 100
results <- sapply(1:100, compare_quantiles)

原因解释 #

其实这是qqnorm()函数的内置规则导致的——它会根据样本量n的大小切换理论分位数的计算方式：

• 当n ≤ 10时，用的是Blom变换：qnorm((rank - 3/8)/(n + 1/4))，不是你手动用的(rank - 0.5)/n；
• 当n > 10时，它才会切换到(rank - 0.5)/n的计算逻辑，这时候就和你的手动结果对上了；
• 至于n=1，因为只有一个数据点，不管哪种公式算出来的分位点都一样，所以结果自然相同。

你可以查看qqnorm()的底层代码（调用stats:::qqnorm.default）验证这一点，里面有明确的分支判断：当n ≤ 10时使用Blom公式，否则使用(i - 0.5)/n。

验证代码 #

如果要手动复现n≤10时qqnorm()的结果，可以改用Blom变换：

theoretical_quantiles_blom <- qnorm((1:length(data) - 3/8)/(length(data) + 1/4))

此时对于n≤10的情况，这个结果会和qqnorm()$x完全一致。

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