theta_seq(0,1,length=5000)
par(mfrow=c(5,5))   #examine beta densities
for(i in c(1,5,10,20, 100)) {
  for(j in c(1,5,10,20,100))
    {plot(theta,dbeta(theta,i,j), type="l")
     title(paste("Beta(",as.character(i),",",as.character(j),")", sep=""))}
}

plot(theta, dbeta(theta,.5,.5), type="l")
par(mfrow=c(2,3))   #sample posteriors uniform prior
for(i in 1:6) {
  y_rbinom(1,10*i,.2)
  plot(theta,dbeta(theta,y+1,10*i-y+1), type="l")
  title(paste("theta=0.2, n =", as.character(10*i), "y =", as.character(y)))
}
#repeat

par(mfrow=c(1,1))   #quantiles of the beta
plot(theta,dbeta(theta,2,10), type="l")
abline(v=qbeta(0.5,2,10),col=2)
abline(v=qbeta(0.025,2,10),col=2)
abline(v=qbeta(0.975,2,10),col=2)


# Example
c(306/651,300/705)  #MLEs

theta <- seq(0,1, length=10000)
plot(theta,dbeta(theta,301, 406 ),type="l",col=1, ylab="Posterior Density")
lines(theta,dbeta(theta,307,346),col=2)

d <- rbeta(5000,301,406)-rbeta(5000,307,346)
hist(d,nclass=30,prob=T)
summary(d)
ind <- (d>0)
mean(d>0)
round(quantile(d,prob=c(0.05,0.95)), 4)
round(quantile(d,prob=c(0.025,0.975)), 4)

quantile(d,prob=c(0.05,0.95))