JAGS fits centered funnel just fine

[Bob Carpenter]

I ran the funnel in JAGS to follow along Daniel's section
in our course.

model {
for (n in 1:9) { x[n] ~ dnorm(0, 1 / exp(y / 2)^2) }
y ~ dnorm(0, 1/3^2)
}

and it appears to fit well out of the box (this is using R2jags):

> jags(data=list(), inits=list(list(), list(), list(), list()),
model.file="funnel.jags", parameters.to.save=c("x", "y"),
n.chains=4, DIC=FALSE)

Inference for Bugs model at "funnel.jags", fit using jags,
4 chains, each with 2000 iterations (first 1000 discarded)
n.sims = 4000 iterations saved
mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat n.eff
x[1] 0.011 7.988 -10.298 -0.530 0.006 0.592 9.085 1.011 2100
x[2] -0.003 8.950 -9.615 -0.576 -0.003 0.570 10.751 1.060 2200
x[3] -0.164 8.098 -9.115 -0.603 -0.009 0.528 9.498 1.011 4000
x[4] -0.205 8.004 -9.270 -0.648 -0.017 0.489 9.132 1.002 1400
x[5] -0.180 9.802 -9.485 -0.557 -0.006 0.554 8.847 1.029 4000
x[6] -0.197 8.416 -9.920 -0.627 -0.013 0.561 10.061 1.004 4000
x[7] -0.093 9.944 -9.484 -0.571 -0.003 0.556 10.595 1.031 4000
x[8] -0.020 9.464 -10.909 -0.649 0.001 0.528 9.320 1.008 2800
x[9] -0.026 14.375 -9.932 -0.557 -0.012 0.571 10.596 1.164 4000
y -0.007 3.026 -5.751 -2.067 -0.016 2.052 5.866 1.001 4000

And here's the fit from Stan with the non-centered parameterization:

parameters {
real y_raw;
vector[9] x_raw;
}
transformed parameters {
real y;
vector[9] x;
y <- 3 * y_raw;
x <- exp(y / 2) * x_raw;
}
model {
y_raw ~ normal(0, 1);
x_raw ~ normal(0, 1);
}

mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
y 0.03 0.05 3.01 -5.79 -1.98 -0.01 2.06 5.93 4000 1
x[1] -0.22 0.15 7.89 -10.04 -0.57 -0.01 0.54 9.51 2776 1
x[2] 0.08 0.22 9.15 -9.11 -0.55 0.00 0.61 11.44 1765 1
x[3] -0.07 0.14 8.16 -10.29 -0.62 0.00 0.57 10.01 3217 1
x[4] -0.44 0.16 8.98 -11.45 -0.61 -0.01 0.55 8.70 3037 1
x[5] 0.02 0.16 9.35 -10.16 -0.67 0.00 0.56 8.77 3251 1
x[6] 0.09 0.14 7.70 -10.83 -0.57 0.00 0.58 10.35 2943 1
x[7] 0.15 0.16 8.04 -9.14 -0.58 0.00 0.60 10.56 2435 1
x[8] 0.03 0.21 10.53 -9.33 -0.55 0.00 0.58 9.87 2520 1
x[9] 0.01 0.14 7.66 -8.48 -0.62 0.00 0.59 9.42 3014 1

So it looks like Gibbs does a reasonable job with the funnel out of
the box. Am I missing something here? I thought Michael said it
wouldn't work.

- Bob

[Andrew Gelman]

Bob,
I think it works ok here because your prior on the hyper-variance excludes zero. At least in the 8-schools example, the problem all comes in the neighborhood where the hyper-variance is close to 0.
A

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

[Bob Carpenter]

This is just the centered version of Neal's funnel example.
It may exclude 0 because it's on the log scale, but JAGS gets
down to y = -6 and exp(-6) = .0025.

One of the x still has R-hat (as calculated by Coda) of 1.15,
so it's not perfect.

Does 8 schools (which one, by the way?) not work in JAGS/BUGS?
I found JAGS refreshingly easy to install on the Mac (there's
a package installer for JAGS and then rjags just finds it on
install, then R2jags is easy to install on top of that --- took
about 2 minutes total). But I forgot just how clunky those
R2WinBUGS and R2jags interfaces were.

- Bob

[Andrew Gelman]

I don't mean that it excludes 0 because it's on the log scale, I mean that it excludes 0 because it has a lognormal prior. At least, that's how it looks. I'm getting a bit confused because of the non-standard parameterization.

Gibbs sampler works for 8 schools with Matt trick, does not work so well without Matt trick; we discuss in appendix C of BDA, I think.

[Bob Carpenter]

That's what I meant --- normal prior on the log of
the scale parameter. And it may exclude getting
an absolute zero value, but you get pretty close.

- Bob

[Michael Betancourt]

Can you plot the y distributions against each other to compare the tails?

[Andrew Gelman]

No, the lognormal has what is in practice a hard constraint near zero. It really is different than a uniform or gamma prior.

[Bob Carpenter]

By "hard constraint", I took you to mean that

lognormal(y | mu, tau) -> 0 as y -> 0

whereas the uniform doesn't have that property. Is
there more to it than that, or was the "no" about my
saying that y = 0.0025 was close to zero? It's not
nearly as close as you can get on a computer in double
precision --- that's 1e-300 or so.

- Bob

[Bob Carpenter]

Not sure what kind of plot you want so that the
tails are visible (let me know and I can generate),
but here's the summary for 100K iterations of both algorithms:

Stan non-centered:

mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat

y 0.01 0.01 3.02 -5.91 -2.02 0.00 2.06 5.94 199067 1

> summary(extract(fit)$y)
Min. 1st Qu. Median Mean 3rd Qu. Max.
-12.620000 -2.017000 0.004688 0.012290 2.058000 13.010000


JAGS centered:

mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat n.eff

y 0.006 2.995 -5.871 -2.008 0.009 2.024 5.858 1.001 130000

> summary(y)
V1
Min. :-13.521772
1st Qu.: -2.007790
Median : 0.008951
Mean : 0.005535
3rd Qu.: 2.023694
Max. : 12.901053

- Bob

[Andrew Gelman]

No, it's not merely that the density goes to 0, it's that it goes to 0 fast, so that in effect there's an epsilon>0 for which Pr (y < epsilon) is approx 0.

[Michael Betancourt]

Literally just histograms of both samples plotted on top of each other.
As seen in Figure 3 of http://arxiv.org/pdf/1312.0906.pdf (hard numbers
in Table I) the Metropolis-within-Gibbs did better than Random Walk
but still had some trouble getting down into the neck (and was much
slower than HMC). I’m guessing that you’d see a lighter tail in the
Gibbs samples, but perhaps not much.

[Krzysztof Sakrejda]

Bob, did you check which sampler jags is using? For some specific combos of nodes JAGS will use specialized samplers. They info is available in the interface.... somewhere...

[Bob Carpenter]

Nope --- if someone knows how to do that, I can check.

- Bob

> On Jun 3, 2016, at 8:12 AM, Krzysztof Sakrejda <krzysztof...@gmail.com> wrote:
>
> Bob, did you check which sampler jags is using? For some specific combos of nodes JAGS will use specialized samplers. They info is available in the interface.... somewhere...
>

[Avraham Adler]

On Friday, June 3, 2016 at 4:18:45 PM UTC-4, Bob Carpenter wrote:
> Nope --- if someone knows how to do that, I can check.
>
> - Bob
>
> > On Jun 3, 2016, at 8:12 AM, Krzysztof Sakrejda <krzysztof...@gmail.com> wrote:
> >
> > Bob, did you check which sampler jags is using? For some specific combos of nodes JAGS will use specialized samplers. They info is available in the interface.... somewhere...

Is the attached what you meant?

Avi

[Krzysztof Sakrejda]

list.factories("sampler") after fitting the model gets you the list for what's included. I thought there was a function to check what got used in a specific node but now I'm pretty sure Martyn is the only source of info on that...

[Avraham Adler]

On Friday, June 3, 2016 at 5:16:03 PM UTC-4, Krzysztof Sakrejda wrote:
> list.factories("sampler") after fitting the model gets you the list for what's included. I thought there was a function to check what got used in a specific node but now I'm pretty sure Martyn is the only source of info on that...

If you call JAGS from the command line, I believe you can. See section 3.1.13 of user manual (excerpt attached).

Avi

[Krzysztof Sakrejda]

Yep, that's the ticket. Thanks for coming up with that!

[Bob Carpenter]

I don't know how to run JAGS from the command line
and I don't have the patience to learn how to do it for
this example.

- Bob

> On Jun 3, 2016, at 5:40 PM, Krzysztof Sakrejda <krzysztof...@gmail.com> wrote:
>
> Yep, that's the ticket. Thanks for coming up with that!
>

[Bob Carpenter]

Was that the right paper/ I don't see any histograms in it.

But I'm giving up on this as I don't want to go down the
rabbit hole of figuring out what JAGS is actually doing given
that it seems hard to figure out.

- Bob

[Bob Carpenter]

The mins, maxes, and 2.5% and 97.5% quantiles were the same.

- Bob

> On Jun 2, 2016, at 10:37 PM, Michael Betancourt <betan...@gmail.com> wrote:
>

[Avraham Adler]

I don't know if this makes a difference, but the NEWS file for JAGS shows this change for 4.0:

* glm module
- Refactorization of samplers. The new Outcome class presents an
abstract normal approximation to the GLM engine.

Avi

[Bob Carpenter]

Michael was right. JAGS falls over on the centered parameterization
when the dimensionality increases. With n=100, it doesn't converge
even after 1M draws in four chains, even without diffuse inits.

- Bob