While studying bayesian statistics, an interval estimation called HPD interval comes in, and I am trying to find HPD interval of 95% of beta distribution.
I've tried several ways, but finding an algorithm that finds HPD interval was difficult.
I want to use:
p <- seq(0, 1, length = 100)
density <- dbeta(p, 7, 25)
to find HPD interval.
Another question would be I know that we get density value from computing dbeta(p, 7, 25), but is there an inverse function of dbeta? I was trying to get the p with corresponding density value.
Thank you in advance for your help!
Look at the code of TeachingDemos::hpd. You need the inverse cdf (icdf).
function(posterior.icdf, conf=0.95, tol=0.00000001,...){
conf <- min( conf, 1-conf )
f <- function(x,posterior.icdf,conf,...){
posterior.icdf(1-conf+x,...) - posterior.icdf(x,...)
}
out <- optimize(f, c(0,conf), posterior.icdf = posterior.icdf,
conf=conf, tol=tol, ...)
return( c( posterior.icdf(out$minimum,...),
posterior.icdf(1-conf+out$minimum,...) ) )
}
Make a picture, you will see it's easy.