initial working draft of kameleon mcmc

[karlnapf]

Hi

as discussed a while ago, I wrote a draft for the kernel adaptive metropolis hastings sampler from
http://jmlr.org/proceedings/papers/v32/sejdinovic14.html
based on our codes from
https://github.com/karlnapf/kameleon-mcmc

This implementation is based on oracle samples. The actual implementation needs to use the chain history. I did not implement this yet in order to discuss general points with you first, but soon to follow.

There are several things that are not as clean as they could

• It would be nice to make this sampler a subclass of a metropolis Hastings one with just the proposal being different (and scaling parameters)
• One problem is that the proposal here is not symmetric, so it needs to be incldued in the MH ratio, see Pseudo Marginal MCMC & MH step #592
• In particular, the scale tuning code is duplicate, I dont know however whether you want this to allow for different schedules
• The helper class for Gaussian would be better put somewhere else, see helper class for Gaussians #593
• I wrote the kernel function using theano (using its autodiff for the gradient, which keeps to code a bit cleaner), should this go to some helper class? Though the Gaussian kernel is the one for using this in practice
• I have a seperate class on my computer where I plot the proposal in every nth step (if the target is 2d), not sure whether this belongs into pymc though, thoughts?

Finally, no unit tests yet. If this stuff is supposed to go on and we finalised the structure, I am happy to extensively unit test everything, for now I can put an example?

Let me know any comments and I'll update the code.

[karlnapf]

I am not sure whether you allow this, but compiling my expressions fails otherwise.

[twiecki]

Very exciting, thanks for putting this up. I'll try to take a look soon but probably won't be able to in the next week or two (vacation after defending ;) ). Regarding the non-symmetric proposals, #547 has an implementation of MH that might help here.

[karlnapf]

Thanks!

I am more than happy to use the MH step in #547 if it is exposed in a way that I can inherit from the class. If not, the current design is super simple: I just store the log-likelihood and the proposal representation of an accepted step.

This patch is not meant to be merged yet, more for figuring out how to embed it within pymc. I need to add a few things on using the chain history. But that's a second step. First I am interested how to organise the code, you will see there is lots of extra stuff that has to be put somewhere.

[karlnapf]

travis fails for an unrelated reason I think

[twiecki]

I thought you only need a custom proposal dist? https://github.com/pymc-devs/pymc/pull/547/files#diff-3c696bd218b69d877a316c792c5efa70R31 seems to suggest it's pretty general or why wouldn't that work?

[karlnapf]

The main use-case of this sampler is when no gradient is available. For example: pseudo-marginal mcmc (aka noisy target density evaluation). In this case, HMC (which is better than the kernel-MH) is not an option, but only a purely random walk.

So to use this thing, one will usually define a custom density which does the pseudo-marginal stuff. The implementation I did also takes care of the accept/reject in a correct way. If one wants to do a std random walk as pseudo-marginal, you need to change your metropolis step method as it re-computes the target density in every step (this is a bad idea in general as it might be expensive) -- one should store the likelihood of the last accepted state.

[fonnesbeck]

sample already has a tune argument to take care of this.

[karlnapf]

Thanks for the comments @fonnesbeck
I changed the tuning internals to work with your sample method. Was not aware of that before.

[fonnesbeck]

This looks really good. Thanks for the contribution!

Though the paper does not seem to say so explicitly, I assume this method is only applicable to continuous random variables?

[fonnesbeck]

Feel free to add an example model to the examples folder also.

[karlnapf]

• Yeah only continuous for now (we are working on extensions)
• where to put the Gaussian sampling?
• I will add an example
• Do you do any unit testing? ideas on that?

[fonnesbeck]

We have a set of unit tests in pymc/tests, in which you will see tests for the other step methods.

As for the Gaussian sampling, I started a branch to implement random sampling across all distributions, but have not worked on it in awhile. This is ideally where we would like it, as we will have more general usage for random sampling than just this particular step method.

[twiecki]

@karlnapf Any plans to continue here or should we close this?

[karlnapf]

No, planning to continue.
In fact taking up work on this again from this week, as ICML deadline (on another topic) is over and I have resources to get involved here again. So expect an update within a week or so.

[twiecki]

@karlnapf Good to hear! Good luck with ICML.

[jsalvatier]

@karlnapf Hey Karl! I hope ICML went well! Any news with this? It would be cool to see it production.

[karlnapf]

Thanks for the reminder! :)
Yes, ICML is submitted, NIPS is submitted. So now there is more time. I actually worked on an extension for this, but will polish this first and send it.

[jsalvatier]

Lovely! I look forward to it

On Mon, Jun 8, 2015 at 3:17 AM, Heiko Strathmann [email protected]
wrote:

Thanks for the reminder! :)
Yes, ICML is submitted, NIPS is submitted. So now there is more time. I
actually worked on an extension for this, but will polish this first and
send it.

—
Reply to this email directly or view it on GitHub
#599 (comment).

[springcoil]

Looking forward to this too.

[springcoil]

Any update on this @karlnapf? I think it'd be a great merge. Maybe we can split this up to a few commits...

[karlnapf]

I re-wrote large part of the internals for a recent NIPS paper and another working paper. I will try to put the things in here when I am back next week as I actually got some time at my hand now after the ICML deadline is over.

[twiecki]

Closing for now, feel free to reopen if you make progress.