非常偏斜的分布(例如对数正态分布)不会产生准确的引导置信区间。这是一个示例,显示无论您在 R 中尝试哪种引导方法,左右尾部区域都远离理想的 0.025:
require(boot)
n <- 25
B <- 1000
nsim <- 1000
set.seed(1)
which <- c('basic', 'perc', 'norm', 'bca', 'stud')
mul <- 0; sdl <- 1.65 # on log scale
dist <- c('normal', 'lognormal')[2]
switch(dist, normal = {g <- function(x) x; mu <- mul},
lognormal = {g <- exp; mu <- exp(mul + sdl * sdl / 2)})
count <- matrix(0, nrow=length(which), ncol=2,
dimnames=list(which, c('lower', 'upper')))
stat <- function(x, j) {
## See http://www.psychology.mcmaster.ca/bennett/boot09/percentileT.pdf
x <- x[j]
m <- mean(x)
s <- sd(x)
n <- length(x)
sem <- s / sqrt(n)
m.var <- sem ^ 2
c(m, m.var)
}
for(i in 1 : nsim) {
if(i %% 100 == 0) cat(i, '')
x <- g(rnorm(n, mul, sdl))
b <- boot(x, stat, R=B)
ci <- boot.ci(b, type=which)
for(w in which) {
nam <- switch(w, perc='percent', norm='normal', basic='basic',
stud='student', bca='bca')
z <- rev(rev(ci[[nam]])[1:2])
count[w, 'lower'] <- count[w, 'lower'] + (z[1] > mu)
count[w, 'upper'] <- count[w, 'upper'] + (z[2] < mu)
}
}
cat('\n')
count / nsim
结果如下:
lower upper
basic 0.000 0.329
perc 0.003 0.257
norm 0.000 0.287
bca 0.015 0.185
stud 0.005 0.129
为了单个引导程序仍然不能提供足够准确的覆盖范围:
lower upper
basic 0.001 0.114
perc 0.005 0.093
norm 0.002 0.102
bca 0.017 0.067
stud 0.011 0.058
从对数正态分布中采样时,经验似然也无法提供准确的置信区间。
是否有一种不依赖于预先知道分布的通用方法?有没有人试图通过将数据拟合到 Tukey 广义来获得平均值的置信区间分布(这种分布是高度灵活的)?对 CDF 使用 Kolmogorov-Smirnov 置信带怎么样?计算 CDF 上限和下限的平均值会非常保守吗?如果一种方法具有广泛的适用性,我会接受一些保守主义。
为了重申目标,我正在寻找一种普遍适用的方法来获得总体平均值的置信区间,使得
- 如果原始数据分布不对称,则区间不对称
- 区间在两个尾部都有正确的覆盖(例如,两者都有 0.025 的错误概率)
- 该程序不要求分析师指定任何有关基础分布或使分布对称所需的变换的信息
请注意,中心极限定理在这里是无关紧要的;我有一个固定的小样本量,置信区间必须不对称才能在两个尾部都准确。参数化- 基于对数正态模型下的置信区间和仍然有不好的覆盖率(左尾误差 0.012,右 0.047,当两者都应该是 0.025 时)。
在继续思考这个问题时,有两种广泛的方法来概念化我想讨论的问题。
- 平均值不是一个适合于非参数推理的量,至少在需要推理的准确性时是这样。样本中位数对于任何连续分布都是有意义的,我们有一个简单的中位数精确置信区间。在大小样本中根据正态分布,中位数的置信区间为比确切的更长- 基于均值的置信区间(参见下面的代码)。也许这个 1.28 的因素是为稳健性和完全的分配自由付出的合理代价。
- 即使没有单一的 bootstrap 可以为来自极端偏斜分布的样本提供足够准确的置信限,但双 bootstrap 可以显着提高两个尾部的置信覆盖率。 Nankervis有一些不错的结果,并提供了出色的计算算法。但是我找不到任何软件可以实现这一点。
R代码说明1.以上:
## Exact CI for median from DescTools package SignTest.default
## See also ttp://www.stat.umn.edu/geyer/old03/5102/notes/rank.pdf,
## http://de.scribd.com/doc/75941305/Confidence-Interval-for-Median-Based-on-Sign-Test
cimed <- function(x, alpha=0.05, na.rm=FALSE) {
if(na.rm) x <- x[! is.na(x)]
n <- length(x)
k <- qbinom(p=alpha / 2, size=n, prob=0.5, lower.tail=TRUE)
## Actual CL: 1 - 2 * pbinom(k - 1, size=n, prob=0.5) >= 1 - alpha
sort(x)[c(k, n - k + 1)]
}
n <- 20
m <- 20000
cil <- cilt <- 0
z <- qt(0.975, n - 1)
for(i in 1 : m) {
x <- rnorm(n)
cil <- cil + diff(cimed(x))
cilt <- cilt + 2 * z * sqrt(var(x) / n)
}
cil <- cil / m
cilt <- cilt / m
c(cil, cilt, cilt / cil, cil / cilt)