Introduction

In Bayesian statistics, a credible interval is an interval within which an unobserved parameter falls with a certain probability. While there are many ways to construct credible intervals, the Highest Posterior Density (HPD) interval is particularly important because it represents the shortest possible interval for a given probability level.

Definition of HPD

The Highest Posterior Density (HPD) region of credibility \(100(1-\alpha)\%\) is the subset \(C\) of the parameter space \(\Theta\) such that:

\[P(\theta \in C) = 1-\alpha\]

For any \(\theta_1 \in C\) and \(\theta_2 \notin C\), \(p(\theta_1) \geq p(\theta_2)\)

In other words, the HPD region contains all points with the highest posterior density that collectively have probability \(1-\alpha\).

For a unimodal posterior distribution, the HPD interval satisfies two conditions:

Probability condition: The interval contains exactly \(100(1-\alpha)\%\) of the posterior probability
Density condition: The density at the lower bound equals the density at the upper bound (\(p(L) = p(U)\))

Why HPD Matters

The HPD interval is optimal in the sense that:

It is the shortest credible interval for a given probability level
It provides the most precise inference for the parameter
It focuses on the most probable regions of the parameter space

Implementation

Let’s demonstrate the concept using a Gamma posterior distribution with shape parameter \(\alpha=3\) and rate parameter \(\beta=0.7\), which might arise from a conjugate prior for a Poisson likelihood.

1. Posterior Distribution

First, let’s define our posterior distribution:

# Define the Gamma posterior distribution
shape <- 3
rate  <- 0.7
prob  <- 0.95

# Create a grid of values
x <- seq(0, qgamma(0.9999, shape, rate), length.out = 5000)

# Calculate density
density <- dgamma(x, shape, rate)

posterior <- data.frame(
  theta = x,
  density = density
)

# Plot the posterior
ggplot(posterior, aes(theta, density)) +
  geom_line(linewidth = 1.2) +
  labs(
    title = "Gamma Posterior Distribution",
    x = expression(theta),
    y = "Posterior Density"
  ) +
  theme_bw(base_size = 14)

2. Function to Calculate Interval Probability

# Function to calculate probability of an interval
interval_prob <- function(lower, upper) {
  pgamma(upper, shape, rate) -
    pgamma(lower, shape, rate)
}

3. Creating Arbitrary 95% Intervals

We can create different 95% credible intervals by choosing different lower tail probabilities:

# Function to create 95% intervals
make_95_interval <- function(alpha) {
  lower <- qgamma(alpha, shape, rate)
  upper <- qgamma(alpha + prob, shape, rate)
  
  data.frame(
    lower = lower,
    upper = upper,
    length = upper - lower,
    probability = interval_prob(lower, upper)
  )
}

4. Computing the HPD Interval

For a unimodal posterior, the HPD interval can be found by solving:

\[P(L<\theta<U)=0.95\] \[p(L)=p(U)\]

# Function to find HPD interval
hpd_function <- function(lower) {
  upper <- qgamma(
    pgamma(lower, shape, rate) + prob,
    shape,
    rate
  )
  
  dgamma(lower, shape, rate) -
    dgamma(upper, shape, rate)
}

# Find the HPD lower bound
HPD.lower <- uniroot(
  hpd_function,
  lower = 0,
  upper = qgamma(0.05, shape, rate)
)$root

# Calculate the HPD upper bound
HPD.upper <- qgamma(
  pgamma(HPD.lower, shape, rate) + prob,
  shape,
  rate
)

# Create HPD interval data frame
HPD <- data.frame(
  Type = "HPD (shortest)",
  lower = HPD.lower,
  upper = HPD.upper
)

HPD$length <- HPD$upper - HPD$lower
HPD$probability <- interval_prob(HPD$lower, HPD$upper)

Comparing Different 95% Credible Intervals

Let’s create and compare several 95% credible intervals:

# Create different intervals
Equal <- make_95_interval(0.025)
Left_shift <- make_95_interval(0.005)
Right_shift <- make_95_interval(0.045)

# Add labels
Equal$Type <- "Equal tail"
Left_shift$Type <- "Shifted left"
Right_shift$Type <- "Shifted right"

# Combine all results
results <- bind_rows(
  HPD,
  Equal,
  Left_shift,
  Right_shift
)

# Format results
results <- results %>%
  select(Type, lower, upper, length, probability) %>%
  mutate(
    length = round(length, 3),
    probability = round(probability * 100, 2),
    lower = round(lower, 3),
    upper = round(upper, 3)
  )

# Display results
knitr::kable(results,
  caption = "Comparison of Different 95% Credible Intervals",
  digits = 3,
  align = "c"
)

Comparison of Different 95% Credible Intervals
Type	lower	upper	length	probability
HPD (shortest)	0.434	9.145	8.711	95
Equal tail	0.884	10.321	9.437	95
Shifted left	0.483	9.199	8.717	95
Shifted right	1.118	13.248	12.130	95

6. Visualizing the Intervals

Let’s visualize how these different intervals cover the posterior distribution:

# Prepare data for plotting intervals
interval_data <- list(
  HPD = HPD,
  Equal = Equal,
  Left = Left_shift,
  Right = Right_shift
)

plot_data <- data.frame()

for(i in names(interval_data)) {
  tmp <- interval_data[[i]]
  
  tmp_density <- posterior %>%
    filter(theta >= tmp$lower &
             theta <= tmp$upper)
  
  tmp_density$Type <- tmp$Type
  
  plot_data <- bind_rows(
    plot_data,
    tmp_density
  )
}

# Create the plot
ggplot(posterior, aes(theta, density)) +
  geom_line(linewidth = 1.2) +
  geom_area(
    data = plot_data,
    aes(theta, density),
    fill = "steelblue",
    alpha = 0.35
  ) +
  geom_vline(
    data = results,
    aes(xintercept = lower, colour = Type),
    linewidth = 1
  ) +
  geom_vline(
    data = results,
    aes(xintercept = upper, colour = Type),
    linewidth = 1
  ) +
  facet_wrap(~Type, ncol = 2) +
  labs(
    title = "Same 95% Posterior Probability, Different Interval Lengths",
    subtitle = "HPD interval is the shortest credible interval",
    x = expression(theta),
    y = "Posterior density"
  ) +
  theme_bw(base_size = 12) +
  theme(
    legend.position = "none",
    strip.background = element_rect(fill = "lightgray")
  )

7. Comparing Interval Lengths

The HPD interval should be the shortest. Let’s verify this:

ggplot(results, aes(x = reorder(Type, length), y = length)) +
  geom_col(fill = c("steelblue", "gray70", "gray70", "gray70")) +
  geom_text(
    aes(label = length),
    vjust = -0.5,
    size = 4
  ) +
  labs(
    title = "Length of 95% Credible Intervals",
    subtitle = "HPD provides the shortest credible interval",
    x = "Interval Type",
    y = "Interval Length"
  ) +
  theme_bw(base_size = 14)

Key Takeaways

HPD intervals are the shortest credible intervals for a given probability level
For symmetric unimodal posteriors (like normal), HPD equals the equal-tail interval
For skewed posteriors (like gamma), the equal-tail interval is not HPD
The HPD interval focuses on the most likely regions of the parameter space
HPD intervals are particularly useful when the posterior is skewed or multimodal

Practical Considerations

For unimodal posteriors, the HPD interval can be found by solving the density equality condition
For multimodal posteriors, the HPD region might consist of multiple disjoint intervals
In practice, HPD intervals are often computed using MCMC samples and density estimation
The HPD interval is invariant under monotone transformations of the parameter

Highest Posterior Density (HPD) Intervals

2026-07-08