Implementing tension statistics

README.rst

@@ -2,7 +2,7 @@
 anesthetic: nested sampling post-processing
 ===========================================
 :Authors: Will Handley and Lukas Hergt
-:Version: 2.8.6
+:Version: 2.9.0
 :Homepage: https://github.com/handley-lab/anesthetic
 :Documentation: http://anesthetic.readthedocs.io/
 

anesthetic/_version.py

@@ -1 +1 @@
-__version__ = '2.8.6'
+__version__ = '2.9.0'

anesthetic/read/csv.py

@@ -5,7 +5,7 @@
 
 
 def read_csv(filename, *args, **kwargs):
-    """Read a CSV file into a :class:`Samples` object."""
+    """Read a CSV file into a :class:`anesthetic.samples.Samples` object."""
     filename = Path(filename)
     kwargs['label'] = kwargs.get('label', filename.stem)
     wldf = wl_read_csv(filename.with_suffix('.csv'))

anesthetic/tension.py

@@ -0,0 +1,90 @@
+"""Tension statistics between two datasets."""
+from anesthetic.samples import Samples
+from scipy.stats import chi2
+
+
+def tension_stats(A, B, AB, nsamples=None, beta=None):  # noqa: D301
+    """Compute tension statistics between two samples.
+
+    Using nested sampling we can compute:
+
+    - ``logZ``: Bayesian evidence
+
+      .. math::
+          \\log Z = \\log \\int L \\pi d\\theta
+
+    - ``D_KL``: Kullback--Leibler divergence
+
+      .. math::
+          D_\\mathrm{KL} = \\int P \\log(P / \\pi) d\\theta
+
+    - ``logL_P``: posterior averaged log-likelihood
+
+      .. math::
+          \\langle\\log L\\rangle_P = \\int P \\log L d\\theta
+
+    - ``d_G``: Gaussian model dimensionality
+      (or posterior variance of the log-likelihood)
+
+      .. math::
+          d_\\mathrm{G}/2 = \\mathrm{var}(\\log L)_P
+
+    - ``p``: p-value for the tension between two samples
+
+      .. math::
+          p = \\int_{d_\\mathrm{G} - 2 \\log S}^{\\infty} \\chi^2 (x)dx
+
+    Parameters
+    ----------
+    A : :class:`anesthetic.samples.Samples` or \
+        :class:`anesthetic.samples.NestedSamples`
+        (Nested) Samples from a sampling run using only dataset A.
+
+    B : :class:`anesthetic.samples.Samples` or \
+        :class:`anesthetic.samples.NestedSamples`
+        (Nested) Samples from a sampling run using only dataset B.
+
+    AB : :class:`anesthetic.samples.Samples` or \
+         :class:`anesthetic.samples.NestedSamples`
+        (Nested) Samples from a sampling run using datasets A and B jointly.
+
+    nsamples : int, optional
+        - If nsamples is not supplied, calculate mean value
+        - If nsamples is integer, draw nsamples from the distribution of
+          values inferred by nested sampling
+
+    beta : float, array-like, optional
+        inverse temperature(s) beta=1/kT. Default 1
+
+    Returns
+    -------
+    samples_stats : :class:`anesthetic.samples.Samples`
+        DataFrame containing the following tension statistics:
+        logZ, D_KL, logL_P, d_G, p
+    """
+    statsA = A.stats(nsamples=nsamples, beta=beta)
+    statsB = B.stats(nsamples=nsamples, beta=beta)
+    statsAB = AB.stats(nsamples=nsamples, beta=beta)
+
+    samples_stats = Samples(index=statsA.index)
+
+    logR = statsAB.logZ-statsA.logZ-statsB.logZ
+    samples_stats['logR'] = logR
+    samples_stats.set_label('logR', r'$\log{R}$')
+
+    logI = statsA.D_KL + statsB.D_KL - statsAB.D_KL
+    samples_stats['logI'] = logI
+    samples_stats.set_label('logI', r'$\log{I}$')
+
+    logS = logR-logI
+    samples_stats['logS'] = logS
+    samples_stats.set_label('logS', r'$\log{S}$')
+
+    d_G = statsA.d_G + statsB.d_G - statsAB.d_G
+    samples_stats['d_G'] = d_G
+    samples_stats.set_label('d_G', r'$d_\mathrm{G}$')
+
+    p = chi2.sf(d_G-2*logS, d_G)
+    samples_stats['p'] = p
+    samples_stats.set_label('p', r'$p$')
+    return samples_stats

docs/source/anesthetic.read.rst

@@ -25,6 +25,15 @@ anesthetic.read.cobaya module
    :show-inheritance:
 
 
+anesthetic.read.csv module
+--------------------------
+
+.. automodule:: anesthetic.read.csv
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+
 anesthetic.read.getdist module
 ------------------------------
 

docs/source/anesthetic.rst

@@ -72,6 +72,7 @@ anesthetic.samples module
 .. automodule:: anesthetic.samples
    :members:
    :undoc-members:
+   :show-inheritance:
 
 
 anesthetic.scripts module
@@ -83,6 +84,15 @@ anesthetic.scripts module
    :show-inheritance:
 
 
+anesthetic.tension module
+~~~~~~~~~~~~~~~~~~~~~~~~~
+
+.. automodule:: anesthetic.tension
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+
 anesthetic.testing module
 ~~~~~~~~~~~~~~~~~~~~~~~~~
 

tests/test_tension.py

@@ -0,0 +1,109 @@
+from anesthetic.examples.perfect_ns import correlated_gaussian
+import numpy as np
+from numpy.linalg import inv, solve, slogdet
+from anesthetic.tension import tension_stats
+
+
+def test_tension_stats_compatiable_gaussian():
+    d = 3
+    V = 6.0
+    nlive = 1000
+    bounds = [[-1, 1], [0, 3], [0, 1]]
+
+    meanA = [0.1, 0.3, 0.5]
+    covA = np.array([[.01, 0.009, 0],
+                     [0.009, .01, 0],
+                     [0, 0, 0.1]])*0.01
+    logLmaxA = 0
+    samplesA = correlated_gaussian(nlive, meanA, covA, bounds, logLmaxA)
+
+    meanB = [0.1, 0.3, 0.5]
+    covB = np.array([[.01, -0.009, 0.01],
+                     [-0.009, .01, -0.001],
+                     [0.01, -0.001, 0.1]])*0.01
+    logLmaxB = 0
+    samplesB = correlated_gaussian(nlive, meanB, covB, bounds, logLmaxB)
+
+    covAB = inv(inv(covA) + inv(covB))
+    meanAB = covAB@(solve(covA, meanA)+solve(covB, meanB))
+    dmeanAB = np.array(meanA)-np.array(meanB)
+    logLmaxAB = -1/2 * dmeanAB@solve(covA+covB, dmeanAB) + logLmaxA + logLmaxB
+    samplesAB = correlated_gaussian(nlive, meanAB, covAB, bounds, logLmaxAB)
+
+    nsamples = 1000
+    beta = 1
+    samples_stats = tension_stats(samplesA, samplesB, samplesAB,
+                                  nsamples, beta)
+
+    logR_std = samples_stats.logR.std()
+    logR_mean = samples_stats.logR.mean()
+    logR_exact = (-1/2*dmeanAB@solve(covA+covB, dmeanAB) -
+                  1/2*slogdet(2*np.pi*(covA+covB))[1] + np.log(V))
+    assert logR_mean - 3 * logR_std < logR_exact < logR_mean + 3 * logR_std
+
+    logS_std = samples_stats.logS.std()
+    logS_mean = samples_stats.logS.mean()
+    logS_exact = d/2 - 1/2*dmeanAB@solve(covA+covB, dmeanAB)
+    assert logS_mean - 3 * logS_std < logS_exact < logS_mean + 3 * logS_std
+
+    logI_std = samples_stats.logI.std()
+    logI_mean = samples_stats.logI.mean()
+    logI_exact = -d/2 - 1/2 * slogdet(2*np.pi*(covA+covB))[1] + np.log(V)
+    assert logI_mean - 3 * logI_std < logI_exact < logI_mean + 3 * logI_std
+
+    assert samples_stats.get_labels().tolist() == ([r'$\log{R}$',
+                                                    r'$\log{I}$', r'$\log{S}$',
+                                                    r'$d_\mathrm{G}$', r'$p$'])
+
+
+def test_tension_stats_incompatiable_gaussian():
+    np.random.rand(0)
+    d = 3
+    V = 6.0
+    nlive = 1000
+    bounds = [[-1, 1], [0, 3], [0, 1]]
+
+    meanA = [0.1, 0.3, 0.5]
+    covA = np.array([[.01, 0.009, 0],
+                     [0.009, .01, 0],
+                     [0, 0, 0.1]])*0.01
+    logLmaxA = 0
+    samplesA = correlated_gaussian(nlive, meanA, covA, bounds, logLmaxA)
+
+    meanB = [0.15, 0.25, 0.45]
+    covB = np.array([[.01, -0.009, 0.01],
+                     [-0.009, .01, -0.001],
+                     [0.01, -0.001, 0.1]])*0.01
+    logLmaxB = 0
+    samplesB = correlated_gaussian(nlive, meanB, covB, bounds, logLmaxB)
+
+    covAB = inv(inv(covA) + inv(covB))
+    meanAB = covAB@(solve(covA, meanA)+solve(covB, meanB))
+    dmeanAB = np.array(meanA)-np.array(meanB)
+    logLmaxAB = -1/2 * dmeanAB@solve(covA+covB, dmeanAB) + logLmaxA + logLmaxB
+    samplesAB = correlated_gaussian(nlive, meanAB, covAB, bounds, logLmaxAB)
+
+    nsamples = 1000
+    beta = 1
+    samples_stats = tension_stats(samplesA, samplesB, samplesAB,
+                                  nsamples, beta)
+
+    logR_std = samples_stats.logR.std()
+    logR_mean = samples_stats.logR.mean()
+    logR_exact = (-1/2*dmeanAB@solve(covA+covB, dmeanAB) -
+                  1/2*slogdet(2*np.pi*(covA+covB))[1] + np.log(V))
+    assert logR_mean - 3 * logR_std < logR_exact < logR_mean + 3 * logR_std
+
+    logS_std = samples_stats.logS.std()
+    logS_mean = samples_stats.logS.mean()
+    logS_exact = d/2 - 1/2*dmeanAB@solve(covA+covB, dmeanAB)
+    assert logS_mean - 3 * logS_std < logS_exact < logS_mean + 3 * logS_std
+
+    logI_std = samples_stats.logI.std()
+    logI_mean = samples_stats.logI.mean()
+    logI_exact = -d/2 - 1/2 * slogdet(2*np.pi*(covA+covB))[1] + np.log(V)
+    assert logI_mean - 3 * logI_std < logI_exact < logI_mean + 3 * logI_std
+
+    assert samples_stats.get_labels().tolist() == ([r'$\log{R}$',
+                                                    r'$\log{I}$', r'$\log{S}$',
+                                                    r'$d_\mathrm{G}$', r'$p$'])