## binomial conjugate example
## based on code from KQ
## ADM 9/21/2004

# helper functions
prior <- function(x, a, b){dbeta(x, a, b)}
post  <- function(x, y, n, a, b){dbeta(x, y+a, n-y+b)}
samp <- function(x, y, n){post(x, y, n, 1, 1)}
pi <- seq(0,1,0.01)

## construct some plots
pdf("binomial.pdf", height=6.5, width=6.5)
par(mfrow=c(2,2))

# cell one
n <- 5
y <- 3
a <- 3
b <- 2
plot(pi, post(pi, y, n, a, b), ylab="", type="l", col="black", lwd=2,
   main="n=5, y=3, a=3, b=2")
lines(pi, samp(pi, y, n), col="grey80", lwd=2)
lines(pi, prior(pi, a, b), col="grey50", lwd=2)
legend(0, 2.5, legend=c("prior", "likelihood", "posterior"),
       col=c("grey50", "grey80", "black"), lwd=2, cex=.75)

# cell two
n <- 20
y <- 12
a <- 3
b <- 2
plot(pi, post(pi, y, n, a, b), ylab="", type="l", col="black", lwd=2,
   main="n=20, y=12, a=3, b=2")
lines(pi, samp(pi, y, n), col="grey80", lwd=2)
lines(pi, prior(pi, a, b), col="grey50", lwd=2)
legend(0, 4, legend=c("prior", "likelihood", "posterior"),
       col=c("grey50", "grey80", "black"), lwd=2, cex=.75)

# cell three
n <- 100
y <- 60
a <- 3
b <- 2
plot(pi, post(pi, y, n, a, b), ylab="", type="l", col="black", lwd=2,
   main="n=100, y=60, a=3, b=2")
lines(pi, samp(pi, y, n), col="grey80", lwd=2)
lines(pi, prior(pi, a, b), col="grey50", lwd=2)
legend(0, 8, legend=c("prior", "likelihood", "posterior"),
       col=c("grey50", "grey80", "black"), lwd=2, cex=.75)

# cell four
n <- 500
y <- 300
a <- 3
b <- 2
plot(pi, post(pi, y, n, a, b), ylab="", type="l", col="black", lwd=2,
   main="n=500, y=300, a=3, b=2")
lines(pi, samp(pi, y, n), col="grey80", lwd=2)
lines(pi, prior(pi, a, b), col="grey50", lwd=2)
legend(0, 15, legend=c("prior", "likelihood", "posterior"),
       col=c("grey50", "grey80", "black"), lwd=2, cex=.75)
dev.off()

## other summaries
prob1 <- pbeta(0.25, 3+3, 5-3+2)
prob2 <- pbeta(0.6, 6, 4) - pbeta(0.4, 6, 4)

left <- qbeta(0.025, 6, 4)
right <- qbeta(0.975, 6, 4)