Example Stan converges but BUGS does not

[jt...@st-andrews.ac.uk]

Hi, 

Are there any standard examples in Stan that does not work in BUGS (ie MCMC does not converge or take a miraculously long time to converge)?

[Stephen Martin]

I've *never* seen a case where Stan takes longer than BUGS *in terms of effective samples*.
BUGS/JAGS may take much less time to cycle through 2000 iterations, and it may even say it converges.

However, and this is a *big* however:
1) BUGS/JAGS do not use random inits. This means the chains are more likely to converge to one mode of an inherently multimodal posterior.
Due to the problems of gibbs, random walk metropolis, and slice sampling, they will likely stay in the one mode.
The convergence diagnostics will say "hey, the gibbs sampler has converged across all chains", when in reality the problem is that gibbs found one mode and failed to escape it due to to the conditional sampling that it uses.
This means you can incorrectly conclude that you explored the full posterior, when in reality you only explored one mode of it, and it failed to escape the mode to find other high density areas of the posterior. Your inference will be potentially jeopardized as a result.

2) Stan will take longer than Gibbs. Period. Stan has to evaluate the gradient of the posterior in order to simulate a particle flowing across the multidimensional surface.
You may say "well, my problem is unimodal, so why not use Gibbs, since it's faster"?
The issue is that per-iteration, Gibbs is faster, but in terms of uncorrelated, representative samples, Stan is usually faster.
Gibbs has less to compute per iteration, so it can crawl through 2000 iterations quickly.
Stan has more to compute per iteration, so it will crawl through 2000 iterations slower.
Then you look at the output and see that Gibbs may only have 40 effective samples (relatively independent samples), because due to conditional sampling, Gibbs can get stuck in highly correlated areas of the posterior.
On the other hand, Stan may have 1800 effective samples, because it samples across the posterior in a joint, unconditional manner, naturally flowing across the surface. It doesn't have much of an issue with correlated areas, and won't get 'stuck' as often as Gibbs.

As a concrete example of this, I currently have a very high dimensional model.
It includes a series of mixture models and latent variables.
It takes Gibbs about 20 minutes to go through it, and Stan about 40 minutes to go through it (2000 iterations each, 1000 of them are warmup in my case).
After 20 minutes, Gibbs has about 30 effective samples (yes, you read that correctly, 30 out of 1000 post-adaptation draws are relatively independent).
After 40 minutes, Stan has about 900.
I'll take that 40 minutes: 900/40 >>> 30/20; roughly 15x 'faster' in terms of effective samples from the posterior.

[Andrew Gelman]

Hi, Bugs and Jags, like Stan, can take discrete inits.  But I think they also can do random inits.  It depends on the setting.  Also, Gibbs is not always faster; it depends on the distribution and on what is available closed form.

But I agree on the general point of effective sample size.  Also it can be truly liberating to not have to think at all about conjugacy.

A

--
You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
To post to this group, send email to stan-...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Stephen Martin]

Right, I shouldn't imply that bugs/jags *can't* take inits. They can, but by default, iirc, they do not randomize their inits (I think Bob looked this up for me when I first started using Stan early this year).

What *is* an example of when stan is actually faster, per-iteration, than gibbs? Seems like gibbs, proper, would be faster, but bugs/jags don't always use gibbs, and I'm guessing their other samplers may be slower.

Didn't even think about the lack of conjugacy. This has been a huge benefit for me. Trying to parameterize a complicated model in terms of unintuitive distributions and unintuitive parameters just to satisfy jags was a hassle (e.g., constantly converting between precision and sd was irritating, to say the least). That on top of non-vectorization meant my model in jags was a bunch of spaghetti code.

[Bob Carpenter]

The punchline of (2) is why we built Stan in the first place.

I also wouldn't trust your estimate of a 30 effective sample
size (30 effective draws) if you had 1000 iterations --- that's too
much long-range autocorrelation for the estimators to be very
reliable. (Also, Stan uses slightly more conservative estimates
than JAGS does in that it uses split R-hat and then uses that
to further adjust effective sample size based on R-hat-like properties;
it's described in the manual.)

HMC/NUTS is almost always slower than Gibbs per iteration because
we almost always evaluate more than a single log density and
gradient and the gradients are done with autodiff.
As Andrew was alluding to, Gibbs can get expensive
when you don't have conjugate forms in the conditionals because
you have to start doing things like slice sampling (JAGS) or
adaptive rejective sampling (BUGS). But in everything we've
looked at, JAGS has been faster than Stan per iteration.

A lot of this comes down to implementation details, too.
JAGS can do some neat tricks with GLMs, Stan can do some
neat tricks with vectorization. So all this relative speed
stuff is very model specific. And can even depend on the
data---Stan's not very good with badly misspecified models.
(Goes back to Andrew's folk theorem.)

- Bob

> On Nov 20, 2016, at 9:39 PM, Stephen Martin <hwki...@gmail.com> wrote:
>

[jt...@st-andrews.ac.uk]

Thank you for that answer. I certainly got a lot of insight from that. Perhaps, I should be more specific to my question. Is there a standard example in the Stan Documentation or the Case Studies that show  an example that Stan would be able to but BUGS cannot. For instance, I know Gibbs Sampling performs poorly when your parameters are correlated such as the x and y variable inside a multivariate normal. I am aiming to re-create or duplicate a simple  example where BUGS cannot solve the problem but STAN can, especially in this case for correlated samples. 

[jt...@st-andrews.ac.uk]

Thanks

[Andrew Gelman]

Hi, we have some examples in the Nuts paper:

http://www.stat.columbia.edu/~gelman/research/published/nuts.pdf

I assume that if you implemented these, they'd fail in Bugs but run fine in Stan.

A

[Bob Carpenter]

The paper example was just a highly correlated multivariate normal.
But you can't code it that way in BUGS/JAGS or you'd just be doing forward
simulation and it'd be fast.

You need to find something in high dimensions where the parameters
are highly correlated in the posterior.

- Bob

[Andrew Gelman]

The paper had several examples!

[Michael Betancourt]

Unfortunately this is a somewhat ill-posed question.

MCMC will always get the right answer asymptotically, but we might

not be able to reach that asymptotic regime in practice.  In other

words, our MCMC estimators will converge to the right answer as

the number of iterations goes to infinity, but they will only become

useful if they converge fast enough.

Consequently there are two stages of questions to ask:

1 - Will my sampler be able to explore the target distribution well

enough to be able to reach this asymptotic regime in a finite

amount of time?  Formally the sufficient and usually necessary

condition for this to happen is called geometric ergodicity.

2 - Assuming these nice enough properties, will the convergence

time be fast enough that I can get sufficiently precise estimators

in practice?

Most people focus on (2), but (1) is an extremely important condition

in practice.  The problem is that (1) is really, really hard to verify in theory,

and still tricky to verify in practice.  We know that HMC is geometrically

ergodic over a larger class of models  than many other algorithms, but

it’s hard to construct a nontrivial example of something that is unequivocally 

fine for HMC but unequivocally bad for another algorithm.

So if you want to do a proper comparison you need to follow these steps:

1 - Simulate data from a model with known parameters.

2 - Ensure that all diagnosis are good, including general diagnostics

like Rhat and sampler-specific diagnostics like those we have for

HMC in Stan.

3 - Ensure that your MCMC algorithm recovers the exact parameters within

a few MCMC standard errors.  

Any samplers passing these conditions have not demonstrated any

failures of (1), although that is not a guarantee that (1) holds so be careful.

Now we can continue on to focus on (2).

4 - Compute effective sample sizes.  If effective sample size / iteration < 1e-3

then your sampler probably isn’t converging fast enough to even quantify

its precision well enough to compare.

5 - Finally, rank the samplers by effective sample size / time.  You can

set a target effective sample size and maximum time so that if

(effective sample size / time) * target effect sample size > maximum time

then the sampler is considered impractical 

Any samplers failing (5) can be considered to fail (2).

Anyways, the reason that we don’t have a bunch of examples that we can

quote off the top of our head is that doing these comparisons correctly is

subtle and requires lots of time and we haven’t been able to test everything.

You can see some hierarchical models that I used a few years ago in

https://arxiv.org/abs/1312.0906, although even those tests weren’t as

careful as we now know that they should have been!

[Andrew Gelman]

Examples do come up all the time, though.  I haven't run Bugs/Jags recently, but back when I did, I found those programs to solve lots of problems but to choke sometimes too, particularly when multivariate distributions were involved. And we've had some discussions on the users list where people report models they can't fit in those other programs but do fit in Stan.  I guess what I'm saying is that, yes, everything Mike says is correct.  And, also, sometimes Bugs/Jags are just so clearly off (and so clearly slow) for problems that Stan can solve.

And, just to emphasize, this is no knock on Bugs/Jags, which do pretty good considering they're based on 20-yr-old algorithms.

A

[Michael Betancourt]

Yup — I really should have added  “tl;dr” emphasizing

that when a model fits poorly, it’s really hard to determine

whether it was just too slow given the finite computational

resources, or it never would have fit in a finite time at all.

So there are lots of examples out there that don’t fit in

BUGS/JAGS that we don’t bother analyzing seriously

because we have Stan, which is not only super fast but

also robust and features all kinds of novel diagnostics

that we can use to provide evidence for geometric

ergodicity in practice.

[Bob Carpenter]

But only one was relevant to this question of comparing
NUTS to Gibbs (Fig 7 and section 4.5):

http://www.jmlr.org/papers/volume15/hoffman14a/hoffman14a.pdf

Everything else is comparing NUTS to static HMC.

- Bob

[Andrew Gelman]

Sure, but I'm guessing that Gibbs and Metropolis wouldn't've fared well in these examples!