Highest Density Regions (HDR) misunderstanding



I came across pymc3.stats.hpd() as a good solution for obtaining HDR. I have one question:

  • How should I interpret the output of pymc3.stats.hpd()? For instance, in my use case, a get array([0.34652132, 1.])

The output for HDR mentioned in the question actually comes from finding HDR on beta distribution with a=0.5, b=0.5. Thus, if the output from the question is correct, than this output does not have any sense to me at the moment since the distribution in my case is bimodal. I should get two intervals, not only one (If the output in my question is considered to be starting and ending points of the HDR interval).

Here is my code if someone is interested:

import pymc3
from scipy.stats import beta

beta_dist = beta.rvs(size=100000, a = 0.5, b = 0.5)
print(pymc3.stats.hpd(beta_dist, alpha=0.4))

Thanks in advance!


I think the current implementation only works for single mode marginal, as you observed for bimodal where the modes are far away it will gives incorrect answer.
I remember seeing a version that works for multiple mode from @aloctavodia, maybe he can provide a bit more information.


Thanks @junpenglao for the quick response! :slight_smile:

I hope will get an answer from the @aloctavodia, perhaps he knows a bit more regarding the bimodal/multimodal case, as you said.


This is the one:

It’s from @aloctavodia’s book.


Thanks a lot @junpenglao for finding this one!


Is it on the roadmap? Or is this kind of thing going into Arviz now?