[Draft] Functions to compute confidence intervals from MCMCSamples or NestedSamples

[Stefan-Heimersheim]

Description

Quick demo of #178, only implements this for kde and not tested for posteriors with multiple peaks. Demo works with tests/test_samples.py

ns = NestedSamples(root='./tests/example_data/mn'); ns.confidence_1d("x0", levels=[0.95,0.68], ncompress=100000)

Checklist:

• I have performed a self-review of my own code
• My code is PEP8 compliant (flake8 anesthetic tests)
• My code contains compliant docstrings (pydocstyle --convention=numpy anesthetic)
• New and existing unit tests pass locally with my changes (python -m pytest)
• I have added tests that prove my fix is effective or that my feature works

[codecov]

Codecov Report

Merging #179 (15501ea) into master (405c0ff) will decrease coverage by 1.99%.
The diff coverage is 9.52%.

@@            Coverage Diff             @@
##           master     #179      +/-   ##
==========================================
- Coverage   95.16%   93.16%   -2.00%     
==========================================
  Files          17       17              
  Lines        1758     1800      +42     
==========================================
+ Hits         1673     1677       +4     
- Misses         85      123      +38     

Impacted Files	Coverage Δ
anesthetic/samples.py	89.11% <9.52%> (-10.89%)	⬇️

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 405c0ff...15501ea. Read the comment docs.

[Stefan-Heimersheim]

Closing this in favor of the much neater version discussed here.

import scipy.interpolate as sip
def fastCL(samples, weights, level=0.68):
    # Sort and normalize
    order = np.argsort(samples)
    samples = samples[order]
    weights = weights[order]/np.sum(weights)
    # Compute inverse cumulative distribution function
    CDF = np.append(np.insert(np.cumsum(weights), 0, 0), 1)
    S = np.array([-np.min(samples), *samples, np.max(samples)])
    #cdf = sip.interp1d(S, CDF)
    invcdf = sip.interp1d(CDF, S)
    # Find smallest interval
    distance = lambda a, level=level: invcdf(a+level)-invcdf(a)
    a = sop.minimize(distance, (1-level)/2, bounds=[(0,1-level)], method="Nelder-Mead").x[0]
    interval = np.array([invcdf(a), invcdf(a+level)])
    return interval

fastCL(samples, weights)