numerical instability

[Sebastian Weber]

Dear Stan users,

I am experiencing rather weird numerical issues with Stan as very few chains go wild for my model. I put weakly informative priors on my parameters and let Stan initialize randomly the chains. However, when I start off many chains (12) for my model I regularly end up with one of the chains not converging at all in almost all the parameters. This single chain then leads to completely wrong results as it is totally off the scale of all the other parameters. My question now would be how to handle this situation? I could set initial values to MLE parameter estimates of a simpler model and let Stan refine them, but then each chain would start off with the same initial parameters which does not sound very well to me. So, what strategies apply in this case?

Moreover, is it possible to drop a chain from the stanfit object? If so, how can I identify the invalid chain and then get rid of it?

My model itself is a mixed effect logistic regression of longitudinal data taken from subjects under different treatments.

Thank you very much for your help.

Best,

Sebastian

--
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]

On Wednesday, October 17, 2012 7:50:00 AM UTC-4, Sebastian Weber wrote:

I am experiencing rather weird numerical issues with Stan as very few chains go wild for my model. I put weakly informative priors on my parameters and let Stan initialize randomly the chains. However, when I start off many chains (12) for my model I regularly end up with one of the chains not converging at all in almost all the parameters. This single chain then leads to completely wrong results as it is totally off the scale of all the other parameters. My question now would be how to handle this situation? I could set initial values to MLE parameter estimates of a simpler model and let Stan refine them, but then each chain would start off with the same initial parameters which does not sound very well to me. So, what strategies apply in this case?

Moreover, is it possible to drop a chain from the stanfit object? If so, how can I identify the invalid chain and then get rid of it?

We have seen situations like this, even for simple models, but haven't figured out the cause yet. I don't think there is a convenient way to drop a chain from a stanfit object. However, if you have the sflist2stanfit() function from the thread on parallel chains, you can avoid including a problematic chain in a stanfit object. Just call stan() 12 times (in serial or in parallel) with chains = 1, mark the good chains, and then assemble a stanfit object, like

chains <- mclapply(...)
good <- c(1:2,4,5:9,11:12)
posterior <- sflist2stanfit(chains[good])

Ben

[Sebastian Weber]

Hi!

Hmm, thats a scary situation. So fixing the initial values did not solve the problem?

Sebastian

> --

[Shige]

Does running a very long chain help?

Shige

On Wednesday, October 17, 2012 7:50:00 AM UTC-4, Sebastian Weber wrote:

[Sebastian Weber]

Maybe, I don't know. I had 12 chains running 1000 iterations with 500 warmup steps. All but one chain equilibrated to their final values after 200 steps at most. The one chain which went off was very far from the rest for all parameters. The chain had a tendency towards the main bulk, but it was still after 1000 steps completely off, I think one would have to wait really long to fix this. Running the chain much longer becomes a bit inconvenient as this already took 4 hours or so. It is very odd to see all but one chain converge.

Sebastian

> --

[Ben Goodrich]

Until we figure out what the problem is, this is all speculation. But I don't think the starting values or the number of iterations are particularly relevant. If 11 out of 12 chains work, the 12th should be able to recover from whatever problem afflicted it, but that isn't happening for some reason.

[Daniel Lee]

Can you provide the model, data, and seed that causes the chain to go
off the rails? If there is a bug, perhaps we can hunt it down.


Daniel

> --
>
>

[Ben Goodrich]

On Wednesday, October 17, 2012 11:34:25 AM UTC-4, Daniel Lee wrote:

Can you provide the model, data, and seed that causes the chain to go
off the rails? If there is a bug, perhaps we can hunt it down.

https://groups.google.com/d/msg/stan-users/OjMSj1pITAQ/Tv8Pgy_88BoJ
https://groups.google.com/d/msg/stan-dev/6hRAFUT1Pdk/Ye5OngvwQ74J

Ben

[Bob Carpenter]

I have more notes below, but this may very well be
related to the patch Matt just made for initial adaptation.

We should release a new version soon with a bunch of
these fixes.

More inline below on specific issues.

> On Wednesday, October 17, 2012 11:11:43 AM UTC-4, Sebastian Weber wrote:
>
>
> Maybe, I don't know. I had 12 chains running 1000 iterations with 500 warmup steps. All but one chain equilibrated
> to their final values after 200 steps at most. The one chain which went off was very far from the rest for all
> parameters. The chain had a tendency towards the main bulk, but it was still after 1000 steps completely off, I
> think one would have to wait really long to fix this. Running the chain much longer becomes a bit inconvenient as
> this already took 4 hours or so. It is very odd to see all but one chain converge.

This isn't odd when you think about what's going
on in the implementation (some hints below).

> > On Wednesday, October 17, 2012 7:50:00 AM UTC-4, Sebastian Weber wrote:
> > Dear Stan users,
> >
> > I am experiencing rather weird numerical issues with Stan as very few chains go wild for my model. I put weakly
> informative priors

Stronger priors may not help if it's an init issue.

> on my parameters and let Stan initialize randomly the chains. However, when I start off many
> chains (12) for my model I regularly end up with one of the chains not converging at all in almost all the
> parameters. This single chain then leads to completely wrong results as it is totally off the scale of all the other
> parameters. My question now would be how to handle this situation? I could set initial values to MLE parameter
> estimates of a simpler model and let Stan refine them, but then each chain would start off with the same initial
> parameters which does not sound very well to me. So, what strategies apply in this case?
> >
> > Moreover, is it possible to drop a chain from the stanfit object? If so, how can I identify the invalid chain and
> then get rid of it?

I don't know about the R interface, but
the problem we have with this in general is that
just pruning chains affects the balance of
the algorithm. For instance, if problems arise
from two different modes, it's not technically wrong
to keep them both (though Stan would almost never
jump between modes with a big drop in log probability
in between). Also, if one chain is walking around in
a fat tail and another in a trunk, it can take a long
time (in theory and practice) for them to converge
so that R-hat looks OK.

On 10/17/12 11:26 AM, Ben Goodrich wrote:
> Until we figure out what the problem is, this is all speculation.
> But I don't think the starting values or the number of
> iterations are particularly relevant. If 11 out of 12 chains work,
> the 12th should be able to recover from whatever
> problem afflicted it, but that isn't happening for some reason.

This would be nice in principle, but floating-point arithmetic
is a harsh mistress. If an init causes issues like NaN or underflow
or overflow, it's essentially unrecoverable.

In the command-line code, we check if there are NaN values
produced and restart a chain with new random inits if there is.
But only during the first part of warmup. We should probably
think harder about this.

If arithmetic issues arise after inits, it'll just lead to
Metropolis rejection of all the subsequent samples (that is,
repeats of the same sample from an output perspective).

- Bob

[Ben Goodrich]

On Wednesday, October 17, 2012 12:14:32 PM UTC-4, Bob Carpenter wrote:

I have more notes below, but this may very well be
related to the patch Matt just made for initial adaptation.

Possibly related, but I am essentially testing with master and can replicate the users' problems.

If arithmetic issues arise after inits, it'll just lead to
Metropolis rejection of all the subsequent samples (that is,
repeats of the same sample from an output perspective).

 
I am still not understanding this part. If a proposal causes lp__ to be NaN or -infinity, how does that proposal get accepted? I was guessing that maybe lp__ is finite and good enough to be accepted with positive probability, but there is a numerical issue with its gradient.

Ben

[Sebastian Weber]

I would volunteer to test if the patches make this behavior going away. I already cloned the git repository, but I wasn't able to figure out on how to generate the rstan package which I use. Is there any help I overlooked to make that neat rstan package?

Concerning the model&data in this case - I used a large data set and when trying to make a smaller one the error disappeared for now. However, all my (larger) data sets I tried to fit so far failed. Larger means several thousand terms in the likelihood.

Best,

Sebastian

[Jiqiang Guo]

See below for how to build/install RStan after you have git cloned the stan repo. 

On Wed, Oct 17, 2012 at 12:39 PM, Sebastian Weber <sebasti...@frias.uni-freiburg.de> wrote:

I would volunteer to test if the patches make this behavior going away. I already cloned the git repository, but I wasn't able to figure out on how to generate the rstan package which I use. Is there any help I overlooked to make that neat rstan package?

Suppose STAN_HOME is the directory where you have stan.  Go to STAN_HOME/rstan, where you can do 

$ make install 

to build and install the rstan package.  Currently, the makefile works well only under Unix-like OS such as Linux and Mac. 

--

Jiqiang

 

--

[Bob Carpenter]

On 10/17/12 12:37 PM, Ben Goodrich wrote:
> On Wednesday, October 17, 2012 12:14:32 PM UTC-4, Bob Carpenter wrote:
>
> I have more notes below, but this may very well be
> related to the patch Matt just made for initial adaptation.
>
>
> Possibly related, but I am essentially testing with master and can replicate the users' problems.
>
> If arithmetic issues arise after inits, it'll just lead to
> Metropolis rejection of all the subsequent samples (that is,
> repeats of the same sample from an output perspective).
>
>
> I am still not understanding this part. If a proposal causes lp__ to be NaN or -infinity, how does that proposal get
> accepted?

It won't ever be accepted. The problem's if it's
an init.

> I was guessing that maybe lp__ is finite and good enough to be accepted with positive probability, but there
> is a numerical issue with its gradient.

This could cause the samples not to move very slowly or
not at all. This can be an issue at points near singularities
or undifferentiable points (like caused by fabs()).

- Bob

[Ben Goodrich]

On Wednesday, October 17, 2012 12:46:25 PM UTC-4, Jiqiang Guo wrote:

See below for how to build/install RStan after you have git cloned the stan repo. 

On Wed, Oct 17, 2012 at 12:39 PM, Sebastian Weber <sebasti...@frias.uni-freiburg.de> wrote:

I would volunteer to test if the patches make this behavior going away. I already cloned the git repository, but I wasn't able to figure out on how to generate the rstan package which I use. Is there any help I overlooked to make that neat rstan package?

Suppose STAN_HOME is the directory where you have stan.  Go to STAN_HOME/rstan, where you can do 

$ make install 

to build and install the rstan package.  Currently, the makefile works well only under Unix-like OS such as Linux and Mac. 

 
And it takes a long time. On a multicore machine, you can do

export MAKEFLAGS=-j

or

export MAKEFLAGS=-j8

before make install to parallelize parts of the install to build the help pages / build / install process.

Ben

[Ben Goodrich]

In the example on stan-dev, the model is just Gaussian with unknown mean and standard deviation. So, I am doubting that it is a problem with inits or singularities.

Ben

[Michael Betancourt]

As Bob says, far away from the typical set you can get extreme gradients that cause the numerical integration to diverge and yield mostly rejections.

The same behavior could also be indicative of p(typical point -> point) << 1. Because of detailed balance p(point -> typical point) will also be much less than one: once you have the misfortune of transitioning to that point it will take forever to transition back to the typical set. The more chains you run, and the longer you run them, the more likely you are to hit one of these unlikely regions. Note that, if we reject on NaNs and infinities, then the first behavior can also be thought of in these terms.

> --
>
>

[Daniel Lee]

Ben, can you help me figure out which seeds are causing these issues?
I see the problem by running rstan, I've pulled the seeds using
get_seed(chain[[3]]), but I'm not seeing the problem from the command
line.

[Sebastian Weber]

Hi!

To report back - I now start stan with seeds "0" as indicated in the stan documentation of rstan and the problem with chains going into the wild is gone. So in my case, the seeds seem to be an issue. I am still using rstan 1.0.2. Next, I will try what happens if I let R generate the seeds. I hope that works.

Is there already a planned release date for 1.0.3? It seems there are a lot of urgent fixes in the upcoming release.

Best,

Sebastian

[Shige Song]

Hi Sebastian,

I could not find in stan manual the instruction of setting the random seed to "0". Were you referring to Modeling Language User's Guide and Reference Manual?

Thanks.

Shige

--

[Sebastian Weber]

Hi!

No, I referred to the function call of stan in the rstan package:

init
One of digit 0, string "0" or "random", a function that returns a list, or a list of initial parameter values with which to indicate how the initial values of parameters are specified. "0": initialize all to be zero on the unconstrained support; "random": randomly generated; list: a list of lists equal in length to the number of chains (parameter chains), where each list in the list of lists specifies the initial values of parameters by name for the corresponding chain. function: a function that returns a list for specifying the initial values of parameters for a chain. The function can take an optional parameter chain_id.


Just type in ?stan in R.

Best,

Sebastian

[Ben Goodrich]

On Thursday, October 18, 2012 8:23:16 AM UTC-4, Sebastian Weber wrote:

To report back - I now start stan with seeds "0" as indicated in the stan documentation of rstan and the problem with chains going into the wild is gone. So in my case, the seeds seem to be an issue. I am still using rstan 1.0.2. Next, I will try what happens if I let R generate the seeds. I hope that works.

Is there already a planned release date for 1.0.3? It seems there are a lot of urgent fixes in the upcoming release.

The problem you were having is probably fixed if you make an rstan from the master branch of the git repo. But not everything is fixed yet.

Ben

[Sebastian Weber]

Hmm,

I don't mind running bleeding edge code, but it always feels a bit scary, especially if one wants to publish results. I mean, when going to press I actually prefer to have a deterministic parameter estimation technique which I get when I fix version numbers and seeds. Ok, report a git hash code as version number does the same thing, but this is not really nice…

Best,

Sebastian

[Ben Goodrich]

On Thursday, October 18, 2012 10:38:54 AM UTC-4, Sebastian Weber wrote:

I don't mind running bleeding edge code, but it always feels a bit scary, especially if one wants to publish results. I mean, when going to press I actually prefer to have a deterministic parameter estimation technique which I get when I fix version numbers and seeds. Ok, report a git hash code as version number does the same thing, but this is not really nice…

I agree. I was just implying that we probably need to fix the other issues before we release v1.0.3. But it seems to be a problem with rstan only, so if you need to do something this week, it should work fine with command-line Stan v1.0.2.

Ben

[Eric Brown]

In the general case, it may be beneficial to report (or have available) the initial parameter set instead of just the seed.  While the process may not be 100% deterministic (due to architecture, compiler, or optimization differences not to mention different versions of Stan), ideally a good initial parameter set would be stable for any version on Stan.

 

-Eric

[Ben Goodrich]

On Thursday, October 18, 2012 11:35:21 AM UTC-4, Eric Brown wrote:

In the general case, it may be beneficial to report (or have available) the initial parameter set instead of just the seed.  While the process may not be 100% deterministic (due to architecture, compiler, or optimization differences not to mention different versions of Stan), ideally a good initial parameter set would be stable for any version on Stan.

The inits are always available in R.

inits <- get_inits(stanfit_object)

Ben

[Bob Carpenter]

On 10/18/12 11:35 AM, Eric Brown wrote:
> In the general case, it may be beneficial to report (or have available) the initial parameter set instead of just the
> seed.

I think that's a good idea. We should make it the
first line of the returned output in the CSV file from
the command line. It is available in R already, but
I don't think it is included in the warmup samples
directly.

> While the process may not be 100% deterministic (due to architecture, compiler, or optimization differences not
> to mention different versions of Stan), ideally a good initial parameter set would be stable for any version on Stan.

Yes, that shouldn't change. If anything, we should
get more robust to crazier parameter inits as we go.
We regression test model behavior (in the CS, not stats sense)
continuously on our build server.

> On Thursday, October 18, 2012 9:38:54 AM UTC-5, Sebastian Weber wrote:
>
> I don't mind running bleeding edge code, but it always feels a bit scary, especially if one wants to publish
> results. I mean, when going to press I actually prefer to have a deterministic parameter estimation technique which
> I get when I fix version numbers and seeds. Ok, report a git hash code as version number does the same thing, but

> this is not really nice�

I very much like having at least this level of reproducibility.
Ideally, I'd take the data release (in whatever compressed
format I get it in), have a one-touch script to munge it to
whatever I need, and a one-touch script to generate all the
results. I rarely achieve quite that level of reproducibility
other than for demos.

I think it's easier on everyone to use "official" versions.
We should be pushing out new versions fairly regularly.

There's no guarantee we'll keep our Git hash code numbers
indefinitely, but we should have all the tarball releases.

- Bob