## binomial beta model
## ADM 6/12/2007
library(MCMCpack)

## data
y <- 3
n <- 5

## priors
a <- 3
b <- 2

## posterior mean and standard deviation
cat("posterior mean", (y+a) / (a + n + b), "\n")
cat("posterior standard deviation", sqrt((y+a) * (n - y + b) / ((a + n + b)^2 *
   (n + a + b + 1))), "\n\n")

## compute Pr(pi <= .25 | y)
hold <- pbeta(.25, y+a, n-y+b)
cat("Pr(pi <= .25 | y): ", hold, "\n")

## compute Pr(.4 < pi <= .6 | y)
hold <- pbeta(.6, y+a, n-y+b) - pbeta(.4, y+a, n-y+b)
cat("Pr(.4 < pi <= .6 | y): ", hold, "\n")

## compute central 95% credible interval
left <- qbeta(0.025, y+a, n-y+b)
right <- qbeta(0.975, y+a, n-y+b) 
cat("central 95% credible interval: ", left, ", ", right, "\n\n", sep="")

## perform inference using a Monte Carlo simulation
post.sample <- MCbinomialbeta(y, n, a, b, mc=1000)
print(MCbinomialbeta)
print(str(post.sample))
print(summary(post.sample))
plot(post.sample)