I promised my Bayesian students that I’d illustrate the hpd function in the TeachingDemos package.

Suppose I’m interested in estimating p, the proportion of undergraduate BGSU students whose home is 200 miles or further from Bowling Green. Initially, I place a uniform prior on p. I collect data from 16 students — only 1 lives 200+ miles from BG. My posterior for p is beta(2, 16).

I’ll assume that you have already installed the TeachingDemos package in R. Then I’ll simulate 10,000 draws from the beta(2, 16) distribution and use the function emp.hpd to find a 90% interval estimate (you use hpd if you actually have the quantile function of distribution available).

> library(TeachingDemos)

> p=rbeta(10000,2,16)

> emp.hpd(p,conf=0.90)

[1] 0.005806384 0.211324962

You could also find an equal-tails interval estimate:

> qbeta(c(.05,.95),2,16)

[1] 0.02131763 0.25012445

But the length of the equal-tails interval is 0.229 which is substantially larger than the length of the HPD interval (0.206). There will be a difference between the two intervals when the posterior density is asymmetric as in this situation.