Hi,

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!