Доверительный интервал долей Опции

[плав]

Поскольку в другой теме было много споров о разных ДИ для долей и огрномный список цитат, решил - для иллюстрации - провести вычислительный эксперимент.
Итак. Были смоделированы популяции в которой содержится х объектов одного класса и 1-х объектов другого класса (х менялась от 1 до 10%). Из этой популяции брались случайные выборки размером 40 объектов. Оценивалось количество объектов одного и другого класса в выборке и рассчитывались доверительные интервалы по Клопперу-Пирсону, Агрести-Коулу и по распределению Пуассона (значения менее 0 заменялись на нулевые).
Таких выборок бралось 10 000 и затем рассчитывался вероятность покрытия популяционного значения доверительным интервалом и средняя ширина доверительного интервала. Что в результате (это небольшой размер выборки и малая вероятнсть:
Ширина 95%ДИ Покрытие
pi__ КП__ АК__ Пу__ КП__ АК__ Пу__
1% 0,104 0,118 0,110 0,993 0,993 0,993
2% 0,119 0,130 0,126 0,992 0,951 0,992
3% 0,132 0,141 0,141 0,994 0,970 0,994
4% 0,146 0,152 0,156 0,979 0,979 0,979
5% 0,157 0,162 0,169 0,986 0,952 0,986
6% 0,168 0,170 0,181 0,991 0,970 0,991
7% 0,177 0,178 0,192 0,981 0,981 0,981
8% 0,186 0,185 0,203 0,988 0,965 0,988
9% 0,195 0,192 0,214 0,953 0,974 0,970
10% 0,203 0,198 0,223 0,972 0,962 0,972
При малых значениях популяционной вероятности (менее 8%) интервал Клоппера-Пирсона является более узким, при больших - боле узкий интервал Агрести-Коула.
Покрытие прыгает, почти всюду больше номинального уровня 95%, однако среднее покрытие для КП - 98,3%, для Агрести - 97,0% и для Пуассона - 98,5%. Агрести-Коула немного ближе к номинальному уровню.
Однако если смотреть на данные реально, принципиальных различий между этими тремя методами нет. В большинстве случаев они дают одинаковые результаты и, как и показано в других работах, КП немного более консервативен, а АК немного более широкий при малых значениях популяционной вероятности.

[nokh]

Благодарю за готовый код! Теперь когда понадобится - буду использовать. Именно BCa - не критично, просто из модного бутстрепа он более моден . Появится время - сравню ДИ по бутстрепу с ДИ по регресии.

[p2004r]

Благодарю за готовый код! Теперь когда понадобится - буду использовать. Именно BCa - не критично, просто из модного бутстрепа он более моден . Появится время - сравню ДИ по бутстрепу с ДИ по регресии.

оно оказывается вот как выглядит

Код

> bmem.ci.bc
function (par.boot, par0, cl = 0.95)
{
    se.boot <- apply(par.boot, 2, sd, na.rm = TRUE)
    estimate <- par0
    p <- ncol(par.boot)
    ci <- NULL
    for (i in 1:p) {
        ci <- rbind(ci, bmem.ci.bc1(par.boot[, i], par0[i], cl))
    }
    cbind(estimate, se.boot, ci)
}
<environment: namespace:bmem>

> bmem.ci.bc1
function (x, b, cl = 0.95)
{
    n <- length(x)
    z0 <- qnorm(sum(x < b, na.rm = TRUE)/n)
    alpha <- (1 - cl)/2
    alpha <- c(alpha, 1 - alpha)
    alpha <- sort(alpha)
    alpha1 <- alpha
    alpha <- pnorm(2 * z0 + qnorm(alpha))
    dig <- max(2L, getOption("digits"))
    np <- length(alpha)
    qs <- quantile(x, alpha, na.rm = TRUE)
    names(qs) <- paste(if (np < 100)
        formatC(100 * alpha1, format = "fg", width = 1, digits = dig)
    else format(100 * alpha1, trim = TRUE, digits = dig), "%",
        sep = "")
    qs
}
<environment: namespace:bmem>

они имели в виду под результатом бутстрепа что то не из boot

вот такое выходит после загрузки library(bmem)

Код

>  bmem.ci.norm(as.data.frame(result),mean(c(rep(2,10),rep(5,3)*2,rep(10,2)*4)),cl=.95)
       estimate  se.boot     2.5%    97.5%
result 8.666667 1.315171 6.088979 11.24435

>  bmem.ci.bc(as.data.frame(result),mean(c(rep(2,10),rep(5,3)*2,rep(10,2)*4)),cl=.95)
       estimate  se.boot 2.5%    97.5%
result 8.666667 1.315171  6.2 11.33333