convergence in stan

770 views
Skip to first unread message

Linas Mockus

unread,
May 17, 2014, 8:55:17 AM5/17/14
to stan-...@googlegroups.com
Hi,

I have a question about convergence in stan. Rhat of a fit is 1 and looks like mixing is good (warmup=20000, n.draws=10000), but there are regions of almost straight lines between the nicely mixing regions (see the attached, page 1).  Those regions with nearly straight lines create humps in the posterior (page 3).  The predictive posterior (green lines .25/.75 quantiles and median blue line on page 4 looks reasonable.  Red line is MAP).  In order to find the MAP Differential Evolution optimization procedure was needed - 2 compartmental PK model is nonlinear and had multiple modes.  In this model I have used multiplicative error (the stan code attached).  For complete picture r code is attached as well.

Can anybody explain/suggest what is going on?

Thanks a lot,
Linas
population.pdf
2comp_multerr_h.stan
2comp_multerr_h.r

Michael Betancourt

unread,
May 17, 2014, 10:52:30 AM5/17/14
to stan-...@googlegroups.com
You’ll want to focus your attention on the sigma_theta trace — everything stalls when
sigma_theta hits zero.  What’s happening here is that HMC (and just about everything)
is not geometrically ergodic for your model, and getting stuck like this is a manifestation 
of the lack of geometric ergodicity.  You’ll also note that n_divergent becomes non-zero
in this region as well, which is an HMC specific manifestation of the problem.

This is particularly important as geometric ergodicity is what guarantees that Monte Carlo 
estimates converge to the true expectations (via the Monte Carlo Central Limit Theorem), 
and in a case like this you have to be very dubious of the estimates.  Note that algorithms
like Random Walk Metropolis and the Gibbs sampler suffer from the same problem, only
they’re usually not fast enough to be run for enough iterations to see the pathology manifest.

For more on the problem and possible solutions see http://arxiv.org/abs/1312.0906.

--
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.
<population.pdf><2comp_multerr_h.stan><2comp_multerr_h.r>

Linas Mockus

unread,
May 21, 2014, 10:20:09 AM5/21/14
to stan-...@googlegroups.com
Michael, thanks for the insight.  I have used theta=theta_sigma*theta_var+theta_mu and now the model runs very smoothly.

The model is a hierarchical 2 compartmental PK model and data includes blood plasma measurements over some period of time for many individuals.  I wonder if I were to apply non-hierarchical version of the same model to single individual (it has only 8 data points) how many parameters the model can handle?  The 2 compartmental model has 6 parameters including error term.  Can I use the same trick which introduces one more parameter (note that I have only 8 data points) and still get convergence?
Thank you,
Linas
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/jne6RL_kDC4/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

Michael Betancourt

unread,
May 21, 2014, 12:08:52 PM5/21/14
to stan-...@googlegroups.com
>
> The model is a hierarchical 2 compartmental PK model and data includes blood plasma measurements over some period of time for many individuals. I wonder if I were to apply non-hierarchical version of the same model to single individual (it has only 8 data points) how many parameters the model can handle? The 2 compartmental model has 6 parameters including error term. Can I use the same trick which introduces one more parameter (note that I have only 8 data points) and still get convergence?

What do you mean by non-hierarchical? Are you taking the hyper parameter posterior from the
first fit and using that as a prior for the new individual? So long as your prior is well-behaved you
can handle as little data as possible. Also, in terms of parameterizations the non-centered
parameterization (which sped up your model) works best with sparse data so would continue
to be the best choice for a single individual.

Reply all
Reply to author
Forward
0 new messages