problem with bivariate normal!!!

[Andrew Gelman]

Hi, all. This is the problem that I was having the other day, I stripped down the model and stripped it down and stripped it down until I came up with the following two examples. I can't figure out what is going on! It's the simplest Gaussian model, it works fine in one dimension but not in two dimensions. I don't know if it's a bug in Stan, or in Nuts, or something stupid I'm doing, or what??? I'm doing 20,000 iterations so I don't think it is a convergence issue.
See you
A

---

In file gaussian1.stan:

parameters {
real b1;
}
model {
b1 ~ normal (0, 1);
}

---

In file gaussian2.stan:

parameters {
real b1;
real b2;
}
model {
b1 ~ normal (0, 1);
b2 ~ normal (0, 1);
}

---

R script:

library ("rstan")
gaussian1 <- stan (file="gaussian1.stan", iter=20000, chains=4)
sims <- extract(gaussian1)$b1
print (mean (sims))
print (var (sims))
gaussian2 <- stan (file="gaussian2.stan", iter=20000, chains=4)
sims <- extract(gaussian2)$b1
print (mean (sims))
print (var (sims))

---

Output in R console [cutting out the deals on the Stan sampling, which showed no problems]:

> sims <- extract(gaussian1)$b1
> print (mean (sims))
[1] 0.006627426
> print (var (sims))
[1] 0.9961059

> sims <- extract(gaussian2)$b1
> print (mean (sims))
[1] -0.004989011
> print (var (sims))
[1] 0.9606021

---

This happens consistently. It started with a more complicated regression model, then I stripped it down. The bivariate normal consistently converges in Stan to have too low a variance. But the univariate is OK. This is what made me worry that there's some sort of Nuts problem. But I have no idea. Maybe it is some implicit constraint happening within Stan? Or something else? I also noticed the problem when b was a vector; I expressed it above as two scalars (b1 and b2) just to make sure that vectorization wasn't the problem.

[Bob Carpenter]

I verified that this is a bug in 2.0.1. I also verified
that our current development branch (what will soon be 2.1.0)
does NOT have this problem.

We need to write to our user's list ASAP about what's going on.
Does someone other than me want to do this?

Clearly N_eff and R_hat aren't enough for diagnosing problems
like these.

It's very worrisome that we haven't been more careful testing
posterior variance estimates. We need to add this as a test
case before 2.1.0 hits the streets.

Here are the answers I get with Andrew's model (b1 and b2
both unit normal).

CmdStan 2.0.1, default adaptation:

Mean MCSE StdDev 5% 50% 95% N_Eff N_Eff/s R_hat
b1 1.775e-03 3.251e-03 9.515e-01 -1.570e+00 1.433e-03 1.567e+00 8.565e+04 1.161e+04 1.000e+00
b2 -3.264e-03 3.202e-03 9.487e-01 -1.571e+00 -4.142e-03 1.562e+00 8.778e+04 1.189e+04 1.000e+00


CmdStan 2.0.1, adapt delta=0.99:

b1 1.101e-04 3.721e-03 9.226e-01 -1.517e+00 3.211e-03 1.522e+00 6.148e+04 4.121e+03 1.000e+00
b2 -2.978e-03 3.703e-03 9.185e-01 -1.515e+00 -2.584e-03 1.516e+00 6.153e+04 4.124e+03 1.000e+00


CmdStan development branch (i.e., 2.1.0), default adaptation:

Mean MCSE StdDev 5% 50% 95% N_Eff N_Eff/s R_hat
b1 -2.045e-03 4.000e-03 1.001e+00 -1.645e+00 -3.660e-03 1.638e+00 6.267e+04 3.546e+04 1.000e+00
b2 2.257e-03 4.014e-03 9.966e-01 -1.637e+00 5.150e-03 1.640e+00 6.166e+04 3.489e+04 1.000e+00


- Bob

[Michael Betancourt]

This doesn't make any sense -- the only real difference was a tweak in how the average
acceptance probability (and hence step size) was being calculated. But the fact that
increasing the average acceptance probability makes the problem worse goes in the
wrong direction.

The median is correct and the distribution is appropriately symmetric. The only thought
I have is that the integrator might have been diverging past 1.5 sigma, but in that case
increasing the target acceptance probability (and decreasing the step size) would help.
And increasing the dimensionality by one wouldn't have any big effect.

I'll look into it.

> --
> 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/groups/opt_out.

[Andrew Gelman]

It's scary, that's for sure! I thought Nuts was supposed to work for any stepsize, at least for something as well behaved as the Gaussian. Also mysterious that it works in 1-D but not in 2-D.

I'd wait on announcing to users-list until we figure out where the problem is.

[Bob Carpenter]

I thought Michael had said before that he'd found a bug
somewhere in the 2.0.x code. I can't recall the discussion,
though.

If this helps, here's the rest of the bin/print output, rerun with 2.0.1
having target acceptance of 0.8:

2.0.1, adapt delta=0.8

Mean MCSE StdDev 5% 50% 95% N_Eff N_Eff/s R_hat

lp__ -9.029e-01 4.647e-03 9.485e-01 -2.804e+00 -5.939e-01 -4.457e-02 4.166e+04 5.065e+03 1.000e+00
accept_stat__ 8.989e-01 3.957e-04 1.251e-01 6.393e-01 9.458e-01 1.000e+00 1.000e+05 1.216e+04 1.000e+00
stepsize__ 9.896e-01 5.389e-13 3.810e-13 9.896e-01 9.896e-01 9.896e-01 5.000e-01 6.079e-02 1.000e+00
treedepth__ 7.100e-01 1.461e-03 4.553e-01 0.000e+00 1.000e+00 1.000e+00 9.707e+04 1.180e+04 1.000e+00
b1 5.071e-03 3.123e-03 9.496e-01 -1.562e+00 6.271e-03 1.575e+00 9.245e+04 1.124e+04 1.000e+00
b2 -9.406e-04 3.164e-03 9.508e-01 -1.574e+00 -1.834e-03 1.569e+00 9.032e+04 1.098e+04 1.000e+00


dev branch (2.1.0), default adapt params:

Mean MCSE StdDev 5% 50% 95% N_Eff N_Eff/s R_hat

lp__ -1.0e+00 5.4e-03 1.0e+00 -3.0e+00 -6.8e-01 -5.1e-02 3.5e+04 19660 1.0e+00
accept_stat__ 8.9e-01 4.1e-04 1.3e-01 6.2e-01 9.4e-01 1.0e+00 1.0e+05 56581 1.0e+00
stepsize__ 9.8e-01 4.4e-13 3.1e-13 9.8e-01 9.8e-01 9.8e-01 5.0e-01 0.28 1.0e+00
treedepth__ 5.4e-01 1.6e-03 5.0e-01 0.0e+00 1.0e+00 1.0e+00 9.8e+04 55556 1.0e+00
n_divergent__ 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 1.0e+05 56581 nan
b1 -2.0e-03 4.0e-03 1.0e+00 -1.6e+00 -3.7e-03 1.6e+00 6.3e+04 35458 1.0e+00
b2 2.3e-03 4.0e-03 1.0e+00 -1.6e+00 5.1e-03 1.6e+00 6.2e+04 34885 1.0e+00


This is just a single chain of 100K warmup, 100K samples. The acept
stats and step sizes seem very close, though the mean treedepth is much lower
in 2.1.0.

- Bob

[Michael Betancourt]

This is the bug Bob was mentions: c59c02d8f0eab4491912b0e3a6f3e2ef3067a2ad, but it doesn't
seem to be that easy. At that commit (Oct 20th) everything looks fine,

b1 5.2e-03 1.3e-02 9.9e-01 -1.7 4.4e-03 1.7e+00 5485 38334 1.00
b2 -1.7e-02 1.3e-02 9.8e-01 -1.6 -3.1e-02 1.6e+00 5713 39931 1.0

But down the line at, say 2ae82f857d9f4919a8a2056aecd5a9101c5c260d (Nov 1) the problem arises

b1 -1.1e-02 1.1e-02 9.7e-01 -1.6 -7.5e-03 1.6e+00 8298 50451 1.0e+00
b2 1.1e-02 1.0e-02 9.6e-01 -1.6 3.4e-03 1.6e+00 8665 52680 1.0e+00

Perhaps I'm getting my times messed by the merging process?

[Michael Betancourt]

Yeah, I may have screwed up the timing -- those were the commit dates not the merge dates. My guess
is that it's c59c02d8f0eab4491912b0e3a6f3e2ef3067a2ad, although I'm not sure how to do a proper
before and after.

[Bob Carpenter]

It would be nice to know what caused the problem.
Sounds like this is a good candidate for git bisect.

While knowing what went wrong in the past is helpful,
I think it's more important to include a test so we can catch
it if we reintroduce the problem somehow.

I just pulled the latest dev branch (last report was from
last night's) and re-ran and it looks OK:

Mean MCSE StdDev 5% 50% 95% N_Eff N_Eff/s R_hat

lp__ -1.005e+00 5.898e-03 9.974e-01 -2.986e+00 -7.016e-01 -5.223e-02 2.860e+04 3.002e+03 1.000e+00
accept_stat__ 9.254e-01 2.855e-04 9.029e-02 7.384e-01 9.584e-01 1.000e+00 1.000e+05 1.050e+04 1.000e+00
stepsize__ 8.257e-01 3.333e-13 2.357e-13 8.257e-01 8.257e-01 8.257e-01 5.000e-01 5.250e-02 1.000e+00
treedepth__ 7.842e-01 1.998e-03 6.317e-01 0.000e+00 1.000e+00 2.000e+00 9.992e+04 1.049e+04 1.000e+00
n_divergent__ 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 1.000e+05 1.050e+04 nan
b1 7.519e-04 4.181e-03 1.012e+00 -1.666e+00 4.777e-03 1.660e+00 5.855e+04 6.147e+03 1.000e+00
b2 2.316e-03 4.077e-03 9.936e-01 -1.626e+00 2.559e-03 1.628e+00 5.938e+04 6.234e+03 1.000e+00

I'm not clear on how to report hashes for this because most
of them are local merges. My log (I'm in the middle of working
on the new manual) shows:

~/stan(feature/issue-324-manual-2.1)$ git log

commit 7f1359c0ed144e9a1dfa1c8d121534bc45d55fc4
Merge: e88f247 b796536
Author: Bob Carpenter <ca...@alias-i.com>
Date: Tue Dec 17 14:47:41 2013 -0500

Merge branch 'develop' of https://github.com/stan-dev/stan into feature/issue-324-manual-2.1

commit b79653664ee1e61041b7568c72f05f6a4acc6b33
Merge: d203b46 6683aab
Author: Daniel Lee <bea...@alum.mit.edu>
Date: Tue Dec 17 07:56:55 2013 -0800

Merge pull request #450 from stan-dev/feature/issue-449-fix-initialization

Fixes #450.


- Bob

[Andrew Gelman]

It's totally wack that it happened in 2D and not in 1.
But I am proud that I discovered a major bug!

[Michael Betancourt]

> It would be nice to know what caused the problem.

The issue in the noted commit was that the left-hand side of the parent tree wasn't being correctly passed
recursively to the sub trees. This caused premature termination but no obvious bias that I could see, although
I did see some effect on the funnel as noted in the commit comments.

> While knowing what went wrong in the past is helpful,
> I think it's more important to include a test so we can catch
> it if we reintroduce the problem somehow.

I think we are all aware of the need to test sampler validity. Once Dan reorganizes the tests (making it much
easier to run only a few at once) I'll add a rudimentary sampler testing framework. Hopefully it will be
extensive enough to quickly add difficult models as we find them.

> I'm not clear on how to report hashes for this because most
> of them are local merges. My log (I'm in the middle of working
> on the new manual) shows:

The hash is correct for the changes, it just may not reflect the history of our flow. Remember that merging
concatenates the history of the two branches.