[R] how to compute highest density interval?
Ben Bolker
bolker at ufl.edu
Sat Nov 24 23:56:06 CET 2007
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.
More information about the R-help
mailing list