beta HPD
While writing an introductory chapter on Bayesian analysis (in French), I came by the issue of computing an HPD region when the posterior distribution is a Beta B(α,β) distribution… There is no analytic solution and hence I resorted to numerical resolution (provided here for α=117.5, β=115.5):
f=function(p){
# find the symmetric
g=function(x){return(x-p*((1-p)/(1-x))^(115.5/117.5))}
return(uniroot(g,c(.504,.99))$root)}
ff=function(alpha){
# find the coverage
g=function(x){return(x-p*((1-p)/(1-x))^(115.5/117.5))}
return(uniroot(g,c(.011,.49))$root)}
and got the following return:
> ff(.95) [1] 0.4504879 > f(ff(.95)) [1] 0.5580267
which was enough for my simple book illustration… Since (.450,558) is then the HPD region at credible level 0.95.
Xi'an's Og 
October 18, 2013 at 4:48 pm
Here is a more general R implementation for finding highest density intervals of density functions that are built into R:
HDIofICDF = function( ICDFname , credMass=0.95 , tol=1e-8 , … ) {
# Arguments:
# ICDFname is R’s name for the inverse cumulative density function
# of the distribution.
# credMass is the desired mass of the HDI region.
# tol is passed to R’s optimize function.
# Return value:
# Highest density iterval (HDI) limits in a vector.
# Example of use: For determining HDI of a beta(30,12) distribution, type
# HDIofICDF( qbeta , shape1 = 30 , shape2 = 12 )
# Notice that the parameters of the ICDFname must be explicitly named;
# e.g., HDIofICDF( qbeta , 30 , 12 ) does not work.
# Adapted and corrected from Greg Snow’s TeachingDemos package.
incredMass = 1.0 – credMass
intervalWidth = function( lowTailPr , ICDFname , credMass , … ) {
ICDFname( credMass + lowTailPr , … ) – ICDFname( lowTailPr , … )
}
optInfo = optimize( intervalWidth , c( 0 , incredMass ) , ICDFname=ICDFname ,
credMass=credMass , tol=tol , … )
HDIlowTailPr = optInfo$minimum
return( c( ICDFname( HDIlowTailPr , … ) ,
ICDFname( credMass + HDIlowTailPr , … ) ) )
}
Examples of use:
> HDIofICDF( qbeta , shape1 = 30 , shape2 = 12 )
[1] 0.5780654 0.8448405
> HDIofICDF( qbeta , shape1 = 117.5 , shape2 = 115.5 )
[1] 0.4402789 0.5682842
(The function is available in the programs with my book :-)
Best wishes, et je vous envoie mes bien amicales pensées.
–John
October 18, 2013 at 10:48 pm
Merci, John! Pretty neat generic HPD function, I hope readers will take a look at it and jump on the opportunity!
X
October 21, 2013 at 7:04 am
You are assuming the density function is unimodal, which is not necessarily always the case. I just did a quick google search, and found Rob Hyndman’s paper, which looks pretty interesting: https://stat.ethz.ch/pipermail/r-help-es/attachments/20090731/d60e56a5/attachment-0001.pdf
October 21, 2013 at 2:49 pm
a fairly relevant remark, Yihui, thanks! (I hope all is well in Ames!)
October 21, 2013 at 3:47 pm
Right, the algorithm assumes a unimodal distribution (which is often, but not always, the case in real applications). The Hyndman article is cited in my book :-) Thanks for making the assumption explicit here!