gbinom <- function(n, p, low=0, high=n,scale = F) { sd <- sqrt(n * p * (1 - p)) if(scale && (n > 10)) { low <- max(0, round(n * p - 4 * sd)) high <- min(n, round(n * p + 4 * sd)) } values <- low:high probs <- dbinom(values, n, p) plot(c(low,high), c(0,max(probs)), type = "n", xlab = "Possible Values", ylab = "Probability", main = paste("Binomial Distribution \n", "n =", n, ", p =", p)) lines(values, probs, type = "h", col = 2) abline(h=0,col=3) return(invisible()) }