Model optimization leads to "stuck" chain

[Sebastian Weber]

Hi!

I am trying to right now to speed up my mixed effect logistic regression model by applying the location-scale separation trick as suggested in the FAQ. However, when applying the trick it seems to work if I only run 4 chains then the ESS goes up. BUT if I increase the chains to 8 the "fast" model produces some chains which constantly produces errors like

Warning (non-fatal rejection): Error in function stan::prob::normal_log(N4stan5agrad3varE): Random variable is nan:0, but must not be nan!

and the result is that at least one chain does not converge at all. I attached the toy model and plots of the traces. The error only occurs if there are a reasonable number of chains (for me 8 are sufficient) as this phenomena seems to depend on initial conditions. Is the bug now me coding a wrong model or is stan having an issue here?

Best,

Sebastian

For the record:

R version 2.15.1 Patched (2012-10-07 r60893)
Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit)

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/de_DE.UTF-8

attached base packages:
[1] splines stats graphics grDevices utils datasets methods
[8] base

other attached packages:
[1] foreach_1.4.0 rstan_1.0.2 RcppEigen_0.3.1 Matrix_1.0-9
[5] lattice_0.20-10 inline_0.3.9 Rcpp_0.9.14 plyr_1.7.1
[9] nlme_3.1-104 Hmisc_3.9-3 survival_2.36-14 hdf5_1.6.9
[13] ggplot2_0.9.2.1 reshape2_1.2.1

loaded via a namespace (and not attached):
[1] cluster_1.14.2 codetools_0.2-8 colorspace_1.1-1 compiler_2.15.1
[5] dichromat_1.2-4 digest_0.5.2 grid_2.15.1 gtable_0.1.1
[9] iterators_1.0.6 labeling_0.1 MASS_7.3-22 memoise_0.1
[13] munsell_0.4 proto_0.3-9.2 RColorBrewer_1.0-5 scales_0.2.2
[17] stringr_0.6.1 tools_2.15.1


--
Sebastian Weber
Group of Cell Communication and Control
Freiburg Institute for Advanced Studies - FRIAS
School of Life Sciences - LIFENET
Albert-Ludwigs-Universität Freiburg
Albertstr. 19
79104 Freiburg
T.: +49-761-203-97237
Fax:+49-761-203-97334

[Ben Goodrich]

Hi Sebastian,

On Saturday, October 13, 2012 4:48:24 AM UTC-4, Sebastian Weber wrote:

I am trying to right now to speed up my mixed effect logistic regression model by applying the location-scale separation trick as suggested in the FAQ. However, when applying the trick it seems to work if I only run 4 chains then the ESS goes up. BUT if I increase the chains to 8 the "fast" model produces some chains which constantly produces errors like

Warning (non-fatal rejection): Error in function stan::prob::normal_log(N4stan5agrad3varE): Random variable is nan:0, but must not be nan!

and the result is that at least one chain does not converge at all. I attached the toy model and plots of the traces. The error only occurs if there are a reasonable number of chains (for me 8 are sufficient) as this phenomena seems to depend on initial conditions. Is the bug now me coding a wrong model or is stan having an issue here?

I should have said (and will clarify) in the FAQ that the Matt trick essentially presumes proper priors on the higher-level parameters, which in your case are the expectation and standard deviation of the betas. When I do that, I seem to get good performance again. I'm not entirely sure why improper priors don't also undermine the performance of the first version of the model, but I haven't spent much time worrying about that. Here is my version of your model_fast with some additional commented tweaks to speed it up.

data {
   int<lower=0> Nmax;
   int<lower=0> S;
   int<lower=0> N;
   int<lower=0,upper=Nmax> y[N];
   real time[N];
   int<lower=0,upper=S> Cs[N];
}
parameters {
   real alpha;
   real mu_beta;
   real<lower=0> sigma;
   vector[S] e_beta;
}
transformed parameters {
   vector[S] beta;
   beta <- mu_beta + e_beta * sigma; // no loop necessary when beta and e_beta are vectors
}
model {
   real<lower=0,upper=1> p[N]; // local parameter, so not stored
   for(i in 1:N) p[i] <- inv_logit(alpha + beta[Cs[i]] * time[i]);

   mu_beta ~ cauchy(0,5); // proper prior; other choices possible
   sigma ~ cauchy(0,5);   // since sigma > 0, symmetry implies a folded-Cauchy; other choices possible

   e_beta ~ normal(0,1); // no loop necessary
   for(i in 1:N) y[i] ~ binomial(Nmax, p[i]);
}

Ben

[Sebastian Weber]

Hi Ben,

thank you very much. When I tried your code, it first did not compile for me since stan complained about the 0-1 limitation for p in the model block. Are limitations not allowed there? Removing the upper and lower statements solved the problem and the model seems to work fine. However, I tried to get rid of the p variable altogether and wrote

for(i in 1:N) y[i] ~ binomial(Nmax, inv_logit(alpha + beta[Cs[i]] * time[i]) );

in my model. This again produces the error that some chains do not converge at all and I see again the

Warning (non-fatal rejection): Error in function stan::prob::cauchy_log(N4stan5agrad3varE): Random variable is nan:0, but must not be nan!

statement for one of the chains. This looks to me like bug - either in stan, the compiler or wherever; if you agree, I can provide more info.

Some more questions in this context:

- Concerning the inv_logit transform: To me it is simply a transform of a given parameter, so why shouldn't I put it in the transformed parameter block? Should I always put inverse links into the model block? And if yes, which transformations should I then included in transformed parameters and which in the model?

- should I always limit my scale parameters by 0 and some upper bound, lets say 100 in the declarations or is a cauchy prior preferable?

Best,

Sebastian

BTW, the schools example on the stan R intro contains a scale parameter tau, which has no proper prior and there the same shift-and-scale trick is applied:

http://code.google.com/p/stan/wiki/RStanGettingStarted

Maybe, a proper prior there is appropriate as well then.

> --

[Ben Goodrich]

On Sunday, October 14, 2012 4:29:09 AM UTC-4, Sebastian Weber wrote:

thank you very much. When I tried your code, it first did not compile for me since stan complained about the 0-1 limitation for p in the model block. Are limitations not allowed there? Removing the upper and lower statements solved the problem and the model seems to work fine.

Sorry about that. Restrictions on the support of a local variable are not allowed (and are not checked) in the model {} block. But since the variable is on the (0,1) interval by construction, it is fine to remove the lower and upper bounds. I believe that transformations in the transformed parameters block are checked, so it is good to provide the bounds in that case.
 

However, I tried to get rid of the p variable altogether and wrote

for(i in 1:N) y[i] ~ binomial(Nmax, inv_logit(alpha + beta[Cs[i]] * time[i]) );

in my model. This again produces the error that some chains do not converge at all and I see again the

Warning (non-fatal rejection): Error in function stan::prob::cauchy_log(N4stan5agrad3varE): Random variable is nan:0, but must not be nan!

statement for one of the chains. This looks to me like bug - either in stan, the compiler or wherever; if you agree, I can provide more info.

When I do it with a large number of chains and a large number of iterations, sometimes it goes perfectly and sometimes one chain experiences a nan and never recovers. Many functions in Stan are not numerically stable, so things like this happen with some probability.
 

Some more questions in this context:

- Concerning the inv_logit transform: To me it is simply a transform of a given parameter, so why shouldn't I put it in the transformed parameter block? Should I always put inverse links into the model block? And if yes, which transformations should I then included in transformed parameters and which in the model?

The question is really: Do you want to analyze its posterior distribution? If yes, put it into the transformed parameters {} block. If no, put it in the model {} block. I assumed no since your R code did not reference those probabilities, but whichever choice you make "should" have no effect on the sampling, in theory. But in practice, maybe numeric instability issues are handled more gracefully when the probability is in the transformed parameters {} block because if it violates its bounds, then a non-fatal exception is thrown before the log-posterior is evaluated.
 

- should I always limit my scale parameters by 0 and some upper bound, lets say 100 in the declarations or is a cauchy prior preferable?

In principle, it is a question of what you believe about the data-generating process. Andrew's paper recommended either of those, and maybe bounding on both sides is more numerically stable. What Andrew found to be objectionable in theory and problematic in practice were (nearly) improper priors over the entire positive half of the real line.

Ben

[Bob Carpenter]

On 10/14/12 12:39 PM, Ben Goodrich wrote:
> On Sunday, October 14, 2012 4:29:09 AM UTC-4, Sebastian Weber wrote:
>
> thank you very much. When I tried your code, it first did not compile for me since stan complained about the 0-1
> limitation for p in the model block. Are limitations not allowed there? Removing the upper and lower statements
> solved the problem and the model seems to work fine.
>
>
> Sorry about that. Restrictions on the support of a local variable are not allowed (and are not checked) in the model {}
> block. But since the variable is on the (0,1) interval by construction, it is fine to remove the lower and upper bounds.
> I believe that transformations in the transformed parameters block are checked, so it is good to provide the bounds in
> that case.

Yes, anywhere you can declare constraints, the
constraints are enforced at the end of the block
in which the variable is declared.

> When I do it with a large number of chains and a large number of iterations, sometimes it goes perfectly and sometimes
> one chain experiences a nan and never recovers. Many functions in Stan are not numerically stable, so things like this
> happen with some probability.

There are also basic functions, like 0.0 / 0.0 that
produce NaN. In some cases, we can get around that.

If there are known places where you think we're doing
something at boundaries that we could make better
behaved than returning NaN or throwing an exception, let
us know and we can fix it.

> Some more questions in this context:
>
> - Concerning the inv_logit transform: To me it is simply a transform of a given parameter, so why shouldn't I put it
> in the transformed parameter block? Should I always put inverse links into the model block? And if yes, which
> transformations should I then included in transformed parameters and which in the model?
>
>
> The question is really: Do you want to analyze its posterior distribution? If yes, put it into the transformed
> parameters {} block. If no, put it in the model {} block. I assumed no since your R code did not reference those
> probabilities, but whichever choice you make "should" have no effect on the sampling, in theory.

We can work out the behavior based on the definitions
of the error checking in the manual.

Transformed parameters: constraints enforced at the end
of the model block

Parameters to probability functions: domain constraints
enforced when called

Given this behavior, if you have a variable that has to
be between 0 and 1, say because it's the parameter for a Bernoulli,
then declaring a constraint on the transformed parameter
catches the error before any of the model is executed. Defining
it as a local variable in the model block will lead to the error
being caught the first time it's used.

All else being equal, trapping errors early is a good thing.

> But in practice, maybe
> numeric instability issues are handled more gracefully when the probability is in the transformed parameters {} block
> because if it violates its bounds, then a non-fatal exception is thrown before the log-posterior is evaluated.

Please let us know if you run into reproducible numerical
stability issues. We could definitely use better limiting
behavior in many cases.

In some cases, it's just plain old IEEE floating-point
arithmetic, like if you divde 0.0 / 0.0 somewhere.

> - should I always limit my scale parameters by 0 and some upper bound, lets say 100 in the declarations or is a
> cauchy prior preferable?
>
>
> In principle, it is a question of what you believe about the data-generating process. Andrew's paper recommended either
> of those, and maybe bounding on both sides is more numerically stable. What Andrew found to be objectionable in theory
> and problematic in practice were (nearly) improper priors over the entire positive half of the real line.

You can try both and see what effect they have on the posterior.
They will also have different posterior geometries (that is, the
covariance structure across the parameter values), which can make
one or the other more efficient in terms of drawing effective
samples.

- Bob

[Bob Carpenter]

On 10/14/12 4:29 AM, Sebastian Weber wrote:
> Hi Ben,
>
> thank you very much. When I tried your code, it first did not compile for me since stan complained about the 0-1 limitation for p in the model block. Are limitations not allowed there? Removing the upper and lower statements solved the problem and the model seems to work fine. However, I tried to get rid of the p variable altogether and wrote
>
> for(i in 1:N) y[i] ~ binomial(Nmax, inv_logit(alpha + beta[Cs[i]] * time[i]) );
>
> in my model. This again produces the error that some chains do not converge at all and I see again the
>
> Warning (non-fatal rejection): Error in function stan::prob::cauchy_log(N4stan5agrad3varE): Random variable is nan:0, but must not be nan!
>
> statement for one of the chains. This looks to me like bug - either in stan, the compiler or wherever; if you agree, I can provide more info.

No, it's not a bug per se, but it is a problem with
floating point arithmetic.

We should have a binomial_logit like we have a bernoulli_logit,
which is better behaved numerically because we can do some
algebra behind the scenes as the log probabilities are
calculated.

We really need to make this error message clearer.
It's confusing everyone. Maybe something that says
it's not a bug.

The issue is not only the math library, but also the
Hamiltonian Monte Carlo simulations, which if they
aren't discretized at a fine enough grain can lead
to errors in the transforms. But the problem isn't
always there, it's then when those parameters are used.

> Some more questions in this context:
>
> - Concerning the inv_logit transform: To me it is simply a transform of a given parameter, so why shouldn't I put it in the transformed parameter block? Should I always put inverse links into the model block? And if yes, which transformations should I then included in transformed parameters and which in the model?

You could put it there. See my previous message.

In Stan, the decision should be made based on whether
you want the error checking in the transformed parameter
block and if you want printing. It's too bad they're
linked together, because they're really independent
functions.

You don't need to do the BUGS thing, like this:

for (n in 1:N) y_hat[n] <- alpha * x[n] + beta;
for (n in 1:N) y[n] ~ normal(y_hat[n],sigma);

If you declare y_hat[n] as a transformed parameter,
you get n outputs, which you may not care about. They'll
take space and time to write, but it might be negligible
compared to everything else going on depending on the model.

You can just write:

for (n in 1:N) y[n] ~ normal(alpha * x[n] + beta, sigma);

or you can write:

for (n in 1:N) {
real y_hat;
y_hat <- alpha * x[n] + beta;
y[n] ~ normal(y_hat, sigma);
}

These won't output the values.

- Bob

[Ben Goodrich]

On Sunday, October 14, 2012 4:48:29 PM UTC-4, Bob Carpenter wrote:

Please let us know if you run into reproducible numerical
stability issues.  We could definitely use better limiting
behavior in many cases.

If you source("stability.txt") in R, a handful of chains per 100 should exhibit these problems.

Ben