Emulated termination criteria

README.rst

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

anesthetic/_version.py

@@ -1 +1 @@
-__version__ = '2.7.3'
+__version__ = '2.8.0'

anesthetic/samples.py

@@ -1269,6 +1269,43 @@ def truncate(self, logL=None):
         index = np.concatenate([dead_points.index, live_points.index])
         return self.loc[index].recompute()
 
+    def terminate(self, eps=1e-3, logL=None, n=None):
+        """Check if a set of samples has reached a termination criterion.
+
+        Uses the termination criterion of
+        [Handley et al. 2015](https://arxiv.org/abs/1506.00171).
+        computes if the ratio of evidence in the live to dead points is less
+        than some precision.
+
+        Parameters
+        ----------
+        eps : float, optional
+            The precision of the criteria.
+            default: 1e-3
+
+        logL : float or int, optional
+            Loglikelihood or iteration number to truncate run.
+            If not provided, truncate at the last set of dead points.
+            default: None
+
+        n : int, optional
+            Number of samples to draw when computing the volume estimates
+            default: None
+
+        """
+        logL = self.contour(logL)
+        i_live = ((self.logL >= logL) & (self.logL_birth < logL)).to_numpy()
+        i_dead = ((self.logL < logL)).to_numpy()
+        if np.any(i_dead):
+            logZ_dead = self[i_dead].recompute().logZ(n).mean()
+            logX_dead = self[i_dead].recompute().logX(n).iloc[-1]
+        else:
+            # logZ if no dead points
+            logZ_dead = -1e30
+            logX_dead = -1e30
+        logZ_live = self[i_live].recompute().logZ(n).mean() + logX_dead
+        return logZ_live - np.logaddexp(logZ_live, logZ_dead) < np.log(eps)
+
     def posterior_points(self, beta=1):
         """Get equally weighted posterior points at temperature beta."""
         return self.set_beta(beta).compress('equal')

tests/test_samples.py

@@ -1180,6 +1180,21 @@ def test_truncate(cut):
         assert_array_equal(pc, truncated_run)
 
 
+def test_terminate():
+    np.random.seed(4)
+    pc = read_chains("./tests/example_data/pc")
+
+    assert not pc.terminate(logL=0)
+
+    assert not pc.terminate(logL=200)
+    assert not pc.terminate(logL=0.0)
+    assert pc.terminate(logL=None)
+
+    assert pc.terminate(logL=200, eps=1.0)
+    assert pc.terminate(logL=0.0, eps=1.0)
+    assert pc.terminate(logL=None, eps=1.0)
+
+
 def test_hist_range_1d():
     """Test to provide a solution to #89"""
     np.random.seed(3)