我可以假设这个样本的（对数）正态性吗？

机器算法验证 解释 对数正态分布 QQ图
2022-03-22 05:45:47

这是我的样本的 QQ 图（注意对数 Y 轴）；：n=1000

在此处输入图像描述
正如 whuber 所指出的，这表明底层分布是左偏的（右尾较短）。

在 R 中使用shapiro.test（对数转换数据），我得到的 p 值，这意味着我们正式拒绝原假设样本在 95% 的置信水平上W=0.97185.172⋅10−13H0:the sample is normal distributed

我的问题是：这在实践中是否足以进行假设（对数）正态性的进一步分析？特别是，我想使用 Cox 和 Land 的近似方法计算相似样本均值的置信区间（在论文中描述：Zou, GY, cindy Yan Huo and Taleban, J. (2009)。简单的置信区间对数正态平均值及其与环境应用的差异。Environmetrics 20, 172–180)：

ci <- function (x) {
        y <- log(x)
        n <- length(y)
        s2 <- var(y)
        m <- mean(y) + s2 / 2
        z <- qnorm(1 - 0.05 / 2) # 95%
        #z <- qnorm(1 - 0.10 / 2) # 90%
        d <- z * sqrt(s2 / n + s2 * s2 / (2 * (n - 1)))

        return(c(exp(m - d), exp(m + d)))
}

我注意到置信区间往往以略高于实际样本均值的点为中心。例如：

> mean(x)
[1] 82.3076
> y <- log(x)
> exp(mean(y) + var(y) / 2)
[1] 91.22831

下应该是相同的。H0

1个回答

与对数正态分布相比，这些数据具有短尾，与 Gamma 分布不同：

set.seed(17)
par(mfcol=c(1,1))
x <- rgamma(500, 1.9)
qqnorm(log(x), pch=20, cex=.8, asp=1)
abline(mean(log(x)) + .1,1.2*sd(log(x)), col="Gray", lwd=2)

QQ图

然而，由于数据强烈右偏，我们可以预期最大值在估计均值及其置信区间中发挥重要作用。因此，我们应该预期对数正态 (LN) 估计量会倾向于高估均值和两个置信限。

让我们检查一下，为了比较，使用通常的估计量：即样本均值及其正态理论置信区间。请注意，通常的估计器仅依赖于样本均值的近似正态性，而不是数据的近似正态性，并且 - 对于如此大的数据集 - 可以预期效果很好。为此，我们需要对ci函数稍作修改：

ci <- function (x, alpha=.05) {
  z <- -qnorm(alpha / 2)
  y <- log(x); n <- length(y); s2 <- var(y)
  m <- mean(y) + s2 / 2
  d <- z * sqrt(s2 / n + s2 * s2 / (2 * (n - 1)))
  exp(c(mean=m, lcl=m-d, ucl=m+d))
}

这是正常理论估计的并行函数：

ci.u <- function(x, alpha=.05) {
 mean(x) + sd(x) * c(mean=0, lcl=1, ucl=-1) / sqrt(length(x)) * qnorm(alpha/2)
}

应用于这个模拟数据集，输出是

> ci(x)
   mean     lcl     ucl 
2.03965 1.87712 2.21626 
> ci.u(x)
   mean     lcl     ucl 
1.94301 1.81382 2.07219

由 1.9产生的正态理论估计值ci.u看起来更接近的真实平均值，但很难从一个数据集中判断哪个过程往往效果更好。为了找出答案，让我们模拟很多数据集：1.9

trial <- function(n=500, k=1.9) {
  x <- rgamma(n, k)
  cbind(ci(x), ci.u(x))
}
set.seed(17)
sim <- replicate(5000, trial())

我们有兴趣将输出与的真实平均值进行比较。一组直方图揭示了这方面：1.9

xmin <- min(sim)
xmax <- max(sim)
h <- function(i, ...) {
  b <- seq(from=floor(xmin*10)/10, to=ceiling(xmax*10)/10, by=0.1)
  hist(sim[i,], freq=TRUE, breaks=b, col="#a0a0FF", xlab="x", xlim=c(xmin, xmax), ...)
  hist(sim[i,sim[i,] >= 1.9], add=TRUE,freq=TRUE, breaks=b, col="#FFa0a0",
                              xlab="x", xlim=c(xmin, xmax), ...)
}
par(mfcol=c(2,3))
h(1, main="LN Estimate of Mean")
h(4, main="Sample Mean")
h(2, main="LN LCL")
h(5, main="LCL")
h(3, main="LN UCL")
h(6, main="UCL")

直方图

现在很明显，对数正态程序倾向于高估平均值和置信限，而通常的程序做得很好。我们可以估计置信区间程序的覆盖范围：

> sapply(c(LNLCL=2, LCL=5, LNUCL=3, UCL=6), function(i) sum(sim[i,] > 1.9)/dim(sim)[2])
 LNLCL    LCL  LNUCL    UCL 
0.2230 0.0234 1.0000 0.9648

这个计算说：

  • LN 下限在大约 22.3% 的时间（而不是预期的 2.5%）中无法覆盖真实平均值。

  • 通常的下限在大约 2.3% 的情况下无法覆盖真实平均值，接近预期的 2.5%。

  • LN 上限将始终超过真实平均值（而不是按预期在 2.5% 的时间内低于该平均值）。这使其成为双边 100% - (22.3% + 0%) = 77.7% 置信区间，而不是 95% 置信区间。

  • 通常的上限在 100 - 96.5 = 3.5% 的时间里无法覆盖真实平均值。这比预期值 2.5% 略大。因此，通常的限制包括双边 100% - (2.3% + 3.5%) = 94.2% 置信区间，而不是 95% 置信区间。

对数正态区间的标称覆盖率从 95% 减少到 77.7% 是可怕的。对于通常的区间，降低到 94.2% 一点也不差，这可以归因于（原始数据的，而不是对数的）偏度的影响。

我们必须得出结论，对平均值的进一步分析不应假设对数正态性。

当心！某些程序（例如预测限）对偏度的敏感性会比这些均值的置信限更敏感，因此可能需要考虑它们的偏态分布。但是，对于几乎任何预期的分析，对数正态过程似乎不太可能在这些数据上表现良好。

其它你可能感兴趣的问题
上一篇矩阵分解和 PCA 之间的区别 下一篇随机森林的特征选择和参数调整