Calculating the minimum value of a higest density interval from a posterior paramter estimate

Question

I have been working through Doing Bayesian Data Analysis by John K. Kruschke. In this book he has created a function for calculating the minimum and maximum values of the Highest Density Interval (HDI) of a mcmc chain, for use in making decisions about whether a the true parameter value of a given effect is non-zero. I deconstructed this function and am confused about the logic of a few steps.

Instead of running an entire JAGS model I have created a vector that will serve in the place of a true mcmc chain

v <- rnorm(50000, 0.6, .1) 

The next steps I am clear about

vSort <- sort(v) # sort from lowest estimate to highest
ciInd <- .95*length(vSort) # calculate 95% of the total number of estimates
int5 <- length(vSort) - ciInd # 5% of the total number of estimates

Now these are the steps I do not understand. We generate a vector of values as long as each side of the HDI. The 1st value generated by the loop is the value at the 47501st position in the vector of sorted estimates minus the value at the 1st position. The 2nd loop generates a value that is the value at 47502nd position minus the value at the 2nd position etc.

vecHDI <- sapply(1:int5, function (i) vSort[i + ciInd] - vSort[i])

To derive the minimum value of the HDI we find the position of the lowest value in the vector we just created (which represents the lowest difference between estimates 47500 positions apart) and substitute that into the original sorted vector

(HDImin <- vSort[which.min(vecHDI)]) 

To find the maximum HDI value derived we simply add 47500 to the position of the minimum value.

(HDImax <- vSort[which.min(vecHDI) + ciInd]) 

What I am unclear about is why we go through the process of generating the vector of differences to derive the reference for the position of the miminimum HDI? Why do we not just use the 2501st position in the sorted vector of parameter estimates (i.e. vSort[int5 + 1])?

Answer 1

The short answer is that these would (usually) be calculating different things - see the following graph which shows highest posterior density interval (using coda for convenience) in blue and a confidence interval based on a simple quantile as you propose in red:

library('coda')
chain <- mcmc(rgamma(10000, 2, 1))

hpd <- HPDinterval(chain)
ci <- quantile(chain, prob=c(0.025,0.975))

curve(dgamma(x, 2, 1), from=0, to=6)
abline(v=ci, col='red')
abline(a=0.2, b=-0.032, col='red')
abline(v=hpd, col='blue')
abline(h=0.04, col='blue')

The blue lines will always be parallel to the axes whereas the red line joining the lower and upper CI points will not be parallel to the x-axis when the distribution is e.g. skewed (although the two types of confidence intervals will be the same in the case of the normal distribution which may be contributing to your confusion).

It is worth reading up a bit more on the motivation behind highest posterior density intervals - see for example the accepted answer to this question: What is a Highest Density Region (HDR)?