Implemented credibility_interval()

anesthetic/samples.py

@@ -8,6 +8,8 @@
 import pandas
 import copy
 import warnings
+from scipy.interpolate import interp1d
+from scipy.optimize import minimize_scalar
 from collections.abc import Sequence
 from anesthetic.plot import (make_1d_axes, make_2d_axes, fastkde_plot_1d,
                              kde_plot_1d, hist_plot_1d, scatter_plot_2d,

anesthetic/utils.py

@@ -3,6 +3,7 @@
 import pandas
 from scipy import special
 from scipy.interpolate import interp1d
+from scipy.optimize import minimize_scalar
 from scipy.stats import kstwobign
 from matplotlib.tri import Triangulation
 import contextlib
@@ -85,6 +86,69 @@ def quantile(a, q, w=None, interpolation='linear'):
     return quant
 
 
+def credibility_interval(samples, weights=None, level=0.68, method="hpd"):
+    """Compute the credibility interval of weighted samples. Based on
+    linear interpolation of the cumulative density function, thus expected
+    discretization errors on the scale of distances between samples.
+
+    Parameters
+    ----------
+    samples: array
+        Samples to compute the credibility interval of.
+    samples: array, optional
+        Weights corresponding to samples. Default: ones
+    level: float, optional
+        Credibility level (probability, <1). Default: 0.68
+    method: str, optional
+        Which definition of interval to use. Default: 'hpd'
+        * hpd: Calculate highest posterior density (HPD) interval,
+          equivalent to iso-pdf interval (interval with same probability
+          density at each end) if distribution is unimodal.
+        * ll, ul: Lower limit, upper limit. One-sided limits for which
+          `level` fraction of the samples lie above / below the limit.
+        * et: Equal-tailed interval with the same fraction of samples
+          below and above the interval region.
+
+    Returns
+    -------
+    limit: tuple or float
+        Tuple [lower, upper] for hpd/et, or float for ll/ul
+    """
+    if level >= 1:
+        raise ValueError('level must be <1, got {0:.2f}'.format(level))
+    if len(np.shape(samples)) != 1:
+        raise ValueError('Support only 1D arrays for samples')
+    if weights is not None and np.shape(samples) != np.shape(weights):
+        raise ValueError('Shape of samples and weights differs')
+
+    weights = np.ones(len(samples)) if weights is None else weights
+    # Sort and normalize
+    order = np.argsort(samples)
+    samples = np.array(samples)[order]
+    weights = np.array(weights)[order]/np.sum(weights)
+    # Compute inverse cumulative distribution function
+    S = np.array([np.min(samples), *samples, np.max(samples)])
+    CDF = np.append(np.insert(np.cumsum(weights), 0, 0), 1)
+    invCDF = interp1d(CDF, S)
+
+    if method == "hpd":
+        # Find smallest interval
+        def distance(Y, level=level):
+            return invCDF(Y+level)-invCDF(Y)
+        res = minimize_scalar(distance, bounds=(0, 1-level), method="Bounded")
+        return np.array([invCDF(res.x), invCDF(res.x+level)])
+    elif method == "ll":
+        # Get value from which we reach the desired level
+        return invCDF(1-level)
+    elif method == "ul":
+        # Get value to which we reach the desired level
+        return invCDF(level)
+    elif method == "et":
+        return np.array([invCDF((1-level)/2), invCDF((1+level)/2)])
+    else:
+        raise ValueError("Method '{0:}' unknown".format(method))
+
+
 def check_bounds(d, xmin=None, xmax=None):
     """Check if we need to apply strict bounds."""
     if len(d) > 0:

tests/test_samples.py

@@ -8,6 +8,7 @@
 from anesthetic import MCMCSamples, NestedSamples, make_1d_axes, make_2d_axes
 from anesthetic.samples import merge_nested_samples
 from anesthetic.samples import merge_samples_weighted
+from anesthetic.utils import credibility_interval
 from numpy.testing import (assert_array_equal, assert_array_almost_equal,
                            assert_array_less)
 from pandas.testing import assert_frame_equal
@@ -887,6 +888,23 @@ def test_plotting_with_integer_names():
         samples_1['0']
 
 
+def test_credibility_interval():
+    np.random.seed(3)
+    samples = NestedSamples(root='./tests/example_data/pc')
+    assert np.allclose(credibility_interval(samples["x0"], level=0.68,
+                       weights=samples.weights, method="hpd"),
+                       [-0.1, 0.1], atol=0.02)
+    assert np.allclose(credibility_interval(samples["x0"], level=0.95,
+                       weights=samples.weights, method="et"),
+                       [-0.2, 0.2], atol=0.02)
+    assert np.isclose(credibility_interval(samples["x0"], level=0.975,
+                      weights=samples.weights, method="ul"),
+                      0.2, atol=0.02)
+    assert np.isclose(credibility_interval(samples["x0"], level=0.5,
+                      weights=samples.weights, method="ll"),
+                      0, atol=0.001)
+
+
 def test_logL_list():
     np.random.seed(5)
     default = NestedSamples(root='./tests/example_data/pc')