What is the largest model you've used with Stan?

[Arthur B.]

Hello group,

Obviously the mere dimensionality of a problem doesn't fully determine the speed of convergence of a MCMC. The correlation between the variables, the complexity of the model, etc all play important roles.

That said, I'm trying to get a sense of the order of magnitude of the size of problems that have a chance to be tackled.

Is it sensible to use Stan on a model with hundreds of parameters? Tens of thousands? Millions?

Some of the default behaviors of cmdstan (such as printing each parameter for each sample) suggest that the typical models have at most dozens of parameters. Is that so?

Thanks,
Arthur

[Bob Carpenter]

The examples are misleading because most of them are from
textbooks that have kept things simple and small.

We've fit simple models like IRT with millions of parameters,
but it takes a few hours. We regularly fit all kinds of
models with thousands of parameters.

CmdStan saves each sample (but can be set to thin). We're
thinking about adding something to CmdStan like RStan has
that will let you suppress output for some variables.

If the saved output is too big, it can choke our posterior
analysis tools. But that's a different problem than the
sampler not working. You can still easily do things like
compute posterior means and sd using online algorithms like
Welford's and you can do pretty well at getting wide
intervals (50% intervals require saving half the samples online,
whereas 5% intervals only require saving 5% in memory).

If you want to do logistic regression with L1 regularization
and millions of parameters and billions of sparse training instances,
you'll need to use something more specialized like Vowpal Wabbit.
Of course, there's a problem with summarizing output given their
hashing, but again, that's a different problem than being able to
fit the model.

- 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.
> For more options, visit https://groups.google.com/d/optout.
>

[Richard McElreath]

I think the largest number of parameters I've used in a Stan model is about 30k. It worked just fine.

In fact, sometimes I add parameters by taking out e.g. marginalization and it actually fits better. Earlier this year, I worked a hidden Markov model with gamma latency emissions with lots of censoring, and just estimating the censored values as missing with known constraint lead to slightly better mixing than using the gamma ccdf, even though it added lots of parameters to the posterior.

[Luc Coffeng]

"and just estimating the censored values as missing with known constraint lead to slightly better mixing than using the gamma ccdf, even though it added lots of parameters to the posterior."

Hi Richard,

What do you mean by the above? What do you increment the log posterior with when you say "estimate as missing"?

Luc

[Richard McElreath]

On Friday, September 5, 2014 1:10:59 PM UTC-7, Luc Coffeng wrote:

What do you mean by the above? What do you increment the log posterior with when you say "estimate as missing"?

I just checked the Stan manual to see what the team says about this, and luckily section 8.3 (page 65) explains exactly what I mean. (The Stan manual is really becoming quite a nice resource---it's improving faster than I can keep up with it.)

I was working with some data in which I believe gamma ccdf was numerically finicky, so the missing value approach was working a little better. And as a bonus it gives estimates for each of the partially observed intervals.

[Arthur B.]

Thanks all for the answers.
I was checking the Stan manual and it says:  "The last tweak to Stan before version 1.0 was to estimate a diagonal matrix during warmup; this has since been upgraded to a full mass matrix in version 1.2."
How exactly is Stan going to handle millions of parameters if it attempts to build (worse, invert) such a matrix?
 

[Bob Carpenter]

The default for HMC is to estimate a diagonal mass
matrix during warmup. You can override that and have
it estimate a full mass matrix or just take a unit
mass matrix. The default for optimization is to use
L-BFGS which approximates a Hessian using a finite history
of gradients, but never needs to store the whole matrix or
do an inversion.

We haven't found many problems where a dense mass matrix
is helpful.

Feel free to ask on the list about how things work now --- some
of the historical notes are now just that.

- Bob

[Arthur B.]

Thanks! Here's one more question: is there an easy way to use a MAP estimate as a starting point for a chain?

[Bob Carpenter]

It's possible using optimization. It's easier in R/Python than
from the command line.

Having said that, Stan usually converges quickly from diffuse
starting points. Running multiple such chains and getting to
the same end point is a further check that you're really
at convergence.

Also, many models of interest, such as hierarchical models,
do not have posterior modes (density approaches infinity as
hierarchical variance approaches zero and the low-level parameters
approach the hierarchical mean). In those cases, it's possible
to only iterate the L-BFGS solver for optimization a few iterations,
but you're in uncharted waters. Do report back if you go there
and survive and have advice for the rest of us.

- Bob

On Sep 5, 2014, at 9:47 PM, Arthur B. <arthur....@gmail.com> wrote:

> Thanks! Here's one more question: is there an easy way to use a MAP estimate as a starting point for a chain?
>

[Luc Coffeng]

Hi Richard,

Thanks for clarifying. Do you have any idea why adding parameters for censored values yielded better mixing, compared to integrating out the censored values? Could it be a results of potential differences in numerical stability of the PDF and CDF functions in Stan? E.g., the CDF function normal_cdf_log will default to -Inf for (y-mu)/sigma < 8.25, which may happen for censored values if the model's expected value for said value is well above the censoring threshold (assuming left-censoring in this case) and sigma is not too large. The manual does not state anything about the numerical stability of normal_log, but if it's more stable than the CDF (i.e. it defaults to -Inf for z << 8.25), it might explain things.

Or maybe there's an entirely different explanation for the better mixing you report?

I've been having mixing issues with a model for data that are partly left-censored. In this model I ise the lognormal_cdf_log function to integrate out censored values.. It'd be great if using the PDF (effectively truncated by the constraints on the parameter) could help there, though it will be pain to code as the observations in the data have been censored at different values (hopefully we'll be able to use vectors of values for parameters constraints in the near future).

Thanks,

Luc

[Richard McElreath]

Luc:

I think the issue I was having with gamma_ccdf was numerical instability with some very long intervals. I recall seeing a lot of Metropolis rejections. Honestly, I never figured out exactly what the issue was.

I happen to have a script that I used to validate the missing value approach with variable observation windows. Attaching it to this message. It simulates gamma distributed latencies with variable censoring windows. Then there's a Stan model that uses transformed parameters to build the imputations. Hopefully it makes sense.

-Richard

[Luc Coffeng]

Thanks for sharing, Richard. The code made sense to me. :)