What happens before the first sample?

[Julian King]

Hi all,

I have a fairly complicated regression-type model which runs slowly but fine (500 iterations is a few hours). The model has N of about 20,000.

Is there anything that STAN does between when it says "Now sampling..." and when the first sample is generated that is quite computationally demanding?  The model takes more than an hour to generate the first sample, then proceeds at a predictable pace.

I'm curious as to what happens here that could take so long.

Thanks in advance

[Bob Carpenter]

First, the data's read in and the transformed data block is
executed and all constraints are validated. That's about it.

I believe RStan also allocates memory to store the draws
at that point, the speed of which will depend on how
you're running it. Are you sure you're not also waiting
for the first 10% of iterations? The warmup iterations
can be slow when they start adapting but haven't found
the high probability mass volume of the posterior.

- Bob

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

[Julian King]

I have refresh=1 in the call to Stan, so I see what's happening sample-by-sample - so I am seeing the first iteration.

I do have a fair bit of stuff in the transformed data block, but it's basically just mapping input parameters into vectors or matrices for later use, and some standardisation of data (dividing by standard deviation etc).

I don't have many initial parameter constraints - all my parameters are basically defined on (-infinity,infinity) or (0,infinity).  I have priors on everything.   My model does have exponential transformations of predictors though.... which might end up being rounded to infinity for poor initial starting choices.

Could it be that the model is trying to generate parameters which give an infinity somewhere, and is then trying again lots of times until it gets something meaningful?

Thanks

[Bob Carpenter]

It can do up to 100 initializations before giving up when
it can't find an initialization that works. So the first
iteration can do 100 evaluations of the log density and
gradient even if it doesn't go anywhere. Of course, it might
do 1000 evaluations in one iteration under default parameters, too --- depends on the
model and inits. So you should see how many leapfrog steps
were done in that step.

Also, make sure that every value of variables that satisfies
the declared constraints has finite log density in the model.
If not, init has to do a lot of guessing and then sampling can
devolve to a random walk because of boundary rejections.

- Bob

[Julian King]

I've checked with a print() statement that it's not the transformed parameters block.

The parameter constraints should be easily satisfied (my parameters are all defined as reals or real<lower=0)s), but it's possible that the initial values produce infinities somewhere because I have a Poisson model with a bunch of exp terms all over the place (could be getting underflow or overflow somewhere...)

Thanks for the suggestions - I will try to see if I can identify the problem through elimination.  I might try passing better inits in and seeing if that helps (at the moment I'm just going with the default)

[Bob Carpenter]

If you have those interval constraints on parameters, remember
the default inits are around the middle of the range, which can
be very problematic.

- Bob

[Julian King]

I narrowed this down to something to do with the large number of stored samples being recorded - I had lots of n ~ 20,000 vectors in the transformed parameters block and the generated quantities block that were being recorded.  Putting these in a local block (ie not saving their values) massively cut down on the initialisation time.

I guess this might be something to do with memory allocation, although the total memory allocation only ends up being 2 GB or so.  No idea why it would take an hour to allocate 2 GB of memory.

[Bob Carpenter]

Is there a reproducible example? This is an RStan (or PyStan)
issue --- there's no preallocation in CmdStan.

- Bob