par(mfrow=c(1,2))

x <- seq(.01,.99,length=99);  top <- 3.25
plot(x,dbeta(x,.5,.5),type="l",xlab=expression(theta),
  ylab="prior density",ylim=c(0,top))
lines(x,dbeta(x,1,1),lty=2)
lines(x,dbeta(x,2,2),lty=3)
legend(.2,top,cex=1.2,legend =c("Beta(.5,.5)  (Jeffreys prior)",
  "Beta(1,1)  (uniform prior)",
  "Beta(2,2)  (skeptical prior)"),lty=1:3)

x <- seq(0,1,length=100);
plot(x,dbeta(x,13.5,3.5),type="l",lty=1,xlab=expression(theta),
  ylab="posterior density")
lines(x,dbeta(x,14,4),lty=2)
lines(x,dbeta(x,15,5),lty=3)
legend(0.1 ,3.5,cex=1.2,legend =c(
 "Beta(13.5,3.5)", "Beta(14,4)", "Beta(15,5)"),lty=1:3)

a <- 0.025 ; b <- 0.5; d <- 0.975
a1 <- qbeta(a,13.5,3.5); a2 <- qbeta(b,13.5,3.5)
a3 <- qbeta(d,13.5,3.5); a4 <- pbeta(0.6,13.5,3.5); a5 <- 1-a4
b1 <- qbeta(a,14,4); b2 <- qbeta(b,14,4)
b3 <- qbeta(d,14,4); b4 <- pbeta(0.6,14,4); b5 <- 1-b4
e1 <- qbeta(a,15,5); e2 <- qbeta(b,15,5)
e3 <- qbeta(d,15,5); e4 <- pbeta(0.6,15,5); e5 <- 1-e4

cat("q_.025    q_.50     q_.975    P(theta >.6)\n")
prob <- c( a1,a2,a3,a5,b1,b2,b3,b5,e1,e2,e3,e5)
p <- t(matrix(prob,4,3))
print(round(p,3))