4
$\begingroup$

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?

$\endgroup$
5
$\begingroup$

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.

Figure

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)
$\endgroup$
  • $\begingroup$ 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 :( $\endgroup$ – Niniko Jul 7 '15 at 19:33
  • 1
    $\begingroup$ I cannot make sense of that request, because if you are integrating a probability density you cannot obtain a value any greater than $1$. $\endgroup$ – whuber Jul 7 '15 at 20:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.