# How to get multivariate credible interval estimate(s) / highest density regions (HDR) after MCMC

I'm estimating 15 parameters of my model using a Bayesian approach and a Markov Chain Monte Carlo (MCMC) method. My data after running a MCMC chain of 100000 samples is therefore a 100000×15 table of parameter values.

I want to find 15-dimensional Highest Density Regions of my posterior distribution.

My problem: Clustering the samples to assign them to HDRs (example using density-based clustering below) needs a distance matrix of all the samples. For 100000 samples, this matrix would take 37 GiB of RAM, which I do not have, not speaking of the computing time. How can I find my HDRs using reasonable amounts of computing resources? Someone must have had this problem before?!

Edited to add: According to this SO question and DBSCAN Wikipedia page, DBSCAN can be brought down to $\mathcal O(n\log n)$ time complexity and $\mathcal O(n)$ space complexity using a spatial index and avoiding a distance matrix. Still looking for an implementation or a description thereof...

## Multivariate highest density regions using density-based clustering (DBSCAN)

A X% highest density region is the region of a distribution which encompasses X% of the probability mass. As samples drawn by a MCMC method apprear with a frequency (asymptotically) proportional to the searched-for posterior distribution, my X% HDR also encompasses X% of my samples.

I planned to use the density-based clustering algorithm DBSCAN to cluster my samples because density of samples is directly related to peak height of my posterior.

In analogy by the method of Hyndman (1996) (paper, SO question), I planned to increase the maximum distance a single sample may have from a cluster to be regarded a part of it iteratively until X% of my samples are part of some cluster: After that step, I would calculate the range of each cluster in each dimension as a way to present my Highest Density regions.

In this example, you would be able to see that the 80% HDR encloses two distinct regions while the 50% HDR only contains one cluster. I would be able to visualize this like shown below, because the plot above is not applicable to more than 2 dimensions: • When you run an MCMC you observe the value of the target at each iteration (up to the same normalising constant), hence you can identify high density simulations in a much more precise way than by frequencies. You can also identify HPD regions by running a standard clustering algorithm on simulations with the highest posterior density values. – Xi'an Oct 7 '16 at 16:34
• Thank you for your suggestion, @Xi'an. But can you specify what you think is a "standard clustering algorithm"? What algo could you suggest other than DBSCAN (or {single|complete}-linkage clustering, which is the same idea) if my notions of clusters are "points of high density"? Moreover, the most samples are in the regions of highest density, therefore If I want to calculate, say, the 90% HDI, I can only expect to exclude about 10% of my samples. – akraf Oct 11 '16 at 9:45
• Samworth & Wand (2010, Annals of Statistics) may offer an alternative approach. – S. Kolassa - Reinstate Monica Jul 31 '18 at 10:00

I have found a Matlab Wrapper for ANN. ANN is a library for approximate nearest neighbor searches (homepage). Besides the usual parameters of a spatial index region query, it uses an additional error parameter eps which gives the "approximateness" of the search: A returned nearest neighbor will be at most 1+eps times farther from the query point than the true (non-approximate) nearest neighbor. Search for the term "error bound" in the Programmers Manual to find information about eps.