N <- 50
   set.seed(42)
   Y <- sample(c(0,1),N,pr=c(.2,.8),replace=T)
   postbetbin <- function(p, Y, N, alpha, beta) { 
       return(dbinom(sum(Y),N,p)*dbeta(p,alpha,beta))
   } 
   lbinom <- function(p,Y,N) dbinom(Y,N,p) 
   dbeta2 <- function(ab,p) return(unlist(lapply(p,dbeta,shape1=ab[1],shape2=ab[2]))) 
   lines2 <- function(y,x,...) lines(x,y[-1],lty=y[1],...)
  
   x=seq(0,1,l=200)
   alphabeta=matrix(0,nrow=4,ncol=2) 
   alphabeta[1,]=c(1,1) 
   alphabeta[2,]=c(60,60)
   alphabeta[3,]=c(5,5)
   alphabeta[4,]=c(2,5)
   labs <- c("beta(1,1)", "beta(60,60)", "beta(5,5)", "beta(2,5)")
   priors <- apply(alphabeta, 1, dbeta2, p=x)
   
   par(mfrow=c(2,2),lwd=2,mar=rep(3,4) ,cex.axis=.6)
   for(j in 1:4){
     plot(x,unlist(lapply(x,lbinom,Y=sum(Y),N=N)), type="l",
        xlab="p", col="gray",ylab="", main=paste("Prior is",labs[j]),
        ylim=c(0,.3))
     lines(x, unlist(lapply(x,postbetbin,Y=sum(Y),N=N,
         alpha=alphabeta[j,1],beta=alphabeta[j,2])),lty=1)
     par(new=T)
     plot(x,dbeta(x,alphabeta[j,1],alphabeta[j,2]),lty=3,axes=F,
       type='l', xlab="", ylab="",ylim=c(0,9))
       axis(4)
     legend(0, 9, legend=c("Prior", "Likelihood", "Posterior"), 
          lty=c(3, 1, 1), col=c( "black", "gray", "black"),cex=.6)
     mtext("Prior", side=4, outer=F, line=2, cex=.6)
     mtext("Likelihood/posterior", side=2, outer=F, line=2, cex=.6)
   }