[R] how to compute highest density interval?

[Duncan Murdoch]

On 24/11/2007 7:43 AM, gallon li wrote:
> Suppose i want to compute a 95% highest density for a beta distribution
> beta(a,b)
> 
> the two end points x1 and x2 shoudl satisfy the following two equations:
> 
> pbeta(x1,a,b)-pbeta(x2,a,b)=95%
> 
> dbeta(x1,a,b)=dbeta(x2,a,b)
> 
> Is there any fast way to compute x1 and x2 in R?

I don't know if it's "fast", but uniroot can do it:

densitydiff <- function(lower, a, b, level=0.95) {
   plower <- pbeta(lower, a, b)
   pupper <- plower + level
   upper <- qbeta(pupper, a, b)
   return(dbeta(lower, a, b) - dbeta(upper, a, b))
}

a <- 2
b <- 10
x2 <- uniroot(densitydiff, c(0, qbeta(0.05, a,b)), a=a, b=b)$root
(and then calculate the upper limit the same way densitydiff did).

Duncan Murdoch