# qbinom.R # T. J. Finney 2017-03-20 message("Check confidence interval limits for binomial distribution...") # Reset. rm(list=ls(all=TRUE)) # Set parameters. setwd("P:/R/scripts") message("Lower quantile: ", p1 <- 0.025) message("Upper quantile: ", p2 <- 0.975) message("No. of trials: ", size <- 26) message("Prob. of success: ", prob <- 0.25) # Set N (a large number) N <- 10000 # Generate a large number of Bernoulli trials (counting numbers of successes) values <- lapply(1:N, function(x) { sum(rbinom(size, 1, prob)) }) # Compare quantiles of generated data with theoretical values values <- unlist(values) message("Table of values for ", N, " iterations") print(table(values)) message("Confidence interval limits (generated data)") print(quantile(values, c(p1, p2))) message("Confidence interval limits (theory)") theory <- qbinom(c(p1, p2), size, prob) print(theory) message("Corresponding distances") print(round(theory/size, digits=3))