Highest Density Regions (HDR) misunderstanding

Questions
1

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!

2

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.

3

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.

4

This is the one:

It’s from @aloctavodia’s book.

2 Likes
5

Thanks a lot @junpenglao for finding this one!

6

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

7

The code above is truncated. The original is gone.

8

@aloctavodia might have an improved version now, pinging him.

9

I will add it to ArviZ ASAP. this the new link https://github.com/aloctavodia/BAP/blob/02d559a664ee7a55ca5d3e8f5af58e49c40dee75/first_edition/code/Chp1/plot_post.py

1 Like
10

BTW, I think I have a more robust method, but I need to test it to see if actually works as expected.

1 Like