# 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 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. $1-\alpha$ of the probability is in the interval: $F(I) = F(y) - F(x) = 1-\alpha$.

2. 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.

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

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$.

3. 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)  • Thank you very much. Could you please help me once more I have beta distr. with (1,14). and the professor told me in order to find 95% HPD credible set I have to integrate the function from 0 to some value approximately 0.2 which has to give the value 1.924 and I have to do this without computer programms. He told me it cannot be approximated to normal distribution. So I do not understand how can I gain the value 1.924 :( – Niniko Jul 7 '15 at 19:33 • I cannot make sense of that request, because if you are integrating a probability density you cannot obtain a value any greater than$1\$. – whuber Jul 7 '15 at 20:12