My name is MIkkel Lykkegaard and I am in the same academic group as has submitted the PR for the Multilevel Delayed Acceptance sampler: https://github.com/alan-turing-institute/pymc3.
I am working to implement the preconditioned Crank Nicolson (pCN) algorithm in PyMC3: https://arxiv.org/pdf/1202.0709.pdf. It is a very simple modification of Metropolis algorithm, but it has dimension-independent acceptance ratio, and works (in theory) even for sampling problems of infinite dimension. The paper I referred to explains its advantages over Random Walk Metropolis in detail. It comes with two caveats, though:
- First, the prior must be a zero-mean Gaussian, and the pertubation vector of the proposal must come from the prior. In PyMC3 that can be achieved by setting S=cov_prior in the stepping method, where cov_prior is the covariance matrix of the prior. This would be fairly simple to guarantee in the code, too.
- Second, and this is the more crucial bit, the proposal is non-symmetric, and the acceptance ratio is simply L(x’)/L(x_i), that is the ratio between the likelihoods of the proposal and the state, and it does not include the prior probability, like standard Metropolis. In PyMC3, the acceptance ratio of Metropolis is computed using the delta_logp() function, which makes use of the model.logpt() method, which, as far as I can tell, includes the prior probability. I have found two ways to hack this “problem”: (1) using a flat prior and the pm.Potential class for the likelihood or (2) add a modified model.logpt() method to the model, which skips the var.logpt and only includes the potentials. However, I am not really satisfied with any of these “solutions”, since they seem to violate some of the design concepts of PyMC3.
So I would like to ask for some advice on how to implement the pCN algorithm most elegantly within the PyMC3 ecosystem. Would the model.datalogpt() be the correct method to call to evaluate only the likelihood? It appears to be doing that, but I can’t seem to find anything in the documentation, and the PyMC3 vocabulary is a bit different to what I am used to from MCMC theory.
Any advice would be greatly appreciated. Thanks in advance.
With kind regards, Mikkel