## 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)
 0.4504879
> f(ff(.95))
 0.5580267
```

which was enough for my simple book illustration… Since (.450,558) is then the HPD region at credible level 0.95.

### 5 Responses to “beta HPD”

1. John K. Kruschke Says:

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 )
 0.5780654 0.8448405
> HDIofICDF( qbeta , shape1 = 117.5 , shape2 = 115.5 )
 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

• xi'an Says:

Merci, John! Pretty neat generic HPD function, I hope readers will take a look at it and jump on the opportunity!
X

• Yihui Says:

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

• xi'an Says:

a fairly relevant remark, Yihui, thanks! (I hope all is well in Ames!)

• John K. Kruschke Says:

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!

This site uses Akismet to reduce spam. Learn how your comment data is processed.