Monte Carlo part 2

rm(list = ls())
p<-0.15 ; n <- 20; no.rep <- 1000
asym_lt95 <- rep(NA,no.rep)
asym_rt95 <- rep(NA,no.rep)
asym_lt99 <- rep(NA,no.rep)
asym_rt99 <- rep(NA,no.rep)
for(i in 1:no.rep){
  # print(i)
  set.seed(i)
  numSuccess <- rbinom(1,n,p)
  phat = numSuccess / n
  asym_lt95[i] <- phat-qnorm(0.975)*sqrt(phat*(1-phat)/n)
  asym_rt95[i] <- phat+qnorm(0.975)*sqrt(phat*(1-phat)/n)
  asym_lt99[i] <- phat-qnorm(0.995)*sqrt(phat*(1-phat)/n)
  asym_rt99[i] <- phat+qnorm(0.995)*sqrt(phat*(1-phat)/n)
}
asym_lt95 = pmax(asym_lt95, 0)
asym_rt95 = pmin(asym_rt95, 1)
asym_lt99 = pmax(asym_lt99, 0)
asym_rt99 = pmin(asym_rt99, 1)
asym_95_ci_cov_rate = mean((asym_lt95<=p) & (p<=asym_rt95)) # 0.803
asym_99_ci_cov_rate = mean((asym_lt99<=p) & (p<=asym_rt99)) # 0.961

asym_95_ci_len = mean(asym_rt95 - asym_lt95) # 0.2766
asym_99_ci_len = mean(asym_rt99 - asym_lt99) # 0.3373

## pretty code
rm(list = ls())
library(plyr)
library(magrittr)

mc_binomial_f = function(no.rep, n, p, prob = 0.95){
  dat = laply(1:no.rep, function(i){
    set.seed(i)
    rbinom(1,n,p) / n
  })
  sd_dat = sqrt(dat*(1-dat)/n)
  dat %>% {cbind(pmax(0, . - qnorm(prob/2+0.5)*sd_dat),
    pmin(. + qnorm(prob/2+0.5)*sd_dat, 1))}
}
p = 0.15
ci_95 = mc_binomial_f(1000, 20, p)
ci_99 = mc_binomial_f(1000, 20, p, 0.99)

ci_95 %>% {mean((.[,1]<=p) & (p<=.[,2]))}
ci_99 %>% {mean((.[,1]<=p) & (p<=.[,2]))}

ci_95 %>% {mean((.[,2] - .[,1]))}
ci_99 %>% {mean((.[,2] - .[,1]))}