Joint credible regions from MCMC draws

Question

Lets say I have $n$ posterior samples of $\theta_1$ and $\theta_2$. I suppose that any region $R$ which contains exactly $(1-\alpha)n$ of the points will be an approximate $(1-\alpha)\times100$ credible region for $(\theta_1, \theta_2)$. Can anybody make any suggestions, or point me to some literature for methods of constructing a good credible region $R$?

If it matters, this is a simulation study and I would like to evaluate the frequentist coverage of the region $R$ (I know there is no guarantee that the coverage will be nominal).

Off the top of my head, I can imagine an algorithm which checks all possible rectangular regions and gives some criteria for choosing one (smallest area?). It seems to me that we should be able to do better without enforcing a rectangular region. Any other method I can think of would use KDE, which intuitively bothers me since it could be sensitive to the choice of bandwith.

Answer 1

Indeed, the two approaches you mention are the most straightforward; I believe for example that the function HPDregionplot in the R package emdbook uses a kernel density estimator.

Another option (and this is just a suggestion) would be to find the mode or centre of your distribution, compute the distance of each point in your sample to that mode, and choose the $(1-\alpha)n$ points with smallest such distance. The Euclidian distance would give you a circular credible region; you can use another distance (maybe based on the empirical covariance matrix?) to get an elliptical region.

Answer 2

If you don't mind having a rectangular confidence region, and you know things are unimodal and symmetric, you could use the approach taken by credible.region() in the bayesSurv R package:
https://rdrr.io/cran/bayesSurv/src/R/credible.region.R

I believe this is the idea behind how it works:
We have $n$ MCMC samples for each of $p$ parameters. Imagine them plotted in $p$-dimensional space.
Start with the smallest hypercube around all the samples.
Shrink it down so it just barely excludes any point which is a marginal min or max along any axis.
(E.g. with $d=2$, we would remove up to 4 points: the top-most, left-most, right-most, and bottom-most points. Or maybe as few as just 2 points, if e.g. the top-most and right-most values are at the same point and likewise for left/bottom.)

Repeat the whole process on the remaining points.
Continue until only 95% (or whatever desired fraction) of samples remain.

(But the way it's implemented in practice, you only need to sort each column once---no need for a loop the way I described it here.)

This rectangular region will be larger than an ellipsoid, but not as computationally-intensive as relying on KDE or Euclidean distances.

Answer 3

I played around with a few different options, but I figured I'd share the one that I found worked best. Note: In my application, the posterior is well approximated as a multivariate normal. In other applications, the HPDregionplot approach may be more flexible.

There are essentially 3 steps...

1. Rotate the data along it's principal components. Scale each direction by the reciprocal square root of the eigenvalue corresponding to that direction.
2. In this new space, find the euclidean ball with the smallest radius $r$, that contains $(1-\alpha)\times n$ points.
3. Convert this "spherical" region back to an "elliptical" region in the original space.

A quick illustration of this process is given below.

---

Edit: This approach seems (roughly, if not exactly) equivalent to using Mahalanobis distance using the sample precision matrix $S^{-1}$. In other words, let $d_\star$ be the smallest possible $d$ such that $$\sum_{i=1}^nI\{({\bf x}_i - {\bf \hat\mu})^T S^{-1}({\bf x}_i - {\bf \hat\mu}) \leq d\} \geq n(1-\alpha),$$ where $I(b)$ is the indicator function (equal to $1$ when $b$ is true). Then $$R = \{{\bf x} \ | \ ({\bf x} - {\bf \hat\mu})^T S^{-1}({\bf x} - {\bf \hat\mu}) \leq d_\star\}$$ is an approximate $(1-\alpha)\times 100\%$ joint confidence region for ${\bf x}$.