binomialCI-simulation.png

CI-binomial-simulation.R

# Fork of Scott's coin flip simulation
# https://gist.github.com/maptracker/e34a8596a9c2a7f261d6b892e5df42c5

flipCoin <- function( n = 50, pvals = 0.5, nosim = 100, file = NULL ) {
    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)
    })[[1]] # De-listify

    if (!is.null(file)) png(file)
    plot(0:1, c(0, nosim + 1), type = 'n',
         xlab = "95% CI for p", ylab = "Replicate Number")
    abline(v = pvals, lty = 2, col = "blue")
    ## Highlight ranges that do not cover the true mean
    color <- apply(coverage, 1, function (x) {
        ifelse(x[1] <= pvals && pvals <= x[2], "black", "red")
    })

    segments(coverage[,1], 1:nosim, coverage[, 2], 1:nosim, col = color)
    ## Draw a density plot (bar) of overlapping coverages
    alpha <- as.integer(255 / nosim)
    if (alpha < 1) alpha <- 1
    rect(coverage[,1], 0, coverage[, 2], nosim/-20,
         col = rgb(0,0,0,alpha, maxColorValue = 255), border = NA)
    title(sprintf("Simulating %s coin flips, %d replicates\nTrue p = %s",
                  n, nosim, pvals))
    if (!is.null(file)) dev.off()
    coverage
}
c <- flipCoin()

Please give me a pull request. Thanks for the corrections and enhancements.