MCMC in Python: PyMC Step Methods and their pitfalls

There has been some interesting traffic on the PyMC mailing list lately. It seems that there is a common trouble with the “Adaptive Metropolis” step method, and it’s failure to converge. I’ve been quite impressed with this approach, and I haven’t had the problems that others reported, so I started wondering: Have I been lucky? Have I been not looking carefully?

I decided to do some experiments to make Metropolis and Adaptive Metropolis shine, and since I’m an aspiring math filmmaker these days, I made my experiments into movies.

I consider the Metropolis step method the essence of MCMC. You have a particular point in parameter space, and you tentatively perturb it to a new random point, and then decide to accept or reject this new point with a carefully designed probability, so that the stationary distribution of the Markov chain is something you’re interested in. It’s magic, but it’s good, simple magic.

Here is the simplest example I could come up with, sampling from a uniform distribution on the unit square in the plane using Metropolis steps. Any self-respecting step method can sample from the unit square in two dimensions!

It is worth mentioning that this model takes nearly no code in PyMC (making the movie takes a little bit of work, though). Setting up and running the model is just this:

X = mc.Uniform('X', lower=-1, upper=1, value=[0,0])
mod = mc.MCMC({'X': X})
mod.use_step_method(mc.Metropolis, mod.X)
mod.sample(5000)

But the Metropolis approach doesn’t fare so well when the parameters are correlated (to speak in statistician’s terms). Geometrically speaking, Metropolis is bad when the sample space is, for example, a diagonal strip, instead of an axis-aligned square:

This is the setting where Adaptive Metropolis really wins big. By updating its proposal distribution after a bit of exploring, Adaptive Metropolis effectively converts this long, narrow diagonal region back into a square. You can see how the autocorrelation drops after the first batch of steps are used for tuning:

Tuning the proposal distribution is not always a good idea, however. In this X-shaped region, the Metropolis stepper isn’t great, but Adaptive Metropolis is tricked really badly. It thinks it has discovered that it lives in a diagonal, and so it gets stuck when it starts to explore the other diagonal in the X.

The situation is much more complicated in the models that arise in real statistical modeling applications, but I think this is a useful sketch of why Metropolis is sometimes better and sometimes worse than Adaptive Metropolis.

What do you do if you need to sample something analogous to the X? I offer here a “proof of concept” PyMC implementation of a very cool step method, the hit-and-run sampler. It’s a pretty simple concept, and I hope that someone can take it to the industrial strength level required to really be used in PyMC. You pick a random direction, and look forwards and backwards in that direction to figure out what the posterior probability density looks like on that line. Then you hop to some point, with probability proportional to the posterior at each point.

To be explicit,

class HitAndRun(Gibbs):
    def __init__(self, stochastic, proposal_sd=None, verbose=None):
        Metropolis.__init__(self, stochastic, proposal_sd=proposal_sd,
                            verbose=verbose, tally=False)
        self.proposal_tau = self.proposal_sd**-2.

    def step(self):
        x0 = copy(self.stochastic.value)
        dx = rnormal(zeros(len(x0)), self.proposal_tau)

        logp = [self.logp_plus_loglike]
        x_prime = [x0]

        for direction in [-1, 1]:
            for i in xrange(25):
                delta = direction*exp(.1*i)*dx
                try:
                    self.stochastic.value = x0 + delta
                    logp.append(self.logp_plus_loglike)
                    x_prime.append(x0 + delta)
                except ZeroProbability:
                    self.stochastic.value = x0

        i = rcategorical(exp(array(logp) - flib.logsum(logp)))
        self.stochastic.value = x_prime[i]

        if i == 0:
            self.rejected += 1
        else:
            self.accepted += 1

And what does this get you? A step method that cuts through the cross-diagonal region like it was the unit square: