# How to Summarize a 2D Posterior Using a Highest Density Ellipse

**Publishable Stuff**, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Making a slight digression from last month’s Probable Points and Credible Intervals here is how to summarize a 2D posterior density using a *highest density ellipse*. This is a straight forward extension of the highest density interval to the situation where you have a two-dimensional posterior (say, represented as a two column matrix of `samples`

) and you want to visualize what region, containing a given proportion of the probability, that has the most probable parameter combinations. So let’s first have a look at a fictional 2D posterior by using a simple scatter plot:

plot(samples)

Whoa… that’s some serious over-plotting and it’s hard to see what’s going on. Sure, the bulk of the posterior is somewhere in that black hole, but where exactly and how much of it?

A highest posterior density ellipse shows this by covering the area that contains the most probable parameter combinations while containing *p*% of the posterior probability. Like finding the highest density interval corresponds to finding the *shortest* interval containing *p*% of the probability, finding the highest density ellipse corresponds to finding the *smallest* ellipse containing *p*% of the probability a.k.a. the minimum volume ellipse. I have spent a lot of time trying to figure out how compute minimum volume ellipses. Wasted time, it turns out, as it can be easily computed using packages that come with R, you just have to know what you are looking for. If you just want the code skip over the next paragraph, if you want to know the tiny bit of detective work I had to do to figure this out, read on.

To find the points in `sample`

that are included in a minimum volume ellipse covering, say, 75% of the samples you can use `cov.mve(samples, quantile.used = nrow(samples) * 0.75)`

from the MASS package, here `quantile.used`

specifies the number of points in `samples`

that should be inside the ellipse. It uses an approximation algorithm described by Van Aelst, S. and Rousseeuw, P. (2009) that is not guaranteed to find the minimum volume ellipse but that will often be pretty close. A problem is that `cov.mve`

does not return the actual ellipse, it returns a robustly measured covariance matrix, but that’s not really what we are after. It does return an object that contains the indices of the points that are covered by the minimum volume ellipse, if `fit`

is the object returned by `cov.mve`

then these points can be extracted like this: `points_in_ellipse <- samples[fit$best, ]`

. To find the ellipse we are going to use `ellipsoidhull`

from the cluster package on the `points_in_ellipse`

. It returns an object which represents the minimum volume ellipse and by using its `predict`

function we get a two column matrix with points that lie on the hull of the ellipse and that we can finally plot.

That wasn’t too easy to figure out, but it’s pretty easy to do. The code below plots a 75% minimum volume / highest density ellipse:

library(MASS) library(cluster) # Finding the 75% highest density / minimum volume ellipse fit <- cov.mve(samples, quantile.used = nrow(samples) * 0.75) points_in_ellipse <- samples[fit$best, ] ellipse_boundary <- predict(ellipsoidhull(points_in_ellipse)) # Plotting it plot(samples, col = rgb(0, 0, 0, alpha = 0.2)) lines(ellipse_boundary, col="lightgreen", lwd=3) legend("topleft", "50%", col = "lightgreen", lty = 1, lwd = 3)