diff --git a/bilby/core/sampler/dynesty_utils.py b/bilby/core/sampler/dynesty_utils.py index 6d9ee5f7c..3096bb97d 100644 --- a/bilby/core/sampler/dynesty_utils.py +++ b/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, diff --git a/test/core/sampler/dynesty_test.py b/test/core/sampler/dynesty_test.py index ca15eb021..edf25e864 100644 --- a/test/core/sampler/dynesty_test.py +++ b/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"],