2
$\begingroup$

In particular, this is a question about distributions with some skew and/or where the mean, median and mode might all be different.

For any pdf, it's possible to find multiple different domains within the support which contain exactly 95% (or another other %ge) of the integral. For distributions like the normal distribution, it's very intuitive that you want to find $b$ such that $\int _{\mu -b}^{\mu +b}f(x)dx = 0.95$ in order to find the 95% confidence interval.

For other distributions, it's less clear that you want to do this, because the mean $\mu$ of the distribution holds a less special status (in the normal case, it is the mean, median and mode).

It would seem to me, that really, one is trying to find (a,b) in a way that minimises (b-a) subject to the constraint that $\int_{a}^{b}f(x)dx=0.95$

In the case of the normal distribution, this will boil down to exactly what I wrote above.

My questions are:

  1. Is this right?
  2. Presumably this has to be done numerically. Is there a known algorithm for doing this? (naively I'd just vary a and find the corresponding b for every a (for many a there will exist no b) and then find the (a,b) pair that minimises (b-a) )
  3. Is there an open source implementation of this in python? Ideally I would pass two vectors representing a pdf, one which tells it at which points on the x-axis I've sampled the pdf and one which gives the density at each of these points.

I'm interested, for example, of finding the 95% confidence interval of the beta distribution for arbtirary $(\alpha, \beta)$. Of course it would be ideal if one could do this analytically but I strongly suspect this won't be possible.

I appreciate that things get even more complicated for multi-modal distributions, I'm happy with answers that only work for unimodal distributions for the moment.

Improve this question
$\endgroup$
5
  • 1
    $\begingroup$ I think this is what you are looking for: stats.stackexchange.com/questions/148439/…. On your points (1) It is almost correct because minimizing $(b-a)$ works only for unimodal distribution; (2) See this: stats.stackexchange.com/a/383626/239481; (3) Just search for HDI in python - it must be there for some bayesian statistics related package as it is used extensively there (called HPDI). $\endgroup$
    – Dayne
    Commented Dec 20, 2020 at 11:56
  • $\begingroup$ A little more explanation on (1): minimizing $(b-a)$ is forcing the interval to be continuous. For a multimodal distribution it maybe that a discontiguous set has smaller range and yet area under $f(x)$ for this set is 0.95. To get such sets, we find the highest horizontal bar such that area under $f(x)$ between pair(s) of intersection points is 0.95 $\endgroup$
    – Dayne
    Commented Dec 20, 2020 at 12:01
  • $\begingroup$ Thanks, useful pointers. Indeed the problem is substantially more complex for multimodal distributions. Surprisingly, googling "Highest Density Interval Python" has not been hugely fruitful (I found a few implementations that use raw data rather than the pdf itself)... but that's OK, at least this gives me confidence I'm not wasting my time implementing this (for unimodal cases) myself. $\endgroup$
    – gazza89
    Commented Dec 20, 2020 at 13:28
  • $\begingroup$ For python implementation you might find this helpful: stackoverflow.com/questions/22284502/… $\endgroup$
    – Dayne
    Commented Dec 20, 2020 at 13:31
  • $\begingroup$ kinda. The pymc3 example uses raw data and the scipy.stats.ppf example finds the central confidence interval rather than the HDI, so neither are quite what I'm after. I guess a somewhat inelegant solution would be to sample a tonne of data from from pdf and then pass to pymc3.stats.hpd $\endgroup$
    – gazza89
    Commented Dec 20, 2020 at 13:37

1 Answer 1

Reset to default
0
$\begingroup$

You can find an analysis of computation of HDRs of continuous univariate distributions in O'Neill (2022). This paper frames the problem as an optimisation problem that can be solved using standard calculus methods, and this usually requires numercal methods to find the optimising point. As you point out, this problem is relatively simple in the case of a unimodal distribution, but it is complicated for multimodal distributions. You are also correct that the interval may be non-symmetric in cases where you have a non-symmetric distribution. The linked paper has a detailed analysis of the optimisation that arises in the unimodal case and it shows how to construct an algorithm to compute the HDR.

As to computational implementation, this method is presently available in R in the HDR.unimodal function in the stat.extend package. That HDR function requires you to have access to the quantile function for the distribution (which is usually available for all classes of distributions of interest). I am not familiar enough with Python to know if there is any corresponding function in that language, but it could be programmed from scratch using the method in the linked paper.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.