Implementation of Adaptive Metropolis Algorithm by Haario et al., + storing rejected points

[Murillo Vellasco]

Hi everyone!

I've been using emcee for some time, and for the purposes of our work I had to do some modifications to emcee and I would like to ask here if these would be useful and would fit the emcee vision before actually submitting a pull request. 

I am running scans on a very high-dimensional parameter space (SMEFT), and the relevant scale for each parameter can very significantly. We therefore decided to use the Adaptive Metropolis Algorithm by Haario et al., and subsequently implemented it in emcee. In this way, the parameters of the Gaussian MH move can be "auto-tuned" to fit the scales of the parameter space.

Because the Adaptive Metropolis Algorithm (AMA) is not technically MCMC (since the hyperparameters of the proposal density evolve over time), the implementation involved some nontrivial modifications. In particular, in the definition of the State class, I had to include the covariance matrix and mean of the sample at the current iteration to the State attributes. The updating of the algorithm hyperparameters is performed in "emsemble.py", such that the user can have the option to update the hyperparameters using the entire sample or only those points which were generated using an AMA move.

Another modification I have implemented is the ability to save all points proposed by emcee. In our work, the computation time for the likelihood function for any given point in our parameter space is of the order of 10 seconds, such that our statistics are quite limited. Since we are taking a frequentist approach, where we eventually evaluate the profile likelihood using the sample obtained with emcee, we don't actually need the sample to be distributed according to the likelihood function. In the end, we only need the scanning algorithms to sufficiently (and efficiently) explore the parameter space. Thus, in principle, we don't actually need to discard the Burn-in points. Also, by having access to all points proposed by emcee, including both those accepted and rejected, we can in principle increase our statistics by a factor between 2 and 5, assuming an acceptance fraction between 0.2 and 0.5. Given the very large computation times for our likelihood function, such a modification has been quite valuable, especially given that, for our purposes, we need to estimate the full 5\sigma likelihood intervals/contours, not just the usual 95% confidence intervals, and most of the points outside of the typical 1-2\sigma regions are discarded by the MCMC algorithms (this is also why we decided to implement the adaptive algorithm, instead of just using just the Stretch Move).

The way we have performed this modification was to change the "propose" functions of the moves such that they return an additional state, which is the state that the chain would be in if all proposed points had been accepted. We then defined "full" versions of the "chain", "log_prob" and "blobs" in the backend to store these. We can access these simply by setting a "full" option to True in the "get_sample", "get_log_prob" and "get_blobs" functions.

I hope these modifications can be useful to someone else!

Murillo