# Coverage simulation
N <- 10000
n <- 10
pi <- 0.25
x <- rbinom(N, prob=0.25, size=n)
covered <- numeric(N)
for (i in 1:N) {
  ci <- binom.test(x[i], n)$conf
  covered[i] <- (ci[1] < pi) & (pi < ci[2])
}
mean(covered)

# Coverage simulation loop
N <- 10000
n <- 20
pi <- c(0.01, 0.02, 0.04, 0.08, 0.16, 0.32, 0.5, 0.68, 0.84, 0.92, 0.96, 0.98, 0.99)
coverage <- numeric(length(pi))
pb <- txtProgressBar(1, length(pi), style=3)
for (j in 1:length(pi)) {
  x <- rbinom(N, prob=pi[j], size=n)
  covered <- numeric(N)
  for (i in 1:N) {
    ci <- binom.test(x[i], n)$conf
    covered[i] <- (ci[1] < pi[j]) & (pi[j] < ci[2])
  }
  coverage[j] <- mean(covered)
  setTxtProgressBar(pb, j)
}

# Coverage results
plot(pi, coverage, type="l", lwd=3, ylim=c(0.8,1), las=1)
abline(h=0.95, col="red", lwd=3)

# Sourcing ci methods
source("http://myweb.uiowa.edu/pbreheny/5710/f16/labs/binom.bayes.R")
prop.wald <- function(x, n, level=0.95) {
  pi.hat <- x/n
  SE <- sqrt(pi.hat*(1-pi.hat)/n)
  z <- qnorm(c((1-level)/2, 1-(1-level)/2))
  pi.hat + z*SE
}
prop.test(31,39) ## Just to show how prop.test works

# Coverage simulation loop 2
N <- 10000
n <- 20
pi <- c(0.01, 0.02, 0.04, 0.08, 0.16, 0.32, 0.5, 0.68, 0.84, 0.92, 0.96, 0.98, 0.99)
coverage <- matrix(NA, length(pi), 4, dimnames=list(pi, c("CP", "Bayes", "Score", "Wald")))
pb <- txtProgressBar(1, length(pi), style=3)
for (j in 1:length(pi)) {
  x <- rbinom(N, prob=pi[j], size=n)
  covered <- matrix(NA, N, 4)
  for (i in 1:N) {
    ci <- binom.test(x[i], n)$conf
    covered[i,1] <- (ci[1] < pi[j]) & (pi[j] < ci[2])
    ci <- binom.bayes(x[i], n)$ci.hpd
    covered[i,2] <- (ci[1] < pi[j]) & (pi[j] < ci[2])
    ci <- prop.test(x[i], n, correct=FALSE)$conf
    covered[i,3] <- (ci[1] < pi[j]) & (pi[j] < ci[2])
    ci <- prop.wald(x[i], n)
    covered[i,4] <- (ci[1] < pi[j]) & (pi[j] < ci[2])
  }
  coverage[j,] <- apply(covered, 2, mean)
  setTxtProgressBar(pb, j)
}

# Coverage results 2
col <- hcl(seq(15, 375, len = 5), 150, 60)[1:4] ## Setting up some colors
matplot(pi, coverage, type="l", lwd=3, lty=1, col=col)
legend("bottom", legend=colnames(coverage), col=col, lwd=3)
matplot(pi, coverage, type="l", lwd=3, lty=1, col=col, ylim=c(0.8,1), las=1)
legend("bottom", legend=colnames(coverage), col=col, lwd=3)
abline(h=0.95, col="gray", lwd=3)