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.

5 Responses to “beta HPD”

  1. 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

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

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