Explaining a bad chain (teaching example)

[Richard McElreath]

Here's an easy and maybe fun question for the Stan team.

Next, week I'm introducing Stan to my Bayesian stats course. So about 60 students are about to start fighting with Stan. I have a handful of examples of pathological chains that I like to introduce MCMC with. One I like is what happens with flat priors and small amounts of data. Here's an example: 

# bad chain is bad

mc_bad <- '

data{

    int N;

    real y[N];

}

parameters{

    real mu;

    real<lower=0> sigma;

}

model{

    y ~ normal( mu , sigma );

}

'

library(rstan)

m_bad <- stan( model_code=mc_bad , data=list(y=c(-1,1),N=2) )

Trace plot for m_bad is comical. I show the class how to fix it up with weak priors, like mu~normal(0,10); sigma~cauchy(0,1);.

Now here's my question. I think I know what's going on during sampling in such chains. But I don't know enough about NUTS to explain it in a technically correct way, with mass and momentum and all the other cool things. I'm not sure I *want* to give that explanation to the class, but I do want to know it.

So would one of the NUTS wizards mind explaining the erratic sampling behavior in the m_bad chains?

[Andrew Gelman]

Hi, I'm ccing stan-dev because one thing we've been talking about, for the benefit of researchers, is to have an option for Stan to return more of its internals, details regarding what is happening within the Nuts or HMC steps.

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.
For more options, visit https://groups.google.com/groups/opt_out.

[Bob Carpenter]

You don't even need to bother with sampling to see what's
going wrong. Just look at the probabilty mass in the posterior,
which is proportional to the log probability function by Bayes's
rule. I believe in this case, you can calculate the posterior
analytically.

The general problem with the improper prior on sigma and mu
is that most (as in infinitely much of) of the probability mass
is away from the values you want, namely in a neighborhood
around mu = 0, sigma=1.

To get a feeling for what's going on geometrically in Stan, you need
to look at the probability function on the unconstrained space,
which log-transforms sigma (and performs the appropriate Jacobian
adjustment).

Then think of the log probability function on the unconstrained space
as a bowl and the parameter vector as a marble.
Now consider what happens when you thwack the marble in a random
direction and let it go for a fixed amount of time. That's essentially
what HMC is doing. You don't have to worry about NUTS to see the problem.

You might want to contrast the geometry with what you get by duplicating
the data points, starting with (1,1,-1,-1) and continuing. You'll see
the geometry begin to look reasonable around the neighborhood where you
want the estimate to be.

- Bob

[Bob Carpenter]

The useful bit's already there --- you'll see a crazy large step
size. Just look at the stanfit object that gets returned.

- Bob

>> stan-users+...@googlegroups.com <mailto:stan-users+...@googlegroups.com>.

[Bob Carpenter]

You don't even need to bother with sampling to see what's
going wrong. Just look at the probabilty mass in the posterior,
which is proportional to the log probability function by Bayes's
rule. I believe in this case, you can calculate the posterior
analytically.

The general problem with the improper prior on sigma and mu
is that most (as in infinitely much of) of the probability mass
is away from the values you want, namely in a neighborhood
around mu = 0, sigma=1.

To get a feeling for what's going on geometrically in Stan, you need
to look at the probability function on the unconstrained space,
which log-transforms sigma (and performs the appropriate Jacobian
adjustment).

Then think of the log probability function on the unconstrained space
as a bowl and the parameter vector as a marble.
Now consider what happens when you thwack the marble in a random
direction and let it go for a fixed amount of time. That's essentially
what HMC is doing. You don't have to worry about NUTS to see the problem.

You might want to contrast the geometry with what you get by duplicating
the data points, starting with (1,1,-1,-1) and continuing. You'll see
the geometry begin to look reasonable around the neighborhood where you
want the estimate to be.

- Bob

[Richard McElreath]

If there's nothing especially unusual about how NUTS is behaving in this case (compared to basic HMC, as explained in Neal's chapter in the Handbook), then I think I'm good. But I thought maybe the U-turn stuff might have something to do with the details of the chain's path. For example, I was concerned that the way Stan learns during warmup might be relevant.

I should definitely make a slide with an infinite hockey rink on it.

[Bob Carpenter]

We do the same warmup for HMC as for NUTS. It looks
for the largest step size it can with the right acceptance
rate. In this case, it's pretty big :-)

- Bob

>> an email to stan-users+...@googlegroups.com <javascript:>.

[Richard McElreath]

Can I bother you folks with another easy question about Stan/NUTS? This one is really about why sampling is sometimes so much slower. I first noticed that some colleagues were saying Stan was slow, but I was getting better timings. But my models all had weak priors of some kind, while colleagues were using flat priors (omitting priors).

Here's a minimal example to isolate what I'm talking about. I'll use a demonstration chain with unidentifiable parameters in it:

mc_bad2 <- '

data{

    int N;

    real y[N];

}

parameters{

    real a1;

    real a2;

    real<lower=0> sigma;

}

model{

    y ~ normal( a1+a2 , sigma );

}

'

library(rstan)

y <- rnorm(100)

m_bad3 <- stan( model_code=mc_bad2 , data=list(y=y,N=length(y)) )

So of course a1 and a2 wander around weirdly, even though their sum provides a good estimate of the mean. The above takes about 5 seconds (per chain) on my old machine. 

Now, the below chain produces practically identical inferences, but executes much faster: about one second per chain.

mc_bad3 <- '

data{

    int N;

    real y[N];

}

parameters{

    real a1;

    real a2;

    real<lower=0> sigma;

}

model{

    sigma ~ cauchy(0,1);

    a1 ~ normal(0,10);

    a2 ~ normal(0,10);

    y ~ normal( a1+a2 , sigma );

}

'

m_bad3 <- stan( model_code=mc_bad3 , data=list(y=y,N=length(y)) )

All I've done is add weak priors. I get why the chains wander a lot more in m_bad2 than in m_bad3. But what I don't get is why each sample takes more computer time in m_bad2.

Thanks.

[Andrew Gelman]

Hey, people are saying that Stan is slow!!!  Please have them send us some examples.

Thx

A

[Michael Betancourt]

Okay, this is a classic model where priors are hugely important.  Short story: Stan is doing just fine.  Long story: read on.

The big problem with the first model is that it is truly unidentified.  No matter how much data you have you’ll never be

able to identify a1 and a2, only their sum.  That means that the a1-a2 posterior converges to the line y = x, which you

can equally well think about as uniform on a one-dimensional space.  It may look like your chains converged here, but I

promise you it did not if only because x, y = +infinty and x, y = -infinty are both neighborhoods with just as much probability

mass as x, y = 0 and it would take quite a long time for your chains to reach that far.  You’re getting biased here because

the means look okay — but the variances are way off.  The marginal variances for a1 and a2 should be infinite!

The beauty of priors on a1 and a2  is that it suppresses that posterior mass at the infinities (the Cauchy prior helps, but it’s

less of a drastic effect in this case).  Now your chains have only a finite region to explore and Stan will cover that region

quickly (there are still strong correlations that stress the sampler, but HMC shines in these examples where as other

samplers explore in fiery debris).  There’s a huge difference between a finite number (here sigma = 10) and an infinite 

number (sigma = infinity), which is why this model provides such a nice argument for there _always_ being prior information 

to incorporate.

This kind of effect is pretty common in complex models.  In high-dimensions and/or with intricate constraints, the likelihood

mass can be warped and twisted and possibly stretch out towards infinity.  But no matter how complex your likelihood may

be it can always be made to look Gaussian if your priors are tight enough.  For example, checkout Figure 3 in

http://arxiv.org/pdf/0708.3797.pdf.

[Andrew Gelman]

Yup.  And, to elaborate a bit more, it's impossible in general to detect that the model is unidentified (i.e., has an improper posterior distribution), as this is a property of the entire density not of any local part.  This is one good motivation to use proper priors--but, more than that, to use reasonably scaled priors, not the normal(0,1000) things that people have traditionally used in Bugs (and which only work with strong data or because people choose reasonable starting points and don't let their chains run long enough to truly drift).

This will be an important thing to mention in our "Statistical Modeling Using Stan" book!

[Bob Carpenter]

If you have colleagues having trouble fitting models in Stan,
please send them to our users list --- we're happy to help.

And no bother at all. Yours are good examples for us to think
about because we have been gearing up to do some teaching
on Stan ourselves (when we create teaching materials, we'll
make them available).

For Gibbs sampling, it's easy to sample from p(a1|a2,sigma)
and from p(a2|a1,sigma). You almost certainly won't get as much drift
from the the initial values for a1 and a2 in the no-prior model.
You should get infinite drift, so Stan's not going far enough
either. There's an improper prior, so all bets are off on MCMC
working at all!

The key to understanding what's going on in Stan is to
think about the posterior geometry.

In the example with no priors, you get a ridge in the posterior
where both parameters can run off long distances in opposite directions
without changing the log probability. Now with initial momenta being
giving randomly, you won't quite run off in opposite directions, but
it'll be close because we use the gradient. The result is that you
often won't get a u-turn at all and the sampler will apply 1000
leapfrog iterations (tree depth = 10).

Tree depth is currently not printed along with the usual print() output
in R and is hard to fish out (see the stanfit-class object). So
I've just created a feature request to dump these out with the other
diagnostics the way that CmdStan does:

https://github.com/stan-dev/rstan/issues/45

You can also see how far the parameters run by inspecting the traceplots
or the posterior intervals (see below).

In the example with priors, when a1 and a2 try to run off
in opposite directions, the prior turns them around. They still get
pretty far, though, because the prior's weak. Tighter priors would
speed things up even further.

I'm working on the next edition of the manual right now and I'm going
to add a chapter on identifiability and use your example (and cite
you, of course). We need to say more on identifiability anyway, and
this is a good place to start because it's such a nice clean example.

Here are the models I used (I dropped the prior for sigma --- it doesn't
suffer from the same problem because it's identified in the posterior).

noprior.stan
-----------------------------

data {
int N;
real y[N];
}
parameters {
real a1;
real a2;
real<lower=0> sigma;
}
model {
y ~ normal(a1 + a2, sigma);
}

generated quantities {
real a;
a <- a1 + a2;
}

weakprior.stan
-------------------------------

data {
int N;
real y[N];
}
parameters {
real a1;
real a2;
real<lower=0> sigma;
}
model {

a1 ~ normal(0,10);
a2 ~ normal(0,10);
y ~ normal(a1 + a2 , sigma );
}

generated quantities {
real a;
a <- a1 + a2;
}

sim.R
--------------------------------------
library('rstan');
N <- 100;
y <- rnorm(N);

fit_weak <- stan('weakprior.stan', data=c("N","y"), seed=12345);
fit_no <- stan('noprior.stan', data=c("N","y"), seed=12345);

print(fit_weak);
print(fit_no);



Timing
----------------------------------
fit_weak: 1.0s per chain
fit_no: 6.5s per chain


Fits:

weakprior.stan
--------------------------------------------
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
a1 -0.5 0.2 7.1 -14.0 -5.2 -0.5 4.3 12.9 864 1
a2 0.5 0.2 7.1 -13.0 -4.4 0.4 5.1 14.0 864 1
sigma 1.2 0.0 0.1 1.0 1.1 1.2 1.3 1.4 1083 1
a 0.0 0.0 0.1 -0.3 -0.1 0.0 0.1 0.2 3458 1
lp__ -67.3 0.0 1.2 -70.5 -67.8 -67.0 -66.4 -65.9 818 1



noprior.stan
--------------------------------------------------
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
a1 -1888.8 483.6 938.1 -3609.5 -2522.1 -1971.5 -1327.2 290.9 4 1.6
a2 1888.8 483.6 938.1 -291.1 1327.0 1971.3 2522.1 3609.5 4 1.6
sigma 1.2 0.0 0.1 1.1 1.2 1.2 1.3 1.4 40 1.1
a 0.0 0.0 0.1 -0.3 -0.1 0.0 0.1 0.2 3213 1.0
lp__ -66.9 0.1 1.1 -69.9 -67.2 -66.6 -66.1 -65.8 59 1.1


So you see that both get reasonable Rhat and n_eff for a. weakprior.stan gets
reasonable n_eff for sigma, but norprior.stan doesn't. Also look at those intervals
on a1 and a2 in noprior; those are highly dependent on initialization.

The big difference is that weakprior.stan actually identifies a1 and a2! If
you ran a really long time, you should see them both take posterior mean of 0.
As is, you can see a very high posterior sd for their estimates and hence
a fairly high MCSE (which is estimated by posterior sd / sqrt(n_eff)).

Here's what happens at iter=100000:

Inference for Stan model: weakprior.
4 chains, each with iter=1e+05; warmup=50000; thin=1;
post-warmup draws per chain=50000, total post-warmup draws=2e+05.

mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
a1 0.0 0 7.0 -13.8 -4.8 0.0 4.7 13.8 38150 1
a2 0.0 0 7.0 -13.9 -4.8 0.0 4.7 13.8 38144 1
sigma 1.2 0 0.1 1.0 1.1 1.2 1.2 1.4 50701 1
a 0.0 0 0.1 -0.3 -0.1 0.0 0.1 0.2 177514 1
lp__ -67.3 0 1.2 -70.5 -67.9 -67.0 -66.4 -65.9 43360 1


Tighter priors would, of course, continue to improve behavior for weak_prior.

- Bob

[Richard McElreath]

Thanks, gents. I get all the unidentified stuff. Every sampler handles such a model slightly differently, so it's become a ritual of mine to try out this model under different algorithms. 

In that quest, I'm still fuzzy on why the model with flat priors takes so much longer to sample. My understanding from Bob's comments is that the leapfrog tree keeps hitting it's maximum in the flat model. The tree keeps growing, because none of the simulations turn around. As the tree grows, there's a bunch of simulations to keep track off, and that's the overhead cost that makes sampling for the flat model so much more intensive than the model with weak priors.

Does that sound right?

Re people with slow Stan models: So far I've been able to fix the problems of everyone who comes to me. I should be posting after action reports on this list, though, for the benefit of others. On that note, there is one user who is working on a censored gamma latency model, and we're having lots of numerical issues. So when I get time to actually dig into it, I likely will be bringing the problem here.

[Michael Betancourt]

>
> In that quest, I'm still fuzzy on why the model with flat priors takes so much longer to sample. My understanding from Bob's comments is that the leapfrog tree keeps hitting it's maximum in the flat model. The tree keeps growing, because none of the simulations turn around. As the tree grows, there's a bunch of simulations to keep track off, and that's the overhead cost that makes sampling for the flat model so much more intensive than the model with weak priors.
>
> Does that sound right?

Not quite. With flat priors your posterior mass extends from x, y = + infinity to x, y = - infinity. In order to draw an independent sample NUTS is going to pick a direction and chug along to one of the infinities. It would go forever but we have a the max_depth limit to avoid these issues. With weak priors there’s only a finite region to explore and NUTS turns around quickly.

Again, it’s not that Stan is slow in the model without priors but rather that the model is ill-posed and it would take infinite time to generate an independent sample no matter what you did. The overhead for building a bit tree is not trivial but it’s not that important.

[Andrew Gelman]

Right. Any method that _doesn't_ blow up in an improper example actually has a problem.

[Bob Carpenter]

On 2/6/14, 7:06 PM, Richard McElreath wrote:
> Thanks, gents. I get all the unidentified stuff. Every sampler handles such a model slightly differently, so it's become
> a ritual of mine to try out this model under different algorithms.
>
> In that quest, I'm still fuzzy on why the model with flat priors takes so much longer to sample. My understanding from
> Bob's comments is that the leapfrog tree keeps hitting it's maximum in the flat model. The tree keeps growing, because
> none of the simulations turn around. As the tree grows, there's a bunch of simulations to keep track off, and that's the
> overhead cost that makes sampling for the flat model so much more intensive than the model with weak priors.
>
> Does that sound right?

It's why it takes so long --- if you look at the negative log
probability it's a long valley, and Stan moves up and down it.
Moving down the valley involves many log probability and gradient
evals. Each update in JAGS/BUGS is conjugate and only moves
a little way. Basically, Gibbs devolves into a random walk
that barely begins to explore the improper posterior.

But, and this is the thing to keep in mind, neither Stan nor
BUGS nor JAGS is actually properly sampling from this model's
posterior, which is improper. As Michael pointed out, values
of lambda_1 = 1e100 should be as likely as those with lambda_1 = 0.

>
> Re people with slow Stan models: So far I've been able to fix the problems of everyone who comes to me. I should be
> posting after action reports on this list, though, for the benefit of others. On that note, there is one user who is
> working on a censored gamma latency model, and we're having lots of numerical issues. So when I get time to actually dig
> into it, I likely will be bringing the problem here.

Bring it on :-) Though our censoring models haven't really
been optimized at all, nor have they been tuned for numerical
stability.

- Bob

[Bob Carpenter]

Technically, Stan doesn't blow up. It hits that max tree
depth. Without a max tree depth, it'd essentially go into
an infinite loop until it ran out of memory (arguably, the
correct behavior here).

I noticed when exploring this, by the way, that the
max tree depth of 10 given the way Michael's counting allows
2000 leapfrog steps. The improper model averages 1300, which
is why each iteration is so slow.

Also, I'm running into some boundary issues in some of these
models that should be better behaved that I'm going to share
in a bit.

- Bob