16
$\begingroup$

I am new to Bayesian Statistics research. I heard from researchers that Bayesian researchers better implement MCMC by themselves rather than using tools like JAGS/Stan. May I ask what is the benefit of implementing MCMC algorithm by oneself (in a "not quite fast" languages like R), except for learning purpose?

$\endgroup$
3
  • $\begingroup$ Because then you can chose your own proposal distribution yourself, you should chose it such that the Markov Chain resulting from it converges as fast as possible to the posterior. $\endgroup$
    – user83346
    Sep 14, 2016 at 4:21
  • 5
    $\begingroup$ If you're an applied researcher who wants to learn more Bayes by using it for applications, I'd recommend starting with JAGS or Stan and then moving to writing your own MCMC if you find you 'need' to. Bear in mind that JAGS and Stan have slightly different strengths and limitations. $\endgroup$ Sep 14, 2016 at 12:38
  • $\begingroup$ thanks! Yes, I am doing applied research. Would you tell me more about the limitations of JAGS and Stan. I first tried Stan, but I just found it does not have "online monitoring" or "sample until converge" features or add-ons --- this is annoying, I may try JAGS now. $\endgroup$
    – user112758
    Sep 14, 2016 at 22:42

3 Answers 3

29
$\begingroup$

In general, I would strongly suggest not coding your own MCMC for a real applied Bayesian analysis. This is both a good deal of work and time and very likely to introduce bugs in the code. Blackbox samplers, such as Stan, already use very sophisticated samplers. Trust me, you will not code a sampler of this caliber just for one analysis!

There are special cases in which in this will not be sufficient. For example, if you needed to do an analysis in real time (i.e. computer decision based on incoming data), these programs would not be a good idea. This is because Stan requires compiling C++ code, which may take considerably more time than just running an already prepared sampler for relatively simple models. In that case, you may want to write your own code. In addition, I believe there are special cases where packages like Stan do very poorly, such as Non-Gaussian state-space models (full disclosure: I believe Stan does poorly in this case, but do not know). In that case, it may be worth it to implement a custom MCMC. But this is the exception, not the rule!

To be quite honest, I think most researchers who write samplers for a single analysis (and this does happen, I have seen it) do so because they like to write their own samplers. At the very least, I can say that I fall under that category (i.e. I'm disappointed that writing my own sampler is not the best way to do things).

Also, while it does not make sense to write your own sampler for a single analysis, it can make a lot of sense to write your own code for a class of analyses. Being that JAGs, Stan, etc. are black-box samplers, you can always make things faster by specializing for a given model, although the amount of improvement is model dependent. But writing an extremely efficient sampler from the ground up is maybe 10-1,000 hours of work, depending on experience, model complexity etc. If you're doing research in Bayesian methods or writing statistical software, that's fine; it's your job. But if your boss says "Hey, you can you analyze this repeated measures data set?" and you spend 250 hours writing an efficient sampler, your boss is likely to be upset. In contrast, you could have written this model in Stan in, say, 2 hours, and had 2 minutes of run time instead of the 1 minute run time achieved by the efficient sampler.

$\endgroup$
6
  • 3
    $\begingroup$ +1. Also, Stan does not directly handle some problems involving discrete distributions, so you have to know enough to integrate these out which is not in itself simple, so that might be a case where rolling your own might help. I believe JAGS handles such cases directly, though, so if you can keep the differing philosophies of BUGS/JAGS and Stan separated in your mind, it would be best to just switch between them. $\endgroup$
    – Wayne
    Sep 14, 2016 at 14:49
  • $\begingroup$ Moreover, Stan can have problems where a diagonal Euclidean metric is not well-suited to the geometry of the posterior; this is the case inter alia when there is only a narrow, oddly-shaped region of the posterior which has much probability. The result is that sampling the posterior is like trying to ride a bike along the edge of a cliff: you might "fall off" if you take a wrong turn! $\endgroup$
    – Sycorax
    Sep 14, 2016 at 15:05
  • 2
    $\begingroup$ +1. My general recommendation to students is to code it up in JAGS. If that doesn't work well, then code it up in Stan. If that doesn't work well, then start writing your own sampler. There are also certain models, e.g. spatial models, where you may want to use BUGS. And certain models, e.g. non-Gaussian state-space models, where you want to use NIMBLE. The opportunity cost of starting by writing your own sampler is just too high. $\endgroup$
    – jaradniemi
    Sep 14, 2016 at 15:42
  • $\begingroup$ I don't understand the "real time" case - if it's possible to have an "already prepared" own sampler why isn't it as easy to use an already compiled Stan model? I also wonder if any MCMC is fast enough for real time applications. $\endgroup$ Sep 14, 2016 at 17:00
  • 1
    $\begingroup$ And I'm not familiar enough with Stan to know what exactly requires compiling new models, but it's not too hard to imagine that whatever the restriction is, there exists a dynamic model such that as new data comes in, the model would become more complex and so recompiling would be necessary. I think non-parametric methods (in which the parameter space grows with the sample size) would fit that criteria? But there may be clever ways of getting around that. $\endgroup$
    – Cliff AB
    Sep 14, 2016 at 17:14
6
$\begingroup$

This question is primarily opinion based, but I think there is enough here write an answer down. There could be many reasons to code one own's sampler for a research problem. Here are some of them

  1. Proposal: As fcop suggested in their comment, if the sample is M-H, then coding your own sampler lets you play around with proposal distributions to get the best mixing sampler.

  2. Flexibility: In built programs might not give you the flexibility you want. You might want to start at a specific random value, or use a specific seed structure.

  3. Understanding: Coding your own sampler helps you understand the behavior of the sampler, giving insights to the Markov chain process. This is useful for a researcher working on the problem.

  4. Onus: If the data on which I am making all my Bayesian inference comes from a program that I didn't code up, then the onus on the inference is no longer on me. As a researcher, I would like to take full responsibility of the methods/results I present. Using in-built methods does not allow you to do that.

There are probably more reasons, but these are the four that make me code my own samplers.

$\endgroup$
4
  • 7
    $\begingroup$ I'd say that the "trust" reason is disputable: Stan is open-source and has lots of contributors, so multiple persons have looked at it's source code and so it is unlikely that it has serious bugs. On another hand, if you do it by yourself you can always overlook the bug that you made - and everybody makes bugs it is just a matter of number of lines of code you write... $\endgroup$
    – Tim
    Sep 14, 2016 at 12:23
  • $\begingroup$ @Tim I agree. I have changed that point to reflect what I was trying to say. Thanks. $\endgroup$ Sep 14, 2016 at 12:31
  • 5
    $\begingroup$ +1 for the Understanding argument. However, the Onus argument seems a little bit overstated. Almost anything you code up yourself will rely on somebody else's statistical language, linear algebra library, random number generator, etc. so 'taking responsibility' is a matter of degree. $\endgroup$ Sep 14, 2016 at 12:43
  • $\begingroup$ @conjugateprior Absolutely agreed. Which is why my answer on that was in the first person. This was purely my opinion. $\endgroup$ Sep 14, 2016 at 14:55
4
$\begingroup$

I gave a +1 to Cliff AB's answer. To add one little tidbit, if you want to work at a lower level but not down to the code-everything-yourself level, you should poke around for the LaplacesDemon package. The original author was brilliant, but seems to have dropped off the grid, and the package has been taken over by someone else. (It's on Github, I believe.)

It implements an impressive number of algorithms used in MCMC and the included vignettes are worth the read even if you don't use the package. Pretty much any kind of sampler you read about, it has. You code in a different way than BUGS/JAGS or Stan, and it's all in R, but often times it's so efficient that it's competitive.

$\endgroup$
5
  • 1
    $\begingroup$ Shameless plug: you can also use [nimble](r-nimble.org) that allows you to customize your MCMC (i.e. use a slice sampler for this node, block updater for that group of nodes, etc.) without the need to rewrite this samplers every time. And you can also write your own samplers to be directly implemented! Disclosure: I used to work on this project. $\endgroup$
    – Cliff AB
    Sep 14, 2016 at 15:27
  • $\begingroup$ @CliffAB: Sounds similar to LaplacesDemon, if you're familiar with that. Glad to hear about nimble as well. I'll at least download it. (Though the multiple LaplacesDemon vignettes might be worth the download even if you use nimble.) ... Ohhh, just went to the page. If its SMC is easy to use, I'll become a big fan. The only R package I've seen that does SMC is horrifically complex. $\endgroup$
    – Wayne
    Sep 14, 2016 at 16:00
  • $\begingroup$ @CliffAB: Wow, after reading the nimble website, it's pretty impressive. Why have I never heard of it? It looks like a great option for people used to the BUGS/JAGS modeling language. Of course, they'll make the best-possible comparisons on the website, but still I like it so far. (Except that with rstanarm and brms, which use Stan under the hood, the ease-of-use-in-R champion would be Stan.) $\endgroup$
    – Wayne
    Sep 14, 2016 at 16:10
  • $\begingroup$ it's still very new: v0.1 was released I think just over 2 years ago? And SMC was a huge motivator for the project: the PI has done considerable publishing on particle filters and was getting annoyed with writing them up from scratch every time. But I've been a bit with current work to keep up with what the current state of the SMC samplers are; when I left (almost two years ago) we had just put together a very primitive one. $\endgroup$
    – Cliff AB
    Sep 14, 2016 at 16:44
  • 1
    $\begingroup$ Just say this arXiv paper you might be interested in. $\endgroup$
    – Cliff AB
    May 24, 2017 at 16:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.