I'd like to estimate the Highest Posterior Density Interval (HPDI) of a calculated density function, rather than from empirical samples as is normally done (e.g., from an mcmc object). The situation is that I am aggregating two probability functions and I end up with x,y values describing a non-parametric probability distribution. I'd like to then know the 80% HPDI of this aggregate function, rather than a symmetric 80% density interval, which I could calculate from the CDF.

Here is a setup for an example:

I have two theoretical probability distributions:

d1 <- data.frame(x=seq(-10,10), y=dnorm(x = seq(-10,10), mean=-4, sd=2), group="A")
d2 <- data.frame(x=seq(-10,10), y=dnorm(x = seq(-10,10), mean=-1, sd=3), group="B")

df <- rbind(d1,d2)


p <- ggplot(df, aes(x,y, group=group, colour=group)) +
  geom_line() +

To aggregate them, I take the arithmetic mean of the two distributions at each value of x like so, and plot the aggregate function:

df$weight <- 1/length(unique(df$group))


sum <- df %>%
  group_by(x) %>%

p.a <- p + geom_line(data=sum, aes(x,y), inherit.aes=F)

I can access the x,y values of the aggregate function using ggplot_build():

agg <- ggplot_build(p.a)$data[[2]][,1:2]

> head(agg,3)
    x           y
1 -10 0.001846604
2  -9 0.006281406
3  -8 0.017868056

I'm focused heavily on implementation rather than theory, so I'd prefer responses focused on functions to help do this.

A figure of the two distributions and the aggregate (in black) in case it's helpful as a visual:

enter image description here

Thanks for any help you can provide!

  • $\begingroup$ The scope of your question isn't apparent. You seem to describe an equally weighted mixture of two Normal distributions. The sense of "nonparametric" is unclear, because that would mean this mixture comes from a rather large family of distributions, but no family is in evidence. What, then, do you need to compute? HPD intervals for mixtures of Normals? HPD intervals for mixtures generally? HPD intervals for distributions whose densities are given as a sequence of vertices on their graphs? $\endgroup$
    – whuber
    Dec 11 '18 at 21:16
  • $\begingroup$ @whuber, Apologies for my loose use of language here. You're right my example is asking about mixtures of normals, but in my actual application, I can have mixtures of different distribution types. I meant "non-parametric" in the sense that the resulting mixtures can't be described as belonging to any particular family of distributions. In effect, what I'm hoping for is the last of what you describe - HDIs for distributions described by a sequence of vertices. Hopefully that clarifies? $\endgroup$
    – phalteman
    Dec 11 '18 at 21:32
  • $\begingroup$ Do you need the single 80% interval of highest density, or the highest-density set with 80% coverage, which need not be a single interval? $\endgroup$ Dec 17 '18 at 15:57
  • $\begingroup$ @Kodiologist, The single X% interval of highest density is what I'm looking for. $\endgroup$
    – phalteman
    Dec 17 '18 at 19:14
  • $\begingroup$ Must the solution truly be an interval or should it be a region in general? If it needn't be an interval, then there are some extremely efficient solutions available in many cases. $\endgroup$
    – whuber
    Dec 17 '18 at 22:20

Since you have some example points from the PDF rather than a closed-form representation, and you're looking for a highest-density interval rather than a highest-density set, it's easiest to do this by brute force. Consider all $a$ and $b$ among the $x$-coordinates you have, check that $[a, b]$ has at least the desired coverage, and return the interval of greatest density. Here's a simple-minded implementation that (a) is slow for large numbers of points and (b) treats the density of each $x_n$ as if it was the density between $x_{n-1}$ and $x$.

hdi = function(x, x.density, coverage)
   {best = 0
    for (ai in 1 : (length(x) - 1))
        {for (bi in (ai + 1) : length(x))
            {mass = sum(diff(x[ai : bi]) * x.density[(ai + 1) : bi])
             if (mass >= coverage && mass / (x[bi] - x[ai]) > best)
                {best = mass / (x[bi] - x[ai])
                 ai.best = ai
                 bi.best = bi}}}
    c(x[ai.best], x[bi.best])}

An example:


x = seq(0, 1, len = 1000)
x.density = dbeta(x, shape1 = 10, shape2 = 2)

interval = hdi(x, x.density, .8)

qplot(x, x.density) + geom_vline(aes(xintercept = interval))

A plot of the example

  • $\begingroup$ I think that does it! Assumption (b) is a safe approximation for my purposes here, and I'm not particularly worried about (a). Nice, simple, solution! Thanks! $\endgroup$
    – phalteman
    Dec 17 '18 at 21:28
  • $\begingroup$ Assuming the solution is an interval, I believe your $O(n^2)$ algorithm could be implemented in $O(n\log(n))$ time by sorting the densities. $\endgroup$
    – whuber
    Dec 17 '18 at 22:23

Depending on your needs, there are efficient solutions available. The purpose of this answer is to show how you can formulate the problem, encapsulate the key ideas, and implement them using capabilities offered by your programming environment.


Given any density function $f$, all highest posterior density sets are of the form $\{x\mid f(x) \ge h\}$ or $\{x\mid f(x) \gt h\}.$ For floating-point computing purposes we needn't distinguish between them, so let's just call either one $f^{[h]}.$ Writing the indicator function as $\mathcal I,$ the total probability of any such set is

$$p_f(h) = \Pr\left(f^{[h]}\right) = \int_{\mathbb{R}} \mathcal{I}(f(x) \gt h) f(x)\mathrm{d}x.$$

Obtaining the $100\alpha\%$ highest probability density set therefore is a matter of solving

$$0 = p_f(h) - \alpha.$$

This reduces the problem to one of root finding, which has well-established solutions and is widely implemented using accurate, efficient methods.


Here is a naive R solution. ("Naive" means it contains no defenses against difficult-to-evaluate pdfs or bad inputs, such as when x.min and x.max do not include the entire support of the distribution.) It employs the function ifelse to implement $\mathcal I,$ integrate to integrate the pdf, and uniroot to find a root. Its arguments are the amount of probability alpha, a density function df, and the support of df (smallest and largest possible values).

highest_alpha <- function(alpha, df, x.min, x.max, ...) {
  p <- function(h) {
    g <- function(x) {y <- df(x); ifelse(y > h, y, 0)}
    integrate(g, x.min, x.max, ...)$value - alpha
  uniroot(p, c(x.min, x.max), tol=1e-12)$root

The question supposes the pdf is given in a somewhat awkward form as a sequence of vertices along its graph: its "spaghetti representation." No problem: just interpolate. Although in many cases an interpolation of the log pdf would be best, some care is needed to handle areas where the pdf might equal zero, so in this example I just use linear interpolation as implemented by approxfun. I do take measures to ensure the pdf is normalized, though, by first integrating it to find a normalizing constant.

as.pdf <- function(x, y, ...) {
  f <- approxfun(x, y, method="linear", yleft=0, yright=0, rule=2)
  const <- integrate(f, min(x), max(x), ...)$value
  approxfun(x, y/const, method="linear", yleft=0, yright=0, rule=2)

Its arguments are the x and y coordinates, sorted by x, of the spaghetti representation of the pdf.

These six lines of code will perform asymptotically well, typically requiring $O(n\log(n))$ calculations when the pdf is given by $n$ vertices of its graph.


As a test, highest_alpha was applied to a standard Normal pdf as shown in the first figure. The total probability was divided into the highest 1/6, highest 2/6, ... highest 6/6 and the corresponding areas have been colored accordingly. By symmetry, the boundaries must lie on the set of points $\Phi^{-1}(1/12), \Phi^{-1}(2/12), \ldots, \Phi^{-1}(11/12)$ where $\Phi$ is the standard normal CDF, so I have plotted those points as vertical black line segments:

Figure 1

The black dots, of course, show the spaghetti representation of the pdf that was used. The horizontal lines are the values of $h$ found by highest_alpha.

Let's see a solution for a more complex pdf. This one is a mixture of three Normal distributions.

Figure 2

A close look shows that indeed the interpolation is working: many of the region boundaries fall between the spaghetti points.

For completeness, here is the code used for the first example. The first four lines create the spaghetti representation while the last two do the calculations.

i <- c(exp(-(1:10)), 1 - exp(-(1:10)), seq(0, 1, length.out=101))
x <- sort(qnorm(i))
x <- x[!is.infinite(x)]
y <- dnorm(x)
k <- 6
h <- sapply(1:(k-1) / k, function(h) highest_alpha(h, as.pdf(x, y), min(x), max(x)))

The figure is made with ggplot2:

n <- 501
X <- data.frame(x=seq(min(x), max(x), length.out=n))
X$y <- as.pdf(x, y)(X$x)
X$Interval <- factor(rowSums(outer(X$y, h, ">")))
dx <- diff(range(x)) / n
ggplot(X, aes(x, y)) + 
  geom_hline(aes(yintercept=h, color=Density), size=1, show.legend=FALSE,
             data=data.frame(h=h, Density=factor(h, labels=signif(h, 2)))) + 
  geom_vline(xintercept=qnorm(1:(2*k-1)/(2*k))) +
  geom_path(color="gray") + 
  geom_col(aes(fill=Interval, color=NULL), alpha=0.5, width=dx) + 
  geom_point(data=data.frame(x=x, y=y)) + 
  scale_fill_manual(values=terrain.colors(k)) + 
  scale_color_manual(values=terrain.colors(k)) +

I dodged one issue by computing explicitly only the heights corresponding to the highest-density regions. To find those regions, one has to determine where the pdf crosses those heights. Because that's just another root-finding exercise, I won't go into the (redundant) details.

If you really must find an interval (a connected set), then you can formulate the problem as one of minimizing the interval length $\delta$ subject to the highest-probability constraint. That is, given $0 \lt \alpha \lt 1,$ the problem is to find an ordered pair $(x, \delta)$ where $\delta$ is a small as possible subject to the constraints

  1. $\delta \gt 0.$

  2. $x$ and $x+\delta$ are in the domain of the pdf $f.$

  3. $F_\alpha(x,\delta) \ge 0$ where $$F_\alpha(x,\delta) = \int_x^{x+\delta} f(x)\mathrm{d}x - \alpha.$$

This problem is amenable to the same approach: (1) use the capabilities of your programming environment to compute the integral accurately and efficiently and (2) employ an efficient two-dimensional constrained optimization routine. Most of the constraints are linear and the single non-linear constraint is differentiable, implying this approach has a good chance of succeeding and yielding accurate solutions.

  • 1
    $\begingroup$ Thanks for the very detailed response. For my purposes, I do need a true interval, as noted in the comments, but this is an excellent resource for anyone (including myself in the future) who is looking for a more universally relevant solution. $\endgroup$
    – phalteman
    Dec 19 '18 at 19:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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