## Setup
x <- 31
n <- 39
a <- 1  ## Prior alpha
b <- 1  ## Prior beta

## Posterior density
theta <- seq(0, 1, 0.005)
plot(theta, dbeta(theta, x+a, n-x+b), type="l")

## Posterior tail probability
pbeta(0.6, x+a, n-x+b)

## Central interval
qbeta(0.025, x+a, n-x+b)
qbeta(0.975, x+a, n-x+b)
## Or:
qbeta(c(0.025, 0.975), x+a, n-x+b)

## function outline
binom.bayes <- function(x, n, a=1, b=1, level=0.95, plot=FALSE, ...) {  
}

## function v1
binom.bayes <- function(x, n, a=1, b=1, level=0.95, plot=FALSE, ...) {
  theta <- seq(0, 1, 0.005)
  plot(theta, dbeta(theta, x+a, n-x+b), type="l")
  qbeta(c(0.025, 0.975), x+a, n-x+b)  
}

## function v2
binom.bayes <- function(x, n, a=1, b=1, level=0.95, plot=FALSE, ...) {
  if (plot) {
    theta <- seq(0, 1, 0.005)
    plot(theta, dbeta(theta, x+a, n-x+b), type="l")    
  }
  qbeta(c(0.025, 0.975), x+a, n-x+b)  
}

## test drive
binom.bayes(31,39)
binom.bayes(31,39, plot=TRUE)

## dots not working
binom.bayes(31,39, plot=TRUE, col="blue")

## function v3
binom.bayes <- function(x, n, a=1, b=1, level=0.95, plot=FALSE, ...) {
  if (plot) {
    theta <- seq(0, 1, 0.005)
    plot(theta, dbeta(theta, x+a, n-x+b), type="l", ...)
  }
  qbeta(c(0.025, 0.975), x+a, n-x+b)  
}
binom.bayes(31,39, plot=TRUE, col="blue", lwd=3, ylab="Posterior density")

## function v4
binom.bayes <- function(x, n, a=1, b=1, level=0.95, plot=FALSE, ...) {
  if (plot) {
    theta <- seq(0, 1, 0.005)
    plot(theta, dbeta(theta, x+a, n-x+b), type="l", ...)
  }
  lp <- (1-level)/2
  up <- 1-(1-level)/2
  qbeta(c(lp, up), x+a, n-x+b)
}
binom.bayes(31, 39, level=0.9)
binom.bayes(31, 39, level=0.8)

## function v5
binom.bayes <- function(x, n, a=1, b=1, level=0.95, plot=FALSE, ...) {
  if (plot) {
    theta <- seq(0, 1, 0.005)
    plot(theta, dbeta(theta, x+a, n-x+b), type="l", ...)
  }
  lp <- (1-level)/2
  up <- 1-(1-level)/2
  interval <- qbeta(c(lp, up), x+a, n-x+b)
  mode <- (x+a-1) / (x+a+n-x+b-2)
  list(mode=mode, interval=interval)
}
binom.bayes(31, 39)

## monte carlo
theta1 <- rbeta(10000, 32, 9)
theta2 <- rbeta(10000, 11, 11)
mean(theta1 > theta2)