How can I estimate the highest posterior density interval from a set of x,y values describing the PDF?

Question

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)

library(ggplot2)

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

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

library(dplyr)

sum <- df %>%
  group_by(x) %>%
  summarise("y"=sum(y*weight))

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

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)

> 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:

Thanks for any help you can provide!

Answer 1

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:

library(ggplot2)

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

Answer 2

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.

Analysis

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.

Implementation

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.

Examples

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:

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.

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:

library(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)) +
  theme(panel.grid=element_blank())

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.