simple stan float

[Andrew Gelman]

Hi all. The following causes stan to float. It's a stupid example, a model that's empty! But I'm surprised that it floats instead of just running (as a random-walk drift) or reporting an error.

Is it possible that it's trying to do infinite steps in the Nuts trajectories? I thought we'd capped the trajectory lenghth at some number like 256.

Here's the R output:

> simple.fit <- stan (model.code=simple.code, n.iter=100, n.chains=4)

COMPILING MODEL 'anon_model' FROM Stan CODE TO C++ CODE NOW.
COMPILING THE C++ CODE FOR MODEL 'anon_model' NOW.
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).


[and then it floats and I have to "force quite" to Rstudio]

And here's the model:

simple.code <- "
data {
}
parameters {
real a;
}
model {
}
"



By the way, it runs fine if I simply add the line
a ~ normal (0, 1)
to the model.

A

[Bob Carpenter]

> simple.code <- "
> data {
> }
> parameters {
> real a;
> }
> model {
> }

We could put some kind of timing monitor on things
perhaps to kill after a certain amount of time. But
that'll probably require dealing with system-specific
thread libs or a wrapper like Boost. Unless we can do
it in R. The problem would be a default value, as we
don't want to kill long-running models.

Yes, there's an upper limit on steps in the leapfrog
algorithm. It's probably hanging before that during
adaptation. So we should probably cap that somewhere,
too and return with an error message.

Matt?

- Bob

[Andrew Gelman]

If there's an upper limit on L, then the problem of Stan floating must be something else, no? I imagine it is some infinite loop somewhere in Stan? Or are you saying that L is bounded in the production runs but not in the warmup stage?

[Matt Hoffman]

It's hanging during find_reasonable_parameters(), which looks for the
smallest step size that results in poor accuracy and then dials that
back. The problem is that for this empty model *any* step size will
result in perfect accuracy. So it never quits.

I've pushed a change that should address this problem—it checks for
huge or zero step sizes during find_reasonable_parameters() and throws
std::runtime_errors with the message "Posterior is improper. Please
check your model." if the step size goes above 1e300 or "No acceptably
small step size could be found. Perhaps the posterior is not
continuous?" if the step size goes to 0.

Feel free to change the wording.

Matt

PS: We don't have to do guesswork to figure out where the program is
hanging—a quick profiler run was all it took to isolate the problem.

[Andrew Gelman]

Matt:

Don't forget also the mystery of Stan for the 8 schools apparently running less efficiently when it has more time to adapt in the warmup pahse.
A

[Matt Hoffman]

One mystery at a time! "It stops running" is an easier symptom to make
diagnoses based on than "it runs slower than it seems like it should."

[Bob Carpenter]

On 8/10/12 5:07 PM, Matt Hoffman wrote:
> One mystery at a time! "It stops running" is an easier symptom to make
> diagnoses based on than "it runs slower than it seems like it should."

Andrew was measuring in terms of effective sample size per
number of iterations. We really want to measure effective sample
size per unit of time.

Speaking of time, did you (Jiqiang, Ben) figure out how
to compile RStan with -O3.

- Bob

[Jiqiang Guo]

The shared part distributed with rstan (libstan.*) is compiled with -O3 for linux and mac, but not for windows. But the model is compiled with the options R has as the default, which might not be -O3, perhaps -O2.  For now, I have no control on that at all and I can only tell the user to modify $R.home()/etc/.../Makeconf to change that to -O3. 

Jiqiang 

[Andrew Gelman]

Matt:
As a reminder, here's what I found:

n.iter ESS
1000 55
10000 450
100000 990
1000000 7500

The inferences all look basically reasonable

Does the sublinear growth of ESS as a function of n.iter suggest that the Nuts tuning is actually getting _worse_ when we have more time in the warmup phase? Giving Nuts enough rope to hang itself??

A

[Matt Hoffman]

It's not impossible. Another possibility is that our ESS calculations
are optimistic, although that seems less likely if we're using
multiple chains to estimate the relevant means and variances. I'll
look into it.

[Andrew Gelman]

Are we currently using our preferred ESS rules?

[Jiqiang Guo]

I do not know what is preferred.  But I did not change the implementation in mcmc::chains. 

Jiqiang 

[Andrew Gelman]

What's preferred, I think, is what's described in the new BDA, see attached.  I thikn it gets implented using FFT.

It's possible that my description in the (forthcoming) book is out of date and that we're doing something better.  If so, please let me know.  Thx

A

[Bob Carpenter]

On 8/10/12 5:41 PM, Andrew Gelman wrote:
> Are we currently using our preferred ESS rules?

We're using the var+ estimator of variance.

The implementation could be sped up by resizing
and indexing rather than using push_back.

We're just putting interrupts on the to-do list,
which is growing faster than our programming
resources. After the first release, for which I
don't want to add another thin mint of functionality
unless we absolutely have to, I think we should sit
back and prioritize the roughly 200 or 300 things
on our to-do list.

- Bob

[Andrew Gelman]

Does that mean we are using (12.8) on page 355 of the document I just sent?

[Bob Carpenter]

The BDA draft was the basis of our implementation,
which uses discrete FFT to compute the autocovariances
on a per-chain (not split) basis and uses the var+ estimate
that's also used in r-hat calcs to normalize.

You and Kenny said what you did may need to be
revised. When that happens, you'll need to let us know
so we can recode it!

- Bob

> On Aug 10, 2012, at 11:45 PM, Jiqiang Guo wrote:
>
>> I do not know what is preferred. But I did not change the implementation in mcmc::chains.
>>
>> Jiqiang

[Bob Carpenter]

On 8/10/12 5:56 PM, Andrew Gelman wrote:
> Does that mean we are using (12.8) on page 355 of the document I just sent?

Yes, that's what everyone uses.

The key is (12.7), which uses the multi-chain var+ estimator
rather than the per-chain variance estimate.

Ack. I just realized the latest manual is missing page
numbers. LaTeX is such a pain. I tried to get page
numbers off part dividers, but it's a huge hassle.

- Bob