sdchasalow

n <- 50
pvals <- 0.5
nosim <- 100
coverage <- lapply(pvals, function(p) {
phats <- (rbinom(nosim, prob = p, size = n) + 2)/(n + 4)
ll <- phats - qnorm(0.975) * sqrt(phats * (1 - phats)/n)
ul <- phats + qnorm(0.975) * sqrt(phats * (1 - phats)/n)
cbind(lower = ll, upper = ul)
})