plot.beta <- function(h, t) { support <- (0:100)/100 bdensity <- dbeta(support, h+1, t+1) plot(support, bdensity, ylim=c(0,8), type="l") ndensity <- dnorm(support, h/(h+t), sqrt(h*t)/(h+t)^(3/2)) points(support, ndensity, type="l") Sys.sleep(1/(h+t)) } run.model <- function(trials, p) { heads <- cumsum(rbinom(trials, 1, p)) mapply(plot.beta, heads, (1:trials)-heads) } ignore <- run.model(100, 0.6)