Skip to content

Rationalize the ACT running average computation in the custom dynesty sampler.

Will Farr requested to merge will-farr/bilby:nmcmc into master

Apologies for the long comment, but I wanted to explain what I'm doing here (a complete list of changes follows this comment).

The goal of the custom dynesty random walk proposal in bilby is to allow the number of MCMC steps taken in the random walk that replaces a "retiring" live point to adapt to the efficiency / autocorrelation length of the sampling. The existing implementation looks a lot like a method I have implemented before (in my Julia Enesemble sampling package), but has several differences that (a) render the ACL estimates unstable / corrupt and (b) dramatically reduce the efficiency of the sampler. This merge request attempts to "fix" the computation.

I have not yet tested against existing "GW" runs (e.g. GW150914), but I have tested a handful of my own synthetic injections and the sampler seems to work well. If someone can point me at the test cases for PE on past events, I'd be happy to run them and report back / adjust the parameters if necessary to make the sampler robust in these situations.

The updated computation is as follows:

For each retired live point, the MCMC runs for a fixed number of steps, based on the running average of the previous ACL estimates times a factor (currently set to default to 1 ACL; see below). When the MCMC is complete, the ACL for that MCMC is estimated using ACL = 2/P_accept - 1 and accumulated into the running average to set the next MCMC run length.

The quantity 2/P-1 is the autocorrelation length of an ideal, independent draw rejection sampler (i.e. Kombine) with acceptance rate P so it is an underestimate of the ACL of an MCMC that produces correlated samples, but hopefully not by too much. With this approximation, running for N ACLs will produce, on average, 2N accepted MCMC steps; the probability, therefore, of no accepted steps in a live point replacement (so a repeated sample in the live point collection) is ~exp(-2N) ~ 10% for nacl = 1. (Setting nacl = 1 does not bias the algorithm. It turns out that the common assumption that nested sampling requires a statistically independent draw for the live point updates is not correct; see http://arxiv.org/abs/1805.03924 . The only requirement is that the draws are, asymptotically, from the correct distribution of the prior truncated by the likelihood constraint.)

The ACL is averaged using a running-mean, exponential-decay algorithm with decay time constant specified by the new adapt_tscale argument (default: 100 live point replacements). This seems a good compromise between having smooth changes in the number of MCMC steps for each live point replacement (which prefers a longer averaging time constant) and allowing the ACL average to adapt to changing conditions / shapes of the constrained prior (which would prefer a smaller time constant). With adapt_tscale = 100 and nlive = 1000, the ACL over a timescale that corresponds to ~10% of the remaining prior volume.

Now that nc is stable and evolves smoothly throughout the simulation, and now that the number of MCMC steps is not changing during each live point replacement (as in the old code!), simulations seem to work well with nacl = 1. (I think before one needed nacl = 5 and draws randomly from the MCMC chain instead of the end-state in order to somewhat compensate from the bias that is introduced by changing the number of MCMC draws based on the current chain state and the rapid and random fluctuations in the ACL estimate from the buggy formula. In other words, the newly-stabilized "control system" for the number of MCMC walks and the ACL estimate makes everything more efficient.)

Detailed changes follow:

  • No longer issue a warning when there are repeated likelihood values in the sampler output. (Since these happen regularly now with nacl = 1.)

  • Add adapt_tscale argument to the custom dynesty rwalk: the running average ACT estimate will average over this number of livepoint replacements using an exponentially-decaying weight. Default is 100.

  • Get aggressive with nact = 1 by default now that the ACT computation is stabilized. Each MCMC live point replacement will run for one estimated ACT. This means, on average, two accepted MCMC steps, so a ~10% (= exp(-2)) chance of zero acceptances and a repeated live point.
    This does not bias the sampler (see http://arxiv.org/abs/1805.03924).

  • The ACT estimate is an exponentially-weighted running mean of 2/P_accept - 1 over each live point replacement; 2/P_accept - 1 is the ACT of an ideal rejection sampler with acceptance rate P_accept, so is a lower-bound for the ACT of the MCMC.

  • Always return the last point reached in the MCMC (do not randomly sample a point from earlier, as this just increases the amount of correlation with the initial point).

  • No longer update the ACL estimate while sampling is occurring; this can lead to significant bias (in practice, not just theoretically!).

  • Documentation updates describing all these.

Merge request reports