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:
from scipy.stats import beta
beta_dist = beta.rvs(size=100000, a = 0.5, b = 0.5)
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.