Recursive step

Skip to first unread message

Aaron MacNeil

Jan 20, 2011, 1:42:25 PM1/20/11
I'm looking for a little expert help in speeding up a recursive step in PyMC. Namely the process error in the state-space surplus production model example on the pymc homepage:

# Lognormal distribution of P's
Pmean0 = 0.
P_0 = Lognormal('P_0', mu=Pmean0, tau=isigma2, trace=False, value=P_inits[0])
P = [P_0]

# Recursive step
for i in range(1,nyears):
    Pmean = Lambda("Pmean", lambda P=P[i-1], k=k, r=r: log(max(P+r*P*(1-P)-k*catch[i-1],0.01)))
    Pi = Lognormal('P_%i'%i, mu=Pmean, tau=isigma2, value=P_inits[i], trace=False)

I've wanted to speed this up for a while as this looping is tremendously slow, but I'm having difficulty thinking how it might be done. Ideas?

Thanks much,

Flavio Coelho

Jan 20, 2011, 2:42:40 PM1/20/11
Hi Aaron,

One possibility is to look into the multiprocessing module. As most CPUs nowadays are at least dual core, you can cut you time in (at least) half, by creating a pool of processes, turning you loop into a function which you call using map, and passing range(1,nyears) as argument. The documentation has examples that you can easily adapt to your code.

it would look somewhat like this:

from multiprocessing import Pool


def your_loop_core(i):
    ...., range(1,nyears))

Flávio Codeço Coelho

You received this message because you are subscribed to the Google Groups "PyMC" group.
To post to this group, send email to
To unsubscribe from this group, send email to
For more options, visit this group at

Chris Fonnesbeck

Jan 20, 2011, 6:52:45 PM1/20/11

Do you get any mileage out of assigning an AdaptiveMetropolis step method to P?

M.use_step_method(AdaptiveMetropolis, M.P)


M. Aaron MacNeil

Jan 20, 2011, 7:04:34 PM1/20/11
It is distinctly faster, but also distinctly bad at converging - in fact it spirals out of control. The Metropolis steps converge on the known values. This is in line with my general experience with AdaptiveMetropolis - often if something isn't converging well, a forced Metropolis will solve the problem. Why might this be?


David Huard

Jan 21, 2011, 1:36:50 PM1/21/11
On Thu, Jan 20, 2011 at 7:04 PM, M. Aaron MacNeil <> wrote:
> It is distinctly faster, but also distinctly bad at converging - in fact it spirals out of control. The Metropolis steps converge on the known values. This is in line with my general experience with AdaptiveMetropolis - often if something isn't converging well, a forced Metropolis will solve the problem. Why might this be?

The AdaptiveMetropolis sampler performance hinges on a lot of
parameters being set correctly. We tried to make it as simple to use
as possible by setting a lot of parameters by defaults, but those
defaults regularly lead to convergence problems. Maybe one solution
would be to never assign AdaptiveMetropolis as the default sampler and
remove its default parameters. Users would need to work out by
themselves what parameter settings work best.

More fundamentally, the main issue with AM is finding the jump
covariance matrix. This covariance is computed from the last N steps,
hence the "adaptive". If for one reason or another this covariance is
problematic, chances are the next samples won't be good, leading to an
equally problematic update of the covariance. Conclusion is, you need
to start the AM sampler in the right track and make sure it stays in
the right track all along. I don't know if its possible to find a set
of parameters that provide those conditions for all problems. Maybe

One thing that is already implemented and could help is using a past
trace to compute the initial covariance matrix. The strategy is to use
Metropolis on individual parameters for a while, then switch to AM to
speed things up, using the first trace to get a good initial
covariance. Maybe the AM sampler should do this by default.


Seth Flaxman

Jan 22, 2011, 3:51:47 AM1/22/11
Hi David. 

I believe I've been experiencing similar problems, and the approach you describe--running Metropolis for awhile, then switching to AdaptiveMetropolis--has helped in some cases, so I definitely recommend that to others as a starting point.

I'm wondering if you or anyone else can provide some more guidance / references on setting up a covariance matrix (not even sure where to start on this) and more simply determining whether the AM sampler stays on the "right track," as you put it. One approach I tried is setting the verbose flag, so I can keep on an eye on the acceptance rate, but I don't have a good sense of what to expect or hope for.

Thanks for any pointers,

Roban Hultman Kramer

Jan 24, 2011, 8:34:05 AM1/24/11

I don't have any insight into setting up initial covariance matrices,
but I've written some code to help graphically diagnose
AdaptiveMetropolis runs that you might find useful.

It's the plot_AM_correlation function in this module:

Direct link to the module file:

Direct link to example script:

Running the example script will produce a few plots (on screen). The
last plot is an example of plot_AM_correlation, showing the
correlations between variables calculated from the trace of the
Adaptive Metropolis proposal covariance matrix. Attached is an example

What you want to see: The variances (square root of the variance is
shown by the thick lines in the top panel) and covariances (or rather
correlations, shown by the other lines in the top panel) should look
like they converged to their final values early on. Any large
oscillations indicate unstable behavior in the algorithm. For
instance, if you have set shrink_if_necessary=True, and you continue
to see large drops in the variance, it indicates that few jumps are
getting accepted for some reason, and the algorithm keeps shrinking
the variance to try to increase the acceptance ratio.

Let me know if you have trouble running the example or applying it to
your own code, and feel free to email plots for help interpreting the



David Huard

Jan 24, 2011, 9:26:58 AM1/24/11
Hi Roban,

This looks like a neat tool. Thanks for sharing it.


Certainly others can chime in on this, but the acceptability of the
acceptance rate will depend on the dimension of the space you are
sampling from. In 1D, we usually target a 30%, 40% acceptance rate. As
you increase the dimensionality of the problem, that figure will
usually become lower. In my experience, problems occur when there are
less than ~3% jumps accepted during one update iteration of the

As for the initial covariance matrix, certainly you can make a guess
about the expected value of the parameters you are sampling and a
rough estimate of their variance. Use the AM method `cov_from_scales`
to transform these guesses into a covariance matrix.

If users continue to share their experience about what works and what
doesn't, I'm hopeful we'll eventually find a strategy that works most
of the time and can be set as the default.


Chris Fonnesbeck

Jan 24, 2011, 3:01:56 PM1/24/11
to PyMC
The competence value of AdaptiveMetropolis is set to zero in the
current development branch, so PyMC should not be selecting this step
method automatically for any model; it needs to be chosen manually.

It would certainly be useful for this sort of feedback to be formally
incorporated into the StepMethod. It shouldn't be too hard, for
example, to automate the switch from Metropolis to AdaptiveMetropolis
during the course of an MCMC run. The key will be coming up with
intelligent rules about when to pull the trigger.

In the meantime, to help Amac get some results (which was the initial
question in this thread), I am going to try to work up the surplus
example using the TWalk algorithm, which should improve the
performance significantly. Its in the sandbox now, and just needs to
be set up to work with containers properly.


A. Flaxman

Jan 28, 2011, 1:49:09 PM1/28/11
Hi All,

I was inspired by this discussion to do some experiments (and make some movies!) about when AM works and don't work. You can see the results here:



Thomas Wiecki

Jan 28, 2011, 2:07:15 PM1/28/11
Hi Abie,

these look great and are really helpful! It would also be interesting
to read _why_ these samplers have these problems and show this type of
behavior in the scenarios you give. But maybe that's beyond the scope
of your post and requires more background knowledge on Metropolis.


Reply all
Reply to author
0 new messages