Frozen periods on a simple model

[Kevin Van Horn]

One of Stan's selling points is the much faster mixing (per iteration) of the Markov chain, so I was surprised to find instead that I got very slow mixing on a simple model (see attached) due to the chain freezing for an extended period. The model is basically a logistic regression with some random effects. There are only two covariates. To get the ratio of MCSE to posterior SD down to 0.05 I found I had to go out to 2,560,000 iterations (with four chains). Half of those are adaptation, and the other half are sampling.

A plot of the lp__ samples from chain 4 reveals what is happening: the chain freezes for 50 * 1280 = 64,000 iterations, with nothing changing. Am I right that this means that the update was rejected 64,000 times in a row? What can cause that to happen? This is a simple model, after all.

In this example you only see one such frozen period, but I have seen other variations of this that had three out of the four chains having such frozen periods, one chain with three of them.

[Kevin Van Horn]

Some info I left out: the run in question had N=27, K=2, M=16, and mfirst=2.

[Michael Betancourt]

This can happen when you have really bad tails -- the probability to transition out to that tail is very small, but by detailed balance that means that the probability of transitioning back is also really small. In fact, the longer you run the more likely you are to jump to a tail and get stuck for long periods of time (ironically making your autocorrelations larger). More sophisticated algorithms like RHMC can help, and using a variable step size can be enough to help you escape the tails, but the biggest improvement often comes from tightening up your prior distributions to make the tails more well-behaved.

> --
> 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.

[Kevin Van Horn]

Thanks for the explanation, but I don't think it applies here. This is not a case of seeing very little variation in lp__ or some other value -- this is a case of seeing *zero* variation. The sample values remain *identical* for 64,000 iterations. As far as I can tell, the only way this behavior can occur is if, for 64,000 times in a row, it simulates a trajectory and the numerical inaccuracies in the final result are large enough to cause the proposed update to be rejected. Is that correct, or am I missing something here?

[Michael Betancourt]

> Thanks for the explanation, but I don't think it applies here. This is not a case of seeing very little variation in lp__ or some other value -- this is a case of seeing *zero* variation. The sample values remain *identical* for 64,000 iterations. As far as I can tell, the only way this behavior can occur is if, for 64,000 times in a row, it simulates a trajectory and the numerical inaccuracies in the final result are large enough to cause the proposed update to be rejected. Is that correct, or am I missing something here?

That's pretty much the same explanation -- if the acceptance probability from A to B is vanishingly small, then the acceptance probability from B to A will also be vanishingly small. One cause could be hugely different gradients in two regions which would cause high integrator error and small acceptance probabilities. This is why I mentioned a variable step size. Again, we're talking very very small acceptance probabilities here, which is why you only see this behavior when running tens and hundreds of thousands of transitions.

[Bob Carpenter]

I don't have anything else to add on mixing theory,
but you can make the model faster by using sum() instead
of incrementing:

ll <- 0;
for (n in 1:N) {
real eta;
eta <- alpha + X[n] * beta + epsilon[month[n]];
ll <- ll + w0[n] * bernoulli_logit_log(0, eta) + w1[n] * bernoulli_logit_log(1, eta);
}

can be written as

{
reals ll[N];
for (n in 1:N) {
real eta;
eta <- alpha + X[n] * beta + epsilon[month[n]];
lls[n] <- w0[n] * bernoulli_logit_log(0, eta) + w1[n] * bernoulli_logit_log(1, eta);
}
}
ll <- sum(lls);

The brackets are important because they make ll a local variable
so that the entries aren't printed.

It doesn't seem like it'd be faster because of the extra
allocation and assignments, but it cuts down the size
of the expression tree for auto-diff, which is a speed
bottleneck for Stan.

Do you know how much eta varies? It needs to stay roughly
in the -700 to 700 range to not overflow or underflow on
bernoulli_logit_log. That can also cause things to get stuck.

- Bob

[Kevin Van Horn]

Bob,

Thanks for the tip on using sum().

The eta values are all moderate, lying between about -0.7 and 0.4. After some more investigation I found two interesting things about the frozen period:

• The epsilon[m] values are all close to zero (ranging from -0.0175 to 0.0176).
• Sigma is quite small (0.0181), barely above the smallest value it ever attains (0.0165).
  

That is, we can almost perfectly fit the data (or rather, the empirical proportion of 1 responses for each row of X) without introducing the random effects.

Have you run into situations like this before? If so, how have you handled them?

[Andrew Gelman]

What prior distribution are you using for the variance of these "random effects"?  If the data "want" the variance to be 0, but you're using a prior such as lognormal or inverse-gamma that excludes 0, you can get a problem, via the Folk Theorem of Statistical Computing.

[Kevin Van Horn]

The prior for sigma is cauchy(0, 2.5).

[Andrew Gelman]

OK, so that's not the problem!

[Kevin Van Horn]

Based on some of the comments I have gotten, I began to suspect that the problem was a stepsize that was usually OK, but was a bad choice for very small values of sigma. So I changed

    eta <- alpha_ + X[n] * beta_ + epsilon_[month[n]];
    ...
    epsilon ~ normal(0, sigma);

to

    eta <- alpha_ + X[n] * beta_ + sigma * epsilon_[month[n]];
    ...
    epsilon ~ normal(0, 1);

and the problem disappeared.

[Andrew Gelman]

The Matt trick!