Rのカバレッジ確率（forループ）

Question

以下のように計算した信頼区間Iのカバレッジ確率を計算しようとしています。私のコードは動作しません。誰でも助けてくれますか？

p1 <- 0.01 
p2 <- 0.05 
p3 <- 0.1 
p4 <- 0.3 
p5 <- 0.5 
p6 <- 0.7 
p7 <- 0.9 
p8 <- 0.95 
p9 <- 0.99 
p10 <- 0.9999 

N = 1000

for (i in 1:10000){ 
    set.seed (123467) 
    l1 <- (mean(rbinom(1000,1,p1))) - 1.96 * ((mean(rbinom(1000,1,p1))* (1 -mean(rbinom(1000,1,p1))))/1000)^(1/2) 
    l2 <- (mean(rbinom(1000,1,p2))) - 1.96 * ((mean(rbinom(1000,1,p2))* (1 -mean(rbinom(1000,1,p2))))/1000)^(1/2) 
    l3 <- (mean(rbinom(1000,1,p3))) - 1.96 * ((mean(rbinom(1000,1,p3))* (1 -mean(rbinom(1000,1,p3))))/1000)^(1/2) 
    l4 <- (mean(rbinom(1000,1,p4))) - 1.96 * ((mean(rbinom(1000,1,p4))* (1 -mean(rbinom(1000,1,p4))))/1000)^(1/2) 
    l5 <- (mean(rbinom(1000,1,p5))) - 1.96 * ((mean(rbinom(1000,1,p5))* (1 -mean(rbinom(1000,1,p5))))/1000)^(1/2) 
    l6 <- (mean(rbinom(1000,1,p6))) - 1.96 * ((mean(rbinom(1000,1,p6))* (1 -mean(rbinom(1000,1,p6))))/1000)^(1/2) 
    l7 <- (mean(rbinom(1000,1,p7))) - 1.96 * ((mean(rbinom(1000,1,p7))* (1 -mean(rbinom(1000,1,p7))))/1000)^(1/2) 
    l8 <- (mean(rbinom(1000,1,p8))) - 1.96 * ((mean(rbinom(1000,1,p8))* (1 -mean(rbinom(1000,1,p8))))/1000)^(1/2) 
    l9 <- (mean(rbinom(1000,1,p9))) - 1.96 * ((mean(rbinom(1000,1,p9))* (1 -mean(rbinom(1000,1,p9))))/1000)^(1/2) 
    l10 <- (mean(rbinom(1000,1,p10))) - 1.96 * ((mean(rbinom(1000,1,p10))* (1 -mean(rbinom(1000,1,p10))))/1000)^(1/2) 

    u1 <- (mean(rbinom(1000,1,p1))) + 1.96 * ((mean(rbinom(1000,1,p1))* (1 -mean(rbinom(1000,1,p1))))/1000)^(1/2) 
    u2 <- (mean(rbinom(1000,1,p2))) + 1.96 * ((mean(rbinom(1000,1,p2))* (1 -mean(rbinom(1000,1,p2))))/1000)^(1/2) 
    u3 <- (mean(rbinom(1000,1,p3))) + 1.96 * ((mean(rbinom(1000,1,p3))* (1 -mean(rbinom(1000,1,p3))))/1000)^(1/2) 
    u4 <- (mean(rbinom(1000,1,p4))) + 1.96 * ((mean(rbinom(1000,1,p4))* (1 -mean(rbinom(1000,1,p4))))/1000)^(1/2) 
    u5 <- (mean(rbinom(1000,1,p5))) + 1.96 * ((mean(rbinom(1000,1,p5))* (1 -mean(rbinom(1000,1,p5))))/1000)^(1/2) 
    u6 <- (mean(rbinom(1000,1,p6))) + 1.96 * ((mean(rbinom(1000,1,p6))* (1 -mean(rbinom(1000,1,p6))))/1000)^(1/2) 
    u7 <- (mean(rbinom(1000,1,p7))) + 1.96 * ((mean(rbinom(1000,1,p7))* (1 -mean(rbinom(1000,1,p7))))/1000)^(1/2) 
    u8 <- (mean(rbinom(1000,1,p8))) + 1.96 * ((mean(rbinom(1000,1,p8))* (1 -mean(rbinom(1000,1,p8))))/1000)^(1/2) 
    u9 <- (mean(rbinom(1000,1,p9))) + 1.96 * ((mean(rbinom(1000,1,p9))* (1 -mean(rbinom(1000,1,p9))))/1000)^(1/2) 
    u10 <- (mean(rbinom(1000,1,p10))) + 1.96 * ((mean(rbinom(1000,1,p10))* (1 -mean(rbinom(1000,1,p10))))/1000)^(1/2) 

} 



CI1000 <- matrix(c(l1,l2,l3,l4,l5,l6,l7,l8,l9,l10, u1,u2,u3,u4,u5,u6,u7,u8,u9,u10,u1-l1, u2-l2, u3-l3, u4-l4, u5-l5, u6-l6, u7-l7, u8-l8, u9-l9, u10-l10),ncol=3,nrow=10,byrow=F) 
rownames(CI1000) <- c("p=0,01","p=0.05","p=0.1","p=0,3","p=0,5","p=0,7", "p=0,9","p=0.95","p=0.99", "p=0.9999") 
colnames(CI1000) <- c("lower bound","upper bound", "width") 

カバレッジ確率

for (i in 1:10000){ 
    set.seed (123467) 
    badl1 <- sum(ifelse(p1 < (mean(rbinom(1000,1,p1))) - 1.96 * ((mean(rbinom(1000,1,p1))* (1 -mean(rbinom(1000,1,p1))))/1000)^(1/2), 0,1)) 
    badu1 <- sum(ifelse(p1 > (mean(rbinom(1000,1,p1))) + 1.96 * ((mean(rbinom(1000,1,p1))* (1 -mean(rbinom(1000,1,p1))))/1000)^(1/2), 0,1)) 

} 
(badl1 + badu1) /1000 -> bad.frac1 

の信頼区間THANKS！

Answer 1

これははるかに簡単で正しいです。 （widthの必要はありませんが、私はそれを保持している。理由は分からないのです。）

ci <- function(i){ 
    lo <- m[i] - 1.96 * (m[i]*(1 - m[i])/n)^0.5 
    hi <- m[i] + 1.96 * (m[i]*(1 - m[i])/n)^0.5 
    c(lo, hi) 
} 

p <- c(p1, p2, p3, p4, p5, p6, p7, p8, p9, p10) 

m <- length(p) 
n <- 1000 
Runs <- 10000 
CI <- array(NA, dim = c(m, 3, Runs)) 
set.seed (123467) # do this just once outside the loop 
for (i in 1:Runs){ 
    x <- matrix(rbinom(10*n, 1, rep(p, each = n)), ncol = 10) 
    m <- colMeans(x) 
    CI[, 1:2, i] <- t(sapply(seq_along(m), ci)) 
    CI[, 3, i] <- CI[, 2, i] - CI[, 1, i] 
} 
dimnames(CI) <- list(c("p=0,01","p=0.05","p=0.1","p=0,3","p=0,5","p=0,7", "p=0,9","p=0,95","p=0,99", "p=0,9999"), 
        c("lower.bound","upper.bound", "width"), 
        sprintf("R%05d", 1:Runs)) 

CI[,, 200] # example result of loop 

bad <- array(NA, dim = c(m, Runs)) 
for(i in 1:Runs){ 
    bad[, i] <- as.integer(!(CI[, 1, i] < p & p < CI[, 2, i])) 
} 

prop_bad <- rowMeans(bad) 
names(prop_bad) <- c("p=0,01","p=0.05","p=0.1","p=0,3","p=0,5","p=0,7", "p=0,9","p=0,95","p=0,99", "p=0,9999") 
prop_bad 
# p=0,01 p=0.05 p=0.1 p=0,3 p=0,5 p=0,7 p=0,9 p=0,95 
# 0.0742 0.0582 0.0460 0.0498 0.0508 0.0535 0.0436 0.0582 
# p=0,99 p=0,9999 
# 0.0742 0.9027 

注標準間隔が十分に持っているので、最後の結果は、0.9027に等しいbadの割合が、異常ではないこと既知のカバレッジの問題は、真の割合で、端の近くで0と1です。