SAMPLE_HPD_CI

This function approximates a highest posterior density (HPD) credible interval from posterior samples by finding the shortest interval that contains the requested posterior probability mass.

Given sorted samples x_{(1)} \leq \dots \leq x_{(n)} and confidence level c, let k = \lceil c n \rceil. The HPD approximation chooses index i minimizing interval width:

w_i = x_{(i+k-1)} - x_{(i)}

and returns the corresponding bounds (x_{(i)}, x_{(i+k-1)}). This is a practical sample-based approximation often used when analytic HPD limits are unavailable.

Excel Usage

=SAMPLE_HPD_CI(samples, confidence)
  • samples (list[list], required): 2D array of posterior sample draws (unitless).
  • confidence (float, optional, default: 0.95): Credible interval probability level (0 to 1).

Returns (list[list]): 2D table with sample median and approximate HPD lower/upper bounds.

Example 1: HPD approximation for symmetric posterior samples

Inputs:

samples confidence
-3 -2 -1 0 1 2 3 0.8

Excel formula:

=SAMPLE_HPD_CI({-3,-2,-1,0,1,2,3}, 0.8)

Expected output:

Median Lower Upper
0 -3 2
Example 2: HPD approximation for right-skewed posterior samples

Inputs:

samples confidence
0.1 0.2 0.3 0.5 1 2.5 4 0.9

Excel formula:

=SAMPLE_HPD_CI({0.1,0.2,0.3,0.5,1,2.5,4}, 0.9)

Expected output:

Median Lower Upper
0.5 0.1 4
Example 3: HPD approximation on clustered sample draws

Inputs:

samples confidence
1 1.1 1.2 3.8 3.9 4 4.1 0.75

Excel formula:

=SAMPLE_HPD_CI({1,1.1,1.2,3.8,3.9,4,4.1}, 0.75)

Expected output:

Median Lower Upper
3.8 1.1 4.1
Example 4: Multi-row posterior sample layout

Inputs:

samples confidence
2 2.2 0.85
2.5 2.7
3 3.1
3.4 3.6

Excel formula:

=SAMPLE_HPD_CI({2,2.2;2.5,2.7;3,3.1;3.4,3.6}, 0.85)

Expected output:

Median Lower Upper
2.85 2 3.4

Python Code

import numpy as np

def sample_hpd_ci(samples, confidence=0.95):
    """
    Approximate a highest posterior density interval from posterior samples using the narrowest empirical window.

    See: https://numpy.org/doc/stable/reference/generated/numpy.quantile.html

    This example function is provided as-is without any representation of accuracy.

    Args:
        samples (list[list]): 2D array of posterior sample draws (unitless).
        confidence (float, optional): Credible interval probability level (0 to 1). Default is 0.95.

    Returns:
        list[list]: 2D table with sample median and approximate HPD lower/upper bounds.
    """
    try:
        def to2d(x):
            return [[x]] if not isinstance(x, list) else x

        samples = to2d(samples)

        if not isinstance(samples, list) or not all(isinstance(row, list) for row in samples):
            return "Error: Invalid input - samples must be a 2D list"

        if not (0 < confidence < 1):
            return "Error: confidence must be between 0 and 1"

        flat = []
        for row in samples:
            for value in row:
                try:
                    flat.append(float(value))
                except (TypeError, ValueError):
                    continue

        if len(flat) < 2:
            return "Error: samples must contain at least two numeric values"

        sorted_samples = np.sort(np.array(flat, dtype=float))
        sample_count = len(sorted_samples)
        window_size = int(np.ceil(confidence * sample_count))

        if window_size >= sample_count:
            lower = float(sorted_samples[0])
            upper = float(sorted_samples[-1])
        else:
            width_count = sample_count - window_size + 1
            widths = sorted_samples[window_size - 1:] - sorted_samples[:width_count]
            start_index = int(np.argmin(widths))
            lower = float(sorted_samples[start_index])
            upper = float(sorted_samples[start_index + window_size - 1])

        median = float(np.quantile(sorted_samples, 0.5))

        return [
            ["Median", "Lower", "Upper"],
            [median, lower, upper]
        ]
    except Exception as e:
        return f"Error: {str(e)}"

Online Calculator

2D array of posterior sample draws (unitless).
Credible interval probability level (0 to 1).