[R] how to compute highest density interval?

[Ben Bolker]

Duncan Murdoch-2 wrote:
> 
> 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
> 


The tcredint function in the emdbook package (on CRAN) incorporates
a similar approach.

  Ben Bolker

-- 
View this message in context: http://www.nabble.com/how-to-compute-highest-density-interval--tf4865715.html#a13930376
Sent from the R help mailing list archive at Nabble.com.