Help with divergent transitions + noncentered paramterization

68 views
Skip to first unread message

Jökull Snæbjarnarson

unread,
May 30, 2017, 11:10:12 AM5/30/17
to Stan users mailing list
Can anyone of you help me out diagnosing the problem with this simple example model and maybe give some on advice on how to fix it ?
It gives me a bunch of  divergent transitions. Funnel like correlations and a highly curved posterior are apparent but I don't really get where those correlations come from.
Maybe I can implement the noncentered paramterization in a different way?

Thanks in advance,
Jökull

nonCentered <- "
data {
  int<lower=1> Ntotal;
  vector[Ntotal] pred;
  vector[Ntotal] pred2;
  vector[Ntotal] y;
}
parameters {
  real<lower=0.0> sigma;
  real<lower=0.0> sigma_w;
  real<lower=0.0> sigma_w2;
  real mu;
  real w_raw;
  real w_raw2;
}
transformed parameters {
  real w;
  real w2;
  vector[Ntotal] y_mu;
  for (n in 1:Ntotal) y_mu[n] = 0.0;
  w = sigma_w * w_raw + mu;
  w2 = sigma_w2 * w_raw2 + mu;
  y_mu = w * pred + w2 * pred2;
}
model {
  sigma_w ~ normal(1,1);
  sigma_w2 ~ normal(1,1);
  sigma ~ normal(1,1);
  mu ~ normal(1,0.2);
  w_raw ~ normal(0, 1);
  w_raw2 ~ normal(0, 1);
  y ~ normal(y_mu,sigma);
}
"

fit <-
  stan(
    model_code = nonCentered,
    data = list(pred = c(0,2,0,0,2,0,0,2,0),
                pred2 = c(1,0,1,1,0,1,1,0,1),
                Ntotal =9, y = c(3,2,3,3,2,3,3,2,3)),
    model_name =  'example',
    iter = 1200,
    chains = 2
  )




Andrew Gelman

unread,
May 30, 2017, 4:12:59 PM5/30/17
to stan-...@googlegroups.com
Hi, just a quick comment, you can just say 0, you don't need to say 0.0.  Also, this line:
 for (n in 1:Ntotal) y_mu[n] = 0.0;
can become simply this:
  y_mu = 0;
This won't solve your problem but it should make the model easier to read.
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.
To post to this group, send email to stan-...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Bob Carpenter

unread,
May 30, 2017, 6:47:07 PM5/30/17
to stan-...@googlegroups.com

> On May 30, 2017, at 4:12 PM, Andrew Gelman <gel...@stat.columbia.edu> wrote:
>
> Hi, just a quick comment, you can just say 0, you don't need to say 0.0. Also, this line:
> for (n in 1:Ntotal) y_mu[n] = 0.0;
> can become simply this:
> y_mu = 0;
> This won't solve your problem but it should make the model easier to read.

Not quite. You'd have to set

y_mu = rep_vector(0, Ntotal);

As to the bigger problem, how much data do you have? The centered
parameterization is better when there's more data.

- Bob

Jökull Snæbjarnarson

unread,
May 31, 2017, 5:11:18 AM5/31/17
to Stan users mailing list
Hey thanks for the quick replies and comments! 

I posted the data; it's 9 data point. 

Bob Carpenter

unread,
May 31, 2017, 4:22:06 PM5/31/17
to stan-...@googlegroups.com
Unless those 9 data points have very low variance, you're
probably going to be better off with the non-centered
parameterization.

- Bob

stijn masschelein

unread,
Jun 2, 2017, 5:26:06 AM6/2/17
to Stan users mailing list
I think a lot of your troubles come from the fact that your predictors have a perfect linear relationship with y. y = 3 - .5 * pred or y = 2 + pred2. If you would have an intercept in your model and you would be modeling y ~ normal(a0 + a1 pred + a2 pred2, sigma) than the nuts sampler would have to switch between the equally probably modes (a0 = 3, a1 = -.5, a2 = 0) and (a0 = 2, a1 = 0, a2 = 1). This is likely to lead to divergences and weird curvature. 

If I read your model correctly, there is no intercept a0 in your model but I still think that the linear relations in your data (see also pred = 2 - 2*pred2) are causing havoc.

Stijn

Bob Carpenter

unread,
Jun 4, 2017, 8:33:07 PM6/4/17
to stan-...@googlegroups.com
Theres a discussion of just this issue in the problematic posteriors
chapter of the manual. It can be partially tamed with priors, but
it's usually much better to just identify the model if you know there's
a problem.

- Bob
Reply all
Reply to author
Forward
0 new messages