I have a dataset where values are distributed in some experimental way, not following any theoretical distribution. The dataset could be unimodal with a nicely defined single peak or multimodal with dips or various depths. What would be a fast algorithm to approximate the 80-90-95-99% (pick one) credible region? Computational performance is valued SIGNIFICANTLY over accuracy. Accuracy to +/-1% is reasonable, as the data already has significant bias/variance in it. 

If you're not familiar with credible regions, it is the interval (or multiple intervals!!!) where the highest probability mass occurs. The slow way to compute this is `NIntegrate[]` as the sum of the areas under the tallest parts of PDF should add up to the level you are trying to achieve (e.g. 95%).

Visually speaking, the question is "How high/low should we set the bar, so that the area under intersecting intervals is ~0.95(+/-1)%?"

Currently, my top idea is to divide data into a large number of bins (say N/100), then sort bins and keep adding largest bins to a list until the total of selectel bin counts is greater than the desired credibility. Then sort bins by their dependent variable value and group any continuous values into a single interval, thus creating the list of intervals which span the highest density region. Am I overcomplicating things? And is there perhaps a built in function, of which I'm unaware, which would solve this task out-of-the-box?

**Unimodal CI vs HDR:**
[![Unimodal confidence interval vs highest density region](https://i.stack.imgur.com/okCu3.png)](https://i.stack.imgur.com/okCu3.png)

**Multimodal CI vs HDR:**
[![Multimodal confidence interval vs highest density region](https://i.stack.imgur.com/3Taa9.png)](https://i.stack.imgur.com/3Taa9.png)