Pierre Jacob writes regarding convergence diagnostics for Markov chain simulation:
I’ve implemented an example of TV upper bounds for (vanilla) HMC on a model written in Stan, see here and here for a self-contained R script.
Basically, this creates a stan fit object to obtain a target’s pdf and gradient, and then implements a pure R implementation of some of the methods we’re proposing with Niloy for a simple, non-dynamic HMC algorithm. The model is a logistic regression with random effects. The dimension is around 500, coupled HMC chains seem to contract quickly, with meetings occurring exactly in typically less than 1000 iterations. From that, we obtain the proposed TV upper bounds. Of course the target is very simple, I wanted the whole script to run in a few minutes.
Let me also advertise that my paper with John O’Leary and Yves Atchadé on that topic will be a read paper at JRSS B, read in December.
Comments can be sent with a deadline two weeks after the discussion meeting on December 11, that is, December 28. If you know of people who might be interested in writing a comment on that paper, please spread the word.
In passing, I’ve also noted a typo in this page in the leapfrog integrator section, if rho is Normal(0, Sigma), then the second step of the integrator should read theta <- theta + epsilon Sigma^{-1} rho. Also, I believe that "mass matrix" usually refers to the covariance of the momentum variable in the literature, e.g. "M" in "MCMC using Hamiltonian dynamics" by Radford Neal, which is contrary to what you write in the section "Auxiliary Momentum Variable".
New research and found a typo: good news all around.
P.S. Relatedly, here’s our recent paper on convergence diagnostics.