How to calculate confidence intervals for a generic (non-normal) pdf 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:

*

*Is this right?

*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) )

*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.
 A: 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.
