beta_binom<-function(n,y,a=1,b=1,main=""){
    #likelihood: y|theta~binom(n,theta)
    #prior: theta~beta(a,b)
    #posterior: theta|y~beta(a+y,n-y+b)

    theta<-seq(0.001,0.999,0.001)
    prior<-dbeta(theta,a,b)
    if(n>0){likelihood<-dbinom(rep(y,length(theta)),n,theta)}
    if(n>0){posterior<-dbeta(theta,a+y,n-y+b)}

    #standardize!
    prior<-prior/sum(prior)
    if(n>0){likelihood<-likelihood/sum(likelihood)}
    if(n>0){posterior<-posterior/sum(posterior)}

    ylim<-c(0,max(prior))
    if(n>0){ylim<-c(0,max(c(prior,likelihood,posterior)))}

    plot(theta,prior,type="l",lty=2,xlab="theta",ylab="",main=main,ylim=ylim)
    if(n>0){lines(theta,likelihood,lty=3)}
    if(n>0){lines(theta,posterior,lty=1,lwd=2)}
    legend("topright",c("prior","likelihood","posterior"),
           lty=c(2,3,1),lwd=c(1,1,2),inset=0.01,cex=.5)
}
 

par(mfrow=c(2,2))
beta_binom(3,2,2.5,7.5,main="Prior: beta(2.5,7.5), data: 2/3")
beta_binom(3,2,25,75,main="Prior: beta(25,75), data: 2/3")
beta_binom(30,20,2.5,7.5,main="Prior: beta(2.5,7.8), data: 20/30")
beta_binom(300,200,2.5,7.5,main="Prior: beta(2.5,7.5), data: 200/300")