plot.beta <- function(h, t) {
  support <- (0:100)/100
  density <- dbeta(support, h+1, t+1)
  plot(support, density, ylim=c(0,8), type="l")
  Sys.sleep(ifelse(h+t==0, 1, 1/(h+t)))
}
run.model <- function(trials, p) {
  heads <- c(0,cumsum(rbinom(trials, 1, p)))
  mapply(plot.beta, heads, (0:trials)-heads)
}
ignore <- run.model(100, 0.6)