FEAT: discard dynesty mcmc chains that are too short

bilby/core/sampler/dynesty_utils.py

@@ -4,6 +4,7 @@
 import numpy as np
 from dynesty.internal_samplers import InternalSampler, SamplerReturn
 from dynesty.utils import SamplerHistoryItem, apply_reflect, get_random_generator
+from scipy.stats import binom
 
 from ...bilby_mcmc.chain import calculate_tau
 from ..utils.log import logger
@@ -109,6 +110,8 @@ def __init__(self, **kwargs):
         self.sampler_kwargs["walks"] = self.walks
         self.naccept = kwargs.get("naccept", 10)
         self.maxmcmc = kwargs.get("maxmcmc", 5000)
+        self.sampler_kwargs["target"] = self.naccept
+        self.sampler_kwargs["maxmcmc"] = self.maxmcmc
 
     def tune(self, tuning_info, update=True):
         """
@@ -228,6 +231,36 @@ def sample(args):
             "reject": walks - naccept,
         }
 
+        # check that the number of accepted steps is consistent with the target
+        # number of accepted steps, if it is much lower than the target, our
+        # mcmc chain is too short and we may be open to bias. If we're hitting
+        # maxmcmc, this is expected, but it also happens at the beginning of the
+        # run or around phase transitions, so we discard this point and so an
+        # equivalent longer mcmc will happen at the next iteration
+        dist = binom(n=walks, p=args.kwargs["target"] / walks)
+        minimum = dist.ppf(1e-4)
+        undersampled = naccept < minimum
+        if undersampled and walks >= args.kwargs["maxmcmc"]:
+            logger.debug(
+                "Requested %s accepted steps per MCMC chain, "
+                "but only %s were accepted out of the maximum allowed number "
+                "of steps (%s). Increase 'maxmcmc' to improve convergence.",
+                args.kwargs["target"],
+                naccept,
+                walks,
+            )
+        elif undersampled:
+            logger.info(
+                "Only %s accepted steps out of target %s. "
+                "Discarding this point. "
+                "This is expected at the beginning of the run, increase 'walks'. "
+                "This can also happen when the geometry of the likelihood surface "
+                "changes drastically.",
+                naccept,
+                args.kwargs["target"],
+            )
+            logl = args.loglstar - 1
+
         return SamplerReturn(
             u=current_u,
             v=current_v,

test/core/sampler/dynesty_test.py

@@ -410,6 +410,7 @@ def test_proposal_functions_run_new(self):
         args = Dummy(u=np.ones(4) / 2, axes=np.ones((2, 2)))
         args.live_points[0] += 1
         args.kwargs["ncdim"] = 2
+        args.kwargs["target"] = 1
         # args.kwargs["live"][0] += 1
         for proposals in [
             ["diff"],