Skip to content

Instantly share code, notes, and snippets.

@JohannesBuchner
Last active March 19, 2024 17:23
Show Gist options
  • Save JohannesBuchner/2027d0f313521387c2cded2424cdcfeb to your computer and use it in GitHub Desktop.
Save JohannesBuchner/2027d0f313521387c2cded2424cdcfeb to your computer and use it in GitHub Desktop.
Highest Density Interval from posterior samples
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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment