library(VGAM)

# RETURNS: width of 95% posterior credible interval for θ given a beta(α, β)
# prior, with data consisting of n trials with k successes.
#
posterior.width <- function(α, β, n, k) {
 x1 <- qbeta(.975, α+k, β+n-k)
 x2 <- qbeta(.025, α+k, β+n-k)
 x1 - x2
}

# RETURNS: Expected value of width of posterior distribution for θ.
# Expectation is taken over the prior predictive distribution for k, the
# number of successes, assuming n trials and a beta(α, β) prior for θ.
#
expected.posterior.width <- function(α, β, n) {
 sum(dbetabinom.ab(0:n, n, α, β) * posterior.width(α, β, n, 0:n))
}

plot.expected.posterior.width <- function(α, β, n.max) {
 n.vec <- unique(1 + round((n.max - 1) * (0:300)/300))
 y <- vapply(n.vec, function(n)expected.posterior.width(α, β, n), 0)
 plot(n.vec, y, ylab=expression(g(w)), xaxs='i', yaxs='i', xlab=expression(n),
 typ='l', lwd=2, ylim=c(0,y[1]))
}