CI-binomial-simulation.R

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)
})

png("binomialCI-simulation.png")
plot( 0:1, c(0, 101), type = 'n', xlab = "95% CI for p", ylab = "Replicate Number")
abline(v = 0.5, lty = 2)
segments(coverage[[1]][,1], 1:100, coverage[[1]][, 2], 1:100)
title("Simulation of 100 coin flips with true p = 0.5")
dev.off()