## Binom sim
N <- 1000
n <- 15
alpha <- .1
pi <- seq(.01,.08,.01)
Covered <- array(NA, dim=c(N, length(pi), 4), dimnames=list(1:N, pi, c("Wald", "Score", "CP", "Bayes")))
inside <- function(ci, theta) {theta > ci[1] & theta < ci[2]}
for (i in 1:N) {
  for (j in 1:length(pi)) {
    y <- rbinom(1, size=n, pi[j])
    p <- y/n
    Covered[i, j, 1] <- inside(p + c(-1,1)*qnorm(1-alpha/2)*sqrt(p*(1-p)/n), pi[j])
    Covered[i, j, 2] <- inside(prop.test(y, n, conf=1-alpha)$conf, pi[j])
    Covered[i, j, 3] <- inside(binom.test(y,n, conf=1-alpha)$conf, pi[j])
    Covered[i, j, 4] <- inside(binom.bayes(y,n, level=1-alpha)$ci.hpd, pi[j])
  }
  displayProgressBar(i,N)
}
coverage <- apply(Covered, 2:3, mean)
matplot(pi, jitter(coverage, amount=.005), type="l", lwd=3, col=pal(4), lty=1, xlab=expression(pi), ylab="Coverage")
toplegend(colnames(coverage), lwd=3, col=pal(4), ncol=4)
abline(h=.9)