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.
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
