Implementing the pCN algortihm

Development
1

Dear all

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

2 Likes
2

Hi Mikkel,
Thanks again for your contribution of the new sampler to PyMC3.
As for the pCN algorithm, I think you can take a look at pymc3.smc, where it isolate the likelihood term alone: https://github.com/pymc-devs/pymc3/blob/ee57bed9aa5f73dc48081dd3c82267dd4f0f2834/pymc3/smc/smc.py#L127

My recommendation is to follow the pattern in SMC and program it into a specialized sampler, you can initialized a class method that contains the state of the last likelihood.

1 Like
3

Hi @junpenglao. Thank you for the pointer. It does indeed seem like the likelihood function in the SMC sampler makes use of the model.datalogpt() method to compute the likehood, so I will have a go with that.

But as to your other suggestion to build it around the pattern in SMC, I think I might have explained the problem poorly. The pCN algortihm is an approach for MCMC sampling. I did not know that the Metropolis algorithm could also relate to SMC. Anyway I was planning to make the pCN algorithm as an MCMC step method, so I think I will use the Metropolis step method as a mold, if that makes sense.