# Calculate confidence interval for given estimate and size using plug in method
calc_ci = function (phat, n, level = .95) {
  z = qnorm((1-level)/2, lower.tail = FALSE)
  lower = phat - z * sqrt(phat * (1-phat) / n)
  upper = phat + z * sqrt(phat * (1-phat) / n)
  c(lower, upper)
}

# Determine if true value p given falls in CI
in_ci = function(ci, p) {
  ci[1] <= p && p <= ci[2]
}

# Collects data for given m, p and the coverage of the samples
stats <- data.frame(matrix(ncol = 3, nrow = 0))
cnames <- c("m", "p", "coverage")
colnames(stats) <- cnames

# Iterate over the target p and m values
for (p in c(.5, .75, .9)) {
  for (m in c(10, 20, 50, 100, 250, 500)) {
    # Track the number of times p fals in the CI
    num.captured = 0
    iterations = 1000
    for (i in 1:iterations) {
      data.sample = rbinom(m, 1, p)
      phat = mean(data.sample)
      if (in_ci(calc_ci(phat, m), p)) {
        num.captured = num.captured + 1
      }
    }
    stats[nrow(stats) + 1,] = list(m, p, num.captured/iterations)
  }
}