Highest Density Interval from posterior samples

hdi.py

from getdist.mcsamples import MCSamples
import getdist.chains

def highest_density_interval_from_samples(xsamples, xlo=None, xhi=None, probability_level=0.68):
    """
    Compute the highest density interval (HDI) from posterior samples.

    Parameters
    ----------
    xsamples : array_like
        The posterior samples from which to compute the HDI.
    xlo : float or None, optional
        Lower boundary limiting the space. Default is None.
    xhi : float or None, optional
        Upper boundary limiting the space. Default is None.
    probability_level : float, optional
        The desired probability level for the HDI. Default is 0.68.

    Returns
    -------
    tuple
        A tuple containing the maximum a posteriori (MAP) estimate and the lower and upper
        bounds of the HDI.

    Notes
    -----
    The function starts at the highest density point and accumulates neighboring points
    until the specified probability level is reached. If `xlo` or `xhi` is provided,
    the HDI is constrained within these bounds.

    Examples
    --------
    >>> xsamples = np.random.normal(loc=0, scale=1, size=100000)
    >>> hdi = highest_density_interval_from_samples(xsamples)
    >>> print('x = %.1f + %.2f - %.2f' % hdi)
    x = 0.0 + 1.02 - 0.96
    """
    getdist.chains.print_load_details = False
    samples = MCSamples(samples=xsamples, names=['x'],
                           settings = {'mult_bias_correction_order':1}, 
                           ranges={'x':[xlo,xhi]})
    density_bounded = samples.get1DDensityGridData('x')
    x = density_bounded.x
    y = density_bounded.P
    
    # Sort the y values in descending order
    sorted_indices = np.argsort(y)[::-1]
    sorted_y = y[sorted_indices] / np.sum(y)
    sorted_x = x[sorted_indices]
    
    # Initialize variables
    total_probability = sorted_y[0]
    map = sorted_x[0]
    x_lo = map
    x_hi = map
    for i in range(1, len(sorted_y)):
        #print(total_probability, sorted_x[i])
        # Add the current probability to the total
        total_probability += sorted_y[i]
        x_lo = min(x_lo, sorted_x[i])
        x_hi = max(x_hi, sorted_x[i])
        # Check if the total probability exceeds or equals the desired level
        if total_probability >= probability_level:
            break
    
    return map, map - x_lo, x_hi - map