Credible set for beta distribution How can I approximate beta distribution $\alpha=x+1$ and $\beta = 14-x$ to normal distribution? Or, could you please tell me how to calculate HPD credible set for beta distribution?
 A: A level $\alpha$ "Highest Posterior Density" (HPD) interval for a (posterior) distribution $F$ with (continuous) density $f$ and mode $\mu$ is an interval $I=[x,y]$ containing $\mu$ for which


*

*$1-\alpha$ of the probability is in the interval: $F(I) = F(y) - F(x) = 1-\alpha$.

*The densities are the same at either end: $f(x) = f(y)$.
Among the various strategies to find $I$, one that stands out as generally effective is the following.


*

*Choose a reasonable starting value for $x$, such as $F^{-1}(\alpha/2)$.

*Define the "$\alpha$ offset" of $x$ to be the point $y$ for which $[x,y]$ has probability $1-\alpha$. Thus
$$y = F^{-1}(F(x) + 1 - \alpha)$$
provided $F(x) \le \alpha$.

*Search for $x$ in the interval $(-\infty, F^{-1}(\alpha))$ at which $f(x) = f(y)$.  The unimodality of $F$ and the continuity of $f$ guarantee such an $x$ exists and is unique.
The search (3) can be carried out in practice by minimizing $(f(y)-f(x))^2$ plus a penalty term in case the probability of $[x,y]$ is not exactly $1-\alpha$.  (The penalty term is useful in case the search procedure provides a candidate value of $x$ for which $F(x)$ exceeds $\alpha$, in which case a valid offset $y$ cannot be found.)
Applying such a general procedure would be a better idea for a Beta distribution compared to using a Normal approximation, because Betas tend to be skewed (unless their parameters $a$ and $b$ are relatively similar).
For example, the orange region in the figure covers a $1-0.05$ HPD interval for a Beta$(11,4)$ distribution whose density is graphed.  The dashed gray line shows the common value of the density at the endpoints.

Here is the R code that performed the calculation.  It is written to be very general: if you can supply functions to compute $F$, $f$, and $F^{-1}$, it will work.  (An example for Normal distributions has been commented out.)
offset <- function(x, alpha=0.05, F=pbeta, F.Inv=qbeta, ...) {
  q <- F(x, ...)
  y <- F.Inv(min(q + 1-alpha, 1), ...)
}
f <- function(x, alpha=0.05, F=pbeta, F.Inv=qbeta, f.dist=dbeta, ...) {
  y <- offset(x, alpha, F=F, F.Inv=F.Inv, ...)
  mapply(function(u,v) diff(f.dist(c(u,v), ...))^2, x, y) +
    (diff(F(c(x,y), ...)) - (1-alpha))^2
}
#
# Specify the problem.
#
alpha <- 0.05    # Level of the CI (between 0 and 1)
a <- 11; b <- 4  # Parameters
F <- pbeta       # Distribution function
F.Inv <- qbeta   # Inverse distribution function
f.dist <- dbeta  # Density function
x.min <- 0       # Minimum to consider
x.max <- 1       # Maximum to consider
# F <- pnorm
# F.Inv <- qnorm
# f.dist <- dnorm
# x.min <- -2
# x.max <- 5
#
# Find the solution.
#
x.0 <- qbeta(alpha/2, a, b)
x.lim <- qbeta(alpha, a, b)
sol <- optimize(function(x) f(x, alpha, F=F, F.Inv=F.Inv, f.dist=f.dist, a, b), 
                interval=c(x.min, x.lim), tol=1e-10)
x <- sol$minimum
y <- offset(x, alpha, F=F, F.Inv=F.Inv, a, b)
#
# Plot the solution.
#
u <- seq(x, y, length.out=1001)
v <- f.dist(u, a, b)
u <- c(u[1], u, u[length(u)])
v <- c(0, v, 0)
curve(f.dist(x, a, b), xlim=c(x.min, x.max), ylim=c(0,max(v)*1.2), n=1001, xlab="X", ylab="Density")
polygon(u, v, col="#f8e0a0", border=NA)
abline(h = f.dist(x, a, b), col="Gray", lty=2, lwd=2)
curve(f.dist(x, a, b), add=TRUE, lwd=2)

