Announcing multichain_mcmc package for advanced mcmc algorithms and the addition of gradient information to PyMC

[John Salvatier]

Dear PyMC List,

multichain_mcmc (http://pypi.python.org/pypi/multichain_mcmc) is a new Python package build around a branch of PyMC; it contains some advanced samplers as well as a framework for building multiple-chain samplers.

The package contains two samplers: a variant of DREAM_ZS and an AMALAsampler (Adaptive Metropolis adjusted Langevin algorithm).

Langevin algorithms use gradient information about the posterior distribution to focus sampling in the direction of increasing probability which is especially useful for high dimensional problems. AMaLa samplers also involve less manual tuning than a lot of other algorithms (more information on AMALA samplers is available from the sampler doc string). In order  to use the AMALA sampler you must install the gradientBranch PyMC branch.

The AMaLa sampler relies on a branch of PyMC (http://github.com/pymc-devs/pymc/tree/gradientBranch) which adds support for calculating the gradient of the log posterior as well as adding Deterministic-factories for many numpy functions (like sin, cos, exp, log, sum), all capable of being used with gradient calculations. In the future, PyMC should support both gradient calculation and hessian calculation so that Stochastic Newton samplers are possible.

I do not have a tutorial for using the AMALA sampler, but I hope to have one soon. The primary difference between writing models for PyMC samplers and writing models for the multichain_mcmc samplers is that for multichain_mcmc, you must write a zero parameter function which returns your model instead of giving it a list or module. This is so that the sampler can initialize multiple chains.

Several examples are included in the source code of multichain_mcmc. The most convenient place to read the examples is http://github.com/jsalvatier/multichain_mcmc/tree/master/multichain_mcmc/examples/ .

Please note that multichain_mcmc does not (yet) support multithreading even though it supports multiple chains. I have used multiple chains primarily to facilitate adaptation and convergence monitoring. Package name change suggestions and help supporting multiple chains are both welcomed.



Please let me know of any difficulty with either the PyMC gradient branch or the multichain_mcmc package, and I will try assist. Please also let me know if you have comments, questions or suggestions.

Best Regards,
John Salvatier

[Delphine Pessoa]

Good afternoon,

I hope this is not a stupid question, but I tried to do setup.py install, which worked, but when I try to import multichain_mcmc, I get this error:

/home/dpessoa/multichain_mcmc-0.3/multichain_mcmc/multichain.py in <module>()

      8 

      9 from numpy import *

---> 10 from  rand_no_replace import *

     11 from convergence import GRConvergence

     12 import time

ImportError: No module named rand_no_replace

I tried looking for it to install it manually, but I didn't find it anywhere.

What can I do?

I'm running python 2.6.4 on Ubuntu 9.10, pymc 2.1beta, numpy 1.3.0, scipy 0.7.0.

Best regards,

Delphine

--
You received this message because you are subscribed to the Google Groups "PyMC" group.
To post to this group, send email to py...@googlegroups.com.
To unsubscribe from this group, send email to pymc+uns...@googlegroups.com.
For more options, visit this group at http://groups.google.com/group/pymc?hl=en.

[John Salvatier]

Hi Delphine,

Could you send me the output text from the failed install?

John

[John Salvatier]

Hi Nor,

Yes, you are correct D needs a function to compute it's own partial derivative. I don't have an example, unfortunately, but I can describe it here:

If you can express what you want through a series of operator overloaded deterministics and built in stochastics, this is the easiest way to do this. PyMC has a lot of determinstic capabilities, so saying a + b where a an b are variables gives you a new deterministic. This includes indexing into arrays and lots of numpy-like functions sum, sin, exp, log etc. (all in the pymc namespace). If you have trouble with this, post your code and I will try to help.

If you cannot do this in this way, I have not implemented this through the decorator syntax explicitly (thought it might work), so you might need to build a stochastic from Stochastic(). In addition to the normal Stochastic() arguments (check the __doc__ string) you will need to supply the logp_partial_gradients argument with a dictionary of functions that calculate the partial gradients with respect to each of the stochastic's parameters. For example for the lognormal distribution the dictionary is :

lognormal_grad_like = {'value' : flib.lognormal_gradx,
                       'mu' : flib.lognormal_gradmu,
                       'tau' : flib.lognormal_gradtau}

The functions take the same arguments as the stochastic and return an array in the shape of the parameter for which they take the gradient with the gradients for each of the members of the input array. Let me know if this is not clear, and I will try to make it more clear.

Thanks for letting me know about your difficulty, I should get to work on improving the documentation for this.

Best Regards,
John

On Tue, Aug 31, 2010 at 8:46 AM, Nor <nor.p...@gmail.com> wrote:
>
> I have managed to install both the pymc gradient branchand

multichain.
The examples seem to run just fine and I am trying to use this for my
MCMC problem which involves a large number of dimensions. Normal MCMC
does not seem to mix right and usually gets stuck and appears to
converge away from the true minima. I was hoping that the use of the
gradient would help.
I am however getting the following error and I am not sure how to get
around it:

NotImplementedError: <pymc.PyMCObjects.Stochastic 'D' at 0x105e54fd0>
has no gradient function for parameter z

Is there something that object D should provide? This is actually a
function that computes a model and compares it to some observations
and returns a logp probability. Does it need t have a way to compute
its own partial derivative? Is this documented somewhere? I am not
sure how to proceed.

N

> most convenient place to read the examples ishttp://github.com/jsalvatier/multichain_mcmc/tree/master/multichain_m....

[John Salvatier]

I have been looking at a number of ways to implement higher order derivatives in the most functional way. I have looked a little bit at pycppad, and a few other cpp AD packages. There are two things that make them unappealing to me. 1) The C++ ones require boost, which I understand is a bit of a pain, and will make it difficult to redistribute. 2) I do not think that most of the AD packages handle operations of numpy arrays efficiently in the sense that the jacobians for UFuncs on numpy arrays have a special pattern because of broadcasting and because most UFuncs are non-interacting (the input to a Ufunc in one location of an array only affects the output in the same location).

John

On Tue, Aug 31, 2010 at 10:02 AM, A. Flaxman <ab...@uw.edu> wrote:

I’ve been meaning to take a look at your multichain code for a while now, and haven’t found the time.  I didn’t realize that it needed gradients.  Have you seen pycppad?  It is an automatic differentiation package for python, written by my colleague at UW Brad Bell http://www.seanet.com/~bradbell/pycppad/index.xml.  Unfortunately it doesn’t work with many PyMC stochastics by default, because it can’t differentiate functions that call fortran library functions.  But it could be very handy for your methods if you were to reimplement the common stochastics as pure python functions.

 

--Abie

 

 

Abraham D. Flaxman

Assistant Professor

Institute for Health Metrics and Evaluation | University of Washington

2301 5th Avenue, Suite 600 | Seattle, WA 98121| USA

Tel: +1-206-897-2800 | Fax: +1-206-897-2899 UW

ab...@uw.edu | http://healthmetricsandevaluation.org | http://healthyalgorithms.wordpress.com

--

[John Salvatier]

Algopy and Theano both look like they might suit my needs, but I want to learn tensor analysis before I continue.

[John Salvatier]

I haven't tested this kind of thing, but I would think it would do a decent job unless you have some 'stiff' regions. Keep in mind that you may have the problem that doing such a calculation is slow if you have to iterate over all of your parameters. Let me know how it goes, I have considered adding an automatic numerical differentiation for the case when derivatives are not supplied.

On Tue, Aug 31, 2010 at 10:33 AM, Nor Pirzkal <npir...@me.com> wrote:

Hi

Thanks.

My @observed function is rather non trivial and uses some external code to generate models, given the input parameters that pymc is trying to sample. Would I be able to get by with a simple numerical computation of the derivative using something like (D(z)-D(z+dz))/dz?

N

[Delphine Pessoa]

Hello John Salvatier,

I am writing a master thesis and would like to reference you package. How should I do?

Thank you,

Delphine Pessoa

On 21 June 2010 19:11, John Salvatier <jsal...@u.washington.edu> wrote:

[John Salvatier]

That is a good question, and I am not sure what the correct answer is. Perhaps the pypi page (http://pypi.python.org/pypi/multichain_mcmc/)?

Does anyone have any advice?

I am curious what your thesis is about.

John

[Delphine Pessoa]

I'm working on estimating some parameters for pneumococcus transmission in healthy children using bayesian inference.

Referenced the site as a web page using zotero, results in this bibliography entry:

[Salvatier] Salvatier, J.  Python package index : multichain mcmc 0.3. http://pypi.python.org/pypi/multichain mcmc.

Maybe a year would be better. Is it 2010?

Delphine

[John Salvatier]

Sounds neat!

The most recent version was uploaded in 2010, so I think that's probably the best year.

Don't forget the underscore _ in multichain_mcmc.

John

[John Salvatier]

You may also wish to link to the specific version you used, http://pypi.python.org/pypi/multichain_mcmc/0.3 , instead of the most recent version.