good news

[Andrew Gelman]

I got the 8 schools working (pretty much, if not perfect) by randomizing epsilon and L at each step. Not to stop periodicity but to allow HMC to get stuck from places where big steps or long chains will invariably get rejected. Currently I'm setting epsilon and L independent:
epsilon <- runif (1, 0, 2*epsilon.0)
L <- ceiling (2*L.0*runif(1))
But now I'm thinking it could be logical to have them be positively correlated.

[Bob Carpenter]

On 7/21/12 3:44 AM, Andrew Gelman wrote:> I got the 8 schools working (pretty much, if not perfect) by randomizing

epsilon and L at each step.
> Not to stop periodicity but to allow HMC to get stuck from places where big steps or long
> chains will invariably get rejected.

What I don't understand is that there's almost no
rejection.

Currently I'm setting epsilon and L independent:
> epsilon <- runif (1, 0, 2*epsilon.0)
> L <- ceiling (2*L.0*runif(1))
> But now I'm thinking it could be logical to have them be positively correlated.

What's epsilon.0?

Hmm. The BDA report from R2WinBUGS (p. 594) has:

sigma.theta = 6.5, sd = 5.7, n.eff = 23

With MCMC error = sd / sqrt(n.eff) = 1, that's
a pretty fat target. Do we have a better idea of what it should
really be? I'm not sure I trust the n.eff enough
to trust the MCMC se.

Here's Stan's vanilla HMC fit with epsilon = 0.05, L = 128, run
for 2K iterations with random inits in [-2,2]:

fit using Stan 1.0.0,
3 chains, each with 2000 iterations (first 1000 discarded)
n.sims = 3000 iterations saved
mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff
lp__ -15.8 7.9 -26.9 -21.5 -17.5 -11.9 4.8 1 350
theta[1] 11.2 8.7 -2.6 5.6 9.9 15.3 32.6 1 3000
theta[2] 7.6 6.1 -4.6 3.9 7.7 11.4 19.5 1 390
theta[3] 5.8 7.3 -10.9 1.9 6.4 10.3 19.3 1 260
theta[4] 7.6 6.2 -4.9 3.7 7.5 11.5 19.9 1 1100
theta[5] 5.1 6.1 -8.6 1.6 5.6 9.2 16.3 1 780
theta[6] 6.2 6.3 -7.9 2.5 6.3 10.2 18.1 1 610
theta[7] 10.5 6.9 -1.5 5.7 9.9 14.2 26.3 1 460
theta[8] 8.6 7.4 -5.7 4.1 8.3 12.8 25.0 1 370
mu_theta 7.8 4.9 -2.0 4.6 7.8 10.9 17.7 1 330
sigma_theta 6.3 5.4 0.2 2.2 5.1 8.9 19.4 1 350

NUTS is similar, but seems to systematically give higher estimates
for sigma_theta and a wider range of parameters (but roughly the
same mu_theta). For example, here's a run of 50K iterations
with default NUTS settings:

3 chains, each with 50000 iterations (first 25000 discarded), n.thin = 25
n.sims = 3000 iterations saved
mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff
lp__ -19.4 4.4 -27.8 -22.5 -19.3 -16.2 -11.5 1 840
treedepth__ 2.4 0.9 1.0 2.0 3.0 3.0 4.0 1 2300
stepsize__ 4.0 0.1 3.8 3.8 4.0 4.1 4.1 1 1
theta[1] 12.2 8.8 -3.0 6.4 11.1 17.1 32.1 1 3000
theta[2] 8.1 6.7 -4.9 3.7 8.0 12.2 21.2 1 3000
theta[3] 5.8 8.2 -11.8 1.1 6.0 10.8 21.3 1 3000
theta[4] 7.5 7.1 -6.7 2.9 7.5 12.0 21.9 1 660
theta[5] 4.4 6.7 -9.5 0.2 5.0 8.8 16.7 1 1000
theta[6] 5.6 7.1 -9.4 1.5 5.8 10.2 19.5 1 3000
theta[7] 11.3 7.1 -1.4 6.4 10.7 15.7 26.4 1 2500
theta[8] 8.6 8.6 -7.5 3.5 8.3 13.7 26.9 1 3000
mu_theta 8.0 5.4 -2.6 4.6 7.8 11.3 19.0 1 2700
sigma_theta 7.9 5.3 2.3 4.1 6.5 10.3 20.9 1 1300

- Bob

[Andrew Gelman]

See below.

On Jul 21, 2012, at 11:24 PM, Bob Carpenter wrote:

> On 7/21/12 3:44 AM, Andrew Gelman wrote:> I got the 8 schools working (pretty much, if not perfect) by randomizing epsilon and L at each step.
> > Not to stop periodicity but to allow HMC to get stuck from places where big steps or long
> > chains will invariably get rejected.
>
> What I don't understand is that there's almost no
> rejection.

What happens is that in the constrained parameter version (sigma unlogged, no bouncing), the chains go off to negative sigma -> rejection.
When I tried unlogged sigma _with_ bouncing, the chains, if they ever got near zero, would get ricocheted off to huge values of sigma -> rejection. The problem arose because of huge gradients on sigma.
When I tried log sigma, the chains would again sometimes go to extreme values.

The only way to avoid rejection was with tiny step sizes, and then there was no mixing.

[Bob Carpenter]

On 7/22/12 1:57 PM, Andrew Gelman wrote:
> See below.
>
> On Jul 21, 2012, at 11:24 PM, Bob Carpenter wrote:
>
>> On 7/21/12 3:44 AM, Andrew Gelman wrote:> I got the 8 schools working (pretty much, if not perfect) by randomizing epsilon and L at each step.
>>> Not to stop periodicity but to allow HMC to get stuck from places where big steps or long
>>> chains will invariably get rejected.
>>
>> What I don't understand is that there's almost no
>> rejection.
>
>
> What happens is that in the constrained parameter version (sigma unlogged, no bouncing), the chains go off to negative sigma -> rejection.

Right -- Stan would probably do that, too, if you don't
declare the variable properly. Another thing we might
want to try to catch and diagnose automatically.

> When I tried unlogged sigma _with_ bouncing, the chains, if they ever got near zero, would get ricocheted off to huge values of sigma -> rejection.

> The problem arose because of huge gradients on sigma.

So you get a huge gradient toward 0 (?) and then a big
bounce? I'd have thought that as sigma approaches 0,
the total gradient would push away from 0 unless the
coefficients were all very tightly clustered.

Maybe Michael can explain how he gets it to work.

> When I tried log sigma, the chains would again sometimes go to extreme values.

I'll have to check this out -- I wonder if that's what's distorting
the the quantiles.

> The only way to avoid rejection was with tiny step sizes, and then there was no mixing.

I was able to run small step sizes and long chains and get
mixing that was good enough to satisfy R-hat on multiple
chains from diffuse starting points.

- Bob

[Michael Betancourt]

On Jul 22, 2012, at 2:50 PM, Bob Carpenter wrote:

> So you get a huge gradient toward 0 (?) and then a big
> bounce? I'd have thought that as sigma approaches 0,
> the total gradient would push away from 0 unless the
> coefficients were all very tightly clustered.
>
> Maybe Michael can explain how he gets it to work.

Yeah, it sounds like there's not enough data to temper the prior and you're essentially trying to sample from 1/\sigma.

What isn't clear to me is why you're getting such extreme ricochets. As you approach zero from positive values the gradients increase, which increases the momenta. When you finally bounce, then, the reflected momenta will be large. But as you continue back through positive values the gradients should temper that momentum -- modulo the small epsilon you move past the origin the evolution should be symmetric with p -> -p. So unless you started at huge values of sigma I'm not sure how you're able to bounce to them.

Have you checked that the magnitude of the momenta along sigma is equal before and after the bounce?

[Andrew Gelman]

See below.

On Jul 22, 2012, at 8:50 PM, Bob Carpenter wrote:

>> What happens is that in the constrained parameter version (sigma unlogged, no bouncing), the chains go off to negative sigma -> rejection.
>
> Right -- Stan would probably do that, too, if you don't
> declare the variable properly. Another thing we might
> want to try to catch and diagnose automatically.

Yes, but Stan with Nuts should still work sorta ok in such situations.

>
>> When I tried unlogged sigma _with_ bouncing, the chains, if they ever got near zero, would get ricocheted off to huge values of sigma -> rejection.
>
> > The problem arose because of huge gradients on sigma.
>
> So you get a huge gradient toward 0 (?) and then a big
> bounce? I'd have thought that as sigma approaches 0,
> the total gradient would push away from 0 unless the
> coefficients were all very tightly clustered.

What happens is that when sigma is near zero, d(log p)/d(sigma) gets huge. Without bouncing, this sends HMC off to big negative values and rejection. With bouncing, it sends HMC off to big positive values. The problem is the persistence of the momentum variable. It goes like this: sigma is 0.6, say, and then the gradient is huge, then the momentum is huge, then it gets bounced off to sigma = 50 or so, now the gradient is negative but small, so it takes a long long long time for the momentum to get back negative. In the meantime, sigma becomes more and more extreme.

>
>> The only way to avoid rejection was with tiny step sizes, and then there was no mixing.
>
> I was able to run small step sizes and long chains and get
> mixing that was good enough to satisfy R-hat on multiple
> chains from diffuse starting points.
>

Yes, I got that too (with log sigma) but then I was getting 99-100% accept rates, which didn't make me so happy. I preferred the results with random epsilon and L at each step.

[Andrew Gelman]

M
I don't understand the question. When I bounce, I simply set sigma <- - sigma and momentum(sigma) <- - momentum (sigma). Here's the R code below. Do I have a mistake? Here, "th" is the vector of 10 parameters (sigma is theta[10]) and "phi" is momentum.
A

hmc.iteration <- function (th, y, s.y, L, epsilon, M){
d <- length (th)
phi <- rnorm (d, 0, 1)
th.old <- th
log.p.old <- log.p.th (th,y,s.y) - 0.5*sum(phi^2)
phi <- phi + 0.5*epsilon*gradient.th (th, y, s.y)
for (l in 1:L){
th <- th + epsilon*phi
if (th[10] < 0) {
th[10] <- - th[10]
phi[10] <- - phi[10]
}
phi <- phi + (if (l==L) 0.5 else 1)*epsilon*gradient.th (th, y, s.y)
}
phi <- -phi
log.p.star <- log.p.th (th,y,s.y) - 0.5*sum(phi^2)
r <- exp (log.p.star - log.p.old)
if (is.nan(r)) r <- 0
p.jump <- min(r,1)
th.new <- if (runif(1) < p.jump) th else th.old
return (list (th=th.new, p.jump=p.jump))

[Michael Betancourt]

On Jul 22, 2012, at 3:15 PM, Andrew Gelman wrote:

> What happens is that when sigma is near zero, d(log p)/d(sigma) gets huge. Without bouncing, this sends HMC off to big negative values and rejection. With bouncing, it sends HMC off to big positive values. The problem is the persistence of the momentum variable. It goes like this: sigma is 0.6, say, and then the gradient is huge, then the momentum is huge, then it gets bounced off to sigma = 50 or so, now the gradient is negative but small, so it takes a long long long time for the momentum to get back negative. In the meantime, sigma becomes more and more extreme.

Oh, oh, oh. Are you bouncing right before you violate the constraint or afterwards? You want to bounce as soon as you're in violation of the constraint (the first time the spatial update passes q = 0) but not before that. As I mentioned in my response to Bob, you should be seeing roughly symmetry bounces and not huge jumps into regions in which you have not yet been.

See Algo 1 (page 13) here: http://arxiv.org/pdf/1005.0157.pdf. In this case C(x) = x and (grad C / | gradC |) = 1.

[Michael Betancourt]

Ah! Sigma stays the same, only the momentum should change in the bounce (which is really a replacement for the usual momentum update).

If you remove the line "th[10] <- - th[10]" everything should behave much better.

[Andrew Gelman]

Really? But then sigma is negative for that one step. That's ok?
A

[Michael Betancourt]

Yup. If you're unlucky and you end on that step then it's just a rejection. It's one of the reasons I like using accept/reject windows -- if the window size is greater than 1 then you're guaranteed to never reject to due a constraint violation (provided you don't run into pathologies, etc).

[Andrew Gelman]

M
Maybe this doesn't matter because I could just forget the damn bouncing just for simplicity, but let me make sure I'm ok…. Here is the key part of the algorithm as I've programmed it:

phi <- phi + 0.5*epsilon*gradient.th (th, y, s.y)
for (l in 1:L){
th <- th + epsilon*phi
if (th[10] < 0) {
th[10] <- - th[10]
phi[10] <- - phi[10]
}
phi <- phi + (if (l==L) 0.5 else 1)*epsilon*gradient.th (th, y, s.y)
}
phi <- -phi

Are you saying that if I just eliminate the line, "th[10] <- - th[10]", then I'm ok? So it's alright that I first move theta to the illegal zone, then I flip the sign of phi, then I update phi? The difficulty here is that my update is based on the gradient of the updated (negative) theta. Is this right?

A

[Michael Betancourt]

The bounce is a replacement for, and not an addition to, the regular momentum update. You just have to change your code to an if/else

phi <- phi + 0.5*epsilon*gradient.th (th, y, s.y)
for (l in 1:L){
th <- th + epsilon*phi
if (th[10] < 0) {

phi[10] <- - phi[10]
}

else {

phi <- phi + (if (l==L) 0.5 else 1)*epsilon*gradient.th (th, y, s.y)
}
phi <- -phi

[Andrew Gelman]

And all this time I thought I was faithfully following Radford's algorithm…

[Andrew Gelman]

P.S. Still slow to converge. Still gives chains with zero acceptances if I don't randomize epsilon and L. But at least I'm (maybe) doing the right thing. I'm happy we have nuts.

On Jul 22, 2012, at 10:13 PM, Michael Betancourt wrote:

[Michael Betancourt]

Have you tried the correlated randomization of epsilon and L so that epsilon * L ~ const? I'm wondering if that would give you 50% of the performance while being much simpler for pedagogical purposes.

[Andrew Gelman]

M
I'll try, but I actually think that a positive corr of epsilon and L will work better!
A