## beta HPD

**W**hile 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.

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!