# Given a 10D MCMC chain, how can I determine its posterior mode(s) in R?

Question: With a 10 dimensional MCMC chain, let's say I'm prepared to hand you a matrix of the draws: 100,000 iterations (rows) by 10 parameters (columns), how best can I identify the posterior modes? I'm especially concerned with multiple modes.

Background: I consider myself a computationally savvy statistician, but when a colleague asked me this question, I was ashamed that I couldn't come up with a reasonable answer. The primary concern is that multiple modes may appear, but only if at least eight or so of the ten dimensions are considered. My first thought would be to use a kernel density estimate, but a search through R revealed nothing promising for problems of greater than three dimensions. The colleague has proposed an ad-hoc binning strategy in ten dimensions and searching for a maximum, but my concern is that bandwidth might either lead to significant sparsity problems or to a lack of resolution to discern multiple modes. That said, I'd happily accept suggestions for automated bandwidth suggestions, links to a 10 kernel density estimator, or anything else which you know about.

Concerns:

1. We believe that the distribution may be quite skewed; hence, we wish to identify the posterior mode(s) and not the posterior means.

2. We are concerned that there may be several posterior modes.

3. If possible, we'd prefer an R based suggestion. But any algorithm will do as long as it isn't incredibly difficult to implement. I guess I'd prefer not to implement a N-d kernel density estimator with automated bandwidth selection from scratch.

Have you considered using a nearest neighbour approach ?

e.g. building a list of the k nearest neighbours for each of the 100'000 points and then consider the data point with the smallest distance of the kth neighbour a mode. In other words: find the point with the 'smallest bubble' containing k other points around this point.

I'm not sure how robust this is and the choice for k is obviously influencing the results.

• Sometimes I just want to thwack myself upside the head. Excellent suggestion. Oct 5 '10 at 12:28
• I also just thought of using the kmeans function in R. I really shouldn't ask questions between midnight and 4am. Oct 5 '10 at 13:44

This is only a partial answer.

I recently used figtree for multidimensional kernel density estimates. It's a C package and I got it to work fairly easily. However, I only used it to estimate the density at particular points, not calculate summary statistics.

If you keep the log likelihoods, you can just select the one with the highest value. Also, if your interest is primarily the mode, just doing an optimization to find the point with the highest log likelihood would suffice.

• This is the most relevant answer, at least the first part! In many MCMC simulations, the (log-)likelihoods are computed for all proposals and can thus be stored. Or the highest value so far and its argument can be stored. Provided the MCMC algorithm has converged over the number of simulations you ran, this is a valid approach. Jan 12 '12 at 20:22

Have you considered 'PRIM / bump hunting' ? (see e.g. Section 9.3. of 'The Elements of Statistical Learning' by Tibshirani et al. or ask your favourite search engine). Not sure whether that's implemented in R though.

[ As far as I understood are you trying to find the mode of the probability density from which your 100'000 rows are drawn. So your problem would be partially solved by finding an appropriate density estimation method ].

• Yes, there's a prim package, with an R vignette: Using prim for bump hunting. It's not obvious to me how it will work in this case, though.
– chl
Nov 4 '10 at 22:36