Non-linear model: horribly mixing chains since rstan v2.1

[Luc Coffeng]

Dear Stan-folks,

I've been tinkering at some code for a non-linear regression model (which I have posted about before), and have gotten stuck with horribly mixing chains. This problem has appeared since I upgraded to rstan 2.1.0 and 2.2.0.

Sometime in August 2013, I got the model running successfully and managed to fit it to a dataset with the latest version of rstan available then (can't tell the version from what I've saved). Back then, not all distributions had yet been vectorized, so upon the release of rstan 2.1.0 I thought I might try and vectorize everything and see if things would speed up. Coding up the vectorization was no problem, but as it turned out, since rstan 2.1.0 (I've tried 2.2.0 as well) the sampler simply does not seem to converge to the right mode anymore, churns out a handfull of effective sample sizes (N < 20) no matter how much brute force applied, and is   v e r y   slow at sampling. When sequentially sampling multiple chains, upon completion of the first chain, R sometimes even freezes completely (even with as few as 1K iterations). Reverting back to the unvectorized code did not change these findings. I have also tried implementing the "Matt trick" for the log-normally distributed random effects (defining unbounded raw random effects, sampling these from standard normal distribution, and transforming these in the transformed parameters block into lognormally distributed random effects - no correction for posterior curvature required, right?) - the sampler either slowed down to a grinding halt, or spitted out an endless stream of red Metropolis rejections.

Anyhow, to start afresh, I've coded up a script that simulates data such as the model would predict it, and I have tried to fit the model to these simulated data. The result remains horrible mixing. Any advice that you may have would be highly appreciated. I'm about to bash my head into a wall. For your convenience, I've attached:

1. Script to simulate data and fit the completely vectorized model to simulated data. The choice of using full mass matrix was based on my earlier impression that it really helped effective sampling back in August 2013. Mind the potential warning in the middle of the script about censored observations after first positive observations.
2. PDF with fit of posterior to subset of real data, August 2013 (some version of rstan before 2.0.0); to give you some idea of what things should be looking like. :) [Edit: the dashed line is the detection limit for the response; observations under the detection limit are censored in the model, i.e. integrated out using a cumulative log-normal distribution).]
   

And via google drive:

1. PDF with visualization of the posteriors based on the real data, February 2014 (with rstan 2.2.0, viz produced with package ggmcmc). :(
2. PDF with visualization of the posteriors based on the real data, August 2013 (some version of rstan before 2.0.0, viz produced with package ggmcmc). :)

Big fat thank you in advance to anyone who wants to chip in.

Luc

P.s.: I'm working Mac OSX 10.7.5, R 3.0.2, rstan 2.2.0, Rcpp 0.10.6

sysname "Darwin", release "11.4.2", version "Darwin Kernel Version 11.4.2: Thu Aug 23 16:25:48 PDT 2012; root:xnu-1699.32.7~1/RELEASE_X86_64"

 

[Bob Carpenter]

The code was largely refactored as of 2.x, but the goal was to maintain
the same sampling strategy. Maybe Michael will have some more suggestions.

The only major difference you should find in 2.x versions (not
counting the buggy 2.0.y versions) are the larger acceptance
rate target (2.1.0) and initial step-size tuning based on the
actual initialization point rather than 0.

I don't know how much the full mass matrix estimation changed in
2.x. But you might try turning that off.

In the print() of the results, what is the tree depth and acceptance
rate you're getting?

How are you initializing? I would try init=0 or some value
of init=z where 0 < z << 2, say z = 0.25. Or even better,
start at some realistic values, perhaps with some noise.

The other thing you can do is try turning the accept rate
target (delta) back down to the old default, 0.65. A larger
target should mean more precise sampling, so an issue might
be that the new version of Stan is actually doing a better job
of sampling --- that is, you might actually have the kind of
multimodal posteriors it looks like you're getting.

Mac OS X Darwin (10.7) is pretty far behind in C++ compilers,
but I don't think that'd be an issue if you managed to compile
and do so with optimization (set_cppo("fast") was set in previous
version before you compiled next version and left on after compile).

The model itself is too complicated to comprehend. The main thing
you have to do is make sure the model has support over any variables
that meet all the declared constraints.

The CDFs are slow, but they should have always been slow.

- Bob

On 2/17/14, 5:59 AM, Luc Coffeng wrote:
> Dear Stan-folks,
>
> I've been tinkering at some code for a non-linear regression model (which I have posted about before

> <https://groups.google.com/d/topic/stan-users/GBLbw8dftec/discussion>), and have gotten stuck with horribly mixing

> chains. This problem has appeared since I upgraded to rstan 2.1.0 and 2.2.0.
>
> Sometime in August 2013, I got the model running successfully and managed to fit it to a dataset with the latest version
> of rstan available then (can't tell the version from what I've saved). Back then, not all distributions had yet been
> vectorized, so upon the release of rstan 2.1.0 I thought I might try and vectorize everything and see if things would
> speed up. Coding up the vectorization was no problem, but as it turned out, since rstan 2.1.0 (I've tried 2.2.0 as well)
> the sampler simply does not seem to converge to the right mode anymore, churns out a handfull of effective sample sizes
> (N < 20) no matter how much brute force applied, and is v e r y slow at sampling. When sequentially sampling
> multiple chains, upon completion of the first chain, R sometimes even freezes completely (even with as few as 1K
> iterations). Reverting back to the unvectorized code did not change these findings. I have also tried implementing the
> "Matt trick" for the log-normally distributed random effects (defining unbounded raw random effects, sampling these from
> standard normal distribution, and transforming these in the transformed parameters block into lognormally distributed
> random effects - no correction for posterior curvature required, right?) - the sampler either slowed down to a grinding
> halt, or spitted out an endless stream of red Metropolis rejections.
>
> Anyhow, to start afresh, I've coded up a script that simulates data such as the model would predict it, and I have tried
> to fit the model to these simulated data. The result remains horrible mixing. Any advice that you may have would be
> highly appreciated. I'm about to bash my head into a wall. For your convenience, I've attached:
>

> 1. Script to simulate data and fit the completely vectorized model to simulated data. The choice of using full mass

> matrix was based on my earlier impression that it really helped effective sampling back in August 2013. Mind the
> potential warning in the middle of the script about censored observations after first positive observations.

> 2. PDF with fit of posterior to subset of real data, August 2013 (some version of rstan before 2.0.0); to give you some

> idea of what things should be looking like. :)
>
>

> And via google drive:
>
> 1. PDF <https://drive.google.com/file/d/0BxTB9IYjuodjQ1BfV2Rzb0huSDA/edit?usp=sharing> with visualization of the
> posteriors based on the real data, February 2014 (with rstan 2.2.0, viz produced with package /ggmcmc/). :(
> 2. PDF <https://drive.google.com/file/d/0BxTB9IYjuodjejh6LWxhc3FTVjg/edit?usp=sharing> with visualization of the

> posteriors based on the real data, August 2013 (some version of rstan before 2.0.0, viz produced with package

> /ggmcmc/). :)

>
>
> Big fat thank you in advance to anyone who wants to chip in.
>
> Luc
>
> P.s.: I'm working Mac OSX 10.7.5, R 3.0.2, rstan 2.2.0, Rcpp 0.10.6
>
> sysname "Darwin", release "11.4.2", version "Darwin Kernel Version 11.4.2: Thu Aug 23 16:25:48 PDT 2012;
> root:xnu-1699.32.7~1/RELEASE_X86_64"
>

> --
> 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/groups/opt_out.

[Bob Carpenter]

I hope Michael has a chance to weigh in on what might be going
on with the sampler that's different between Stan 1.x and 2.x.

It may just be the new sampler's doing a better job of exploring
the posterior (he says hopefully). Two questions on that:

* Do you know if the earlier version of Stan explored across what
now look like separated modes?

* Does the probabilty mass around the modes look similar in the
current exploration?

If your tree depth's going that deep, then you probably don't have
a problem with support --- they'd almost certainly kick out early
if that were the case.

As to initialization, it works in RStan like in
R2WinBUGS. Here's the doc from the stan() command in R:

list: a list of lists equal in length to the number of chains (parameter chains),
where each named list (an element of the list of lists) specifies the initial values
for parameters by names for a 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 for which
the chain_id if specified or the integers from 1 to chains will be supplied respectively
to the function for generating initial values. See the examples below of defining such
functions and using a list of lists for specify initial values. Additionally, see below notes
on that RStan treats a vector of length 1 as a scalar.

So the simplest thing to do if you have the list of inits init_vals, is to

init_fun <- function(x) { return(init_vals); }

and then in the call to stan(), use

init=init_fun

The inits are on the constrained scale as specified in the model.

- Bob

On 2/17/14, 7:44 PM, Luc Coffeng wrote:
> Dear Bob,
>
> Thanks for your input! Apologies for the lack of model introduction on my part.
>
> I'm working on a hierarchical model for blood parasite levels in experimentally infected subjects, assuming a log-normal
> distribution for blood parasite levels. The blood parasites appear in cycles of generations (called waves in the model,
> parameter T representing the timing of the center of the first wave), each cycle being modeled by two scaled normal
> cdf's: one for the build-up of parasites, another for the disappearance of parasites. Both CDF's are scaled by parameter
> M, which represents the amplitude of the first parasite wave. Due to parasite reproduction (this happens when parasites
> are not visible in the blood), consecutive parasite waves have increasing amplitude (multiplication factor represented
> by parameter R), resulting in a typical pattern (see attached pdf). Timing (T), amplitude (M), and multiplication rate
> (R) are all assumed to vary between persons. To anchor the model, I'm explicitly modeling censored observations, and I
> assume that all observations up to the first observation above the detection limit are samples of first generation of
> blood parasites.
>
> *Mass matrix*
> The full mass matrix doesn't seem to be the problem, apart from taking longer to complete. With a diagonal mass matrix,
> N_eff is still low (<20 for 1000 samples in 60 seconds). Applying brute force with 20000 samples (and thin = 10, warmup
> = 2000) doesn't help.
>
> *Acceptance and treedepth*
> With diagonal mass matrix, 2000 warmup, 4000 samples, init = 0, and 0.65 target acceptance rate:
>
> Tree depth:
> -1 0 1 2 3 4 5 6 7 8 9
> 403 643 581 215 63 23 38 20 10 3 1
>
> Average acceptance rate was ~65%, as targeted.
>
> *Multimodal posterior*
> Indeed, for some parameters the posterior seems to be multimodal. However, the multimodality of the posterior results in
> this case seems to be due to chains getting stuck in different areas of the posterior. Somehow, I need to get the
> sampler to behave more flexibly, and my impression is that it was able to explore the posterior much more efficiently
> before; I used to get nice caterpillar plots and chains used to converge to the same area of the posterior. Now, the
> traceplots are very jagged and chains end up in different places :(.
>
> *Initializing parameter values*
> So far, I tried initializing with random inits or init = 0 (before rstan 2.x, I let rstan initialize itself). Today, I
> tried feeding the model a list of inits (and lists of lists of inits for multiple chains), using the actual parameter
> values that I used to simulate the data (also added this to the attached r script):
>
> gen.init <- function() {
> list(
> log_T1 = log_T1_mu,
> S_T1 = log_T1_sd,
> t1 = rlnorm(Nind, log_T1_mu, log_T1_sd),
>
> log_M1 = log_M_mu,
> S_M1 = log_M_sd,
> m1 = rlnorm(Nind, log_M_mu, log_M_sd),
>
> log_R = log_R_mu,
> S_R = log_R_sd,
> r = rlnorm(Nind, log_R_mu, log_R_sd),
>
> t23 = t2 + t3,
> t2of3 = t2 / (t2 + t3),
> s1 = 0.12,
> sigma = 0.5
> )
> }
>
> Unfortunately, stan threw an error when I tried to feed it these inits (trying to run a single chain):
>
> "Error in config_argss(chains = chains, iter = iter, warmup = warmup, thin = thin, :
> initial value list mismatchs number of chains
> error in specifying arguments; sampling not done"
>
> I'm not sure what I'm doing wrong; I used the template in the stan-helpfile available in R to generate a list of initial
> values for each parameter in the parameter block. I vaguely remember that in the past, inits would have to be
> transformed to the unconstrained scale; is this still the case? Couldn't find this information in the latest docs. Also,
> I assume I don't have to provide inits for the parameters in the transformed parameters block, or do I?
>
> *Parameter constraints
> *I'm fairly certain that the parameter support and priors align well. What may be a problem though (I'm not sure), is
> that I highly constrained hyperpriors for R, M, and T, while allowing their children r, m, t vary freely (given some
> standard deviation).
>
> Thanks again for your advice,
>
> Luc

[Michael Betancourt]

The biggest issue here is that the chains aren’t sampling from

the same locations.  Are the PDFs on gDrive from multiple chains

or a single chain?

What would be helpful is to have the step size before and after

as well as the distribution of mass matrix elements.  

[Luc Coffeng]

Hi Bob, Michael,

To answer Bob's questions:

•  The PDFs on gDrive are from multiple chains, so yes, the chains aren't sampling from the same locations.
• The earlier version of stan did not explore the region that the current sampler is exploring. This doesn't bode well for the previous estimates, although they made much more sense (biologically) than the current ones. But I'm probably biased on that part, wanting to see those nice catterpillar plots again.
• Regarding the mass around the modes, I find it difficult to judge (or I am not getting the question right); in general the probability mass is distributed very irregularly in the lastest version of stan.
• Good to hear you confirm that the parameter support is probably not the problem (although I don't understand how you can tell from the treedepth, this is beyond my grasp, although I'd love to understand it).
• The quoted doc on initializing from R is exactly the one I'm using, so I'm not quite sure what I'm doing wrong. My init-generating function (see previous mail) basicly has the format that you describe ( init_fun <- function(x) { return(init_vals); } ). Just to check, I don't have to initialize the transformed parameters, right?

In response to Michaels' q's:

• I can get the step sizes later today when I'm home. Is mean stepsize useful information, or would you like something more detailed? I am sorry but I don't know how to extract the mass matrix from a stan object in R. Note on the side: it would be helpful if the R help-file on stan and/or the manual would list what kind of values are in an object resulting from stan(fit = ..., data = ..., ...) in R, and how this information can be pulled.

Thanks for thinking with me!

Luc

[Bob Carpenter]

On 2/18/14, 12:48 AM, Luc Coffeng wrote:
> Hi Bob, Michael,
>
> To answer Bob's questions:
>

> * The PDFs on gDrive are from multiple chains, so yes, the chains aren't sampling from the same locations.

So that does seem to indicate that the chains are getting stuck in local modes.

> * The earlier version of stan did not explore the region that the current sampler is exploring. This doesn't bode well

> for the previous estimates, although they made much more sense (biologically) than the current ones. But I'm
> probably biased on that part, wanting to see those nice catterpillar plots again.

This is probably where you want to focus debugging.

> * Regarding the mass around the modes, I find it difficult to judge (or I am not getting the question right); in

> general the probability mass is distributed very irregularly in the lastest version of stan.

Look at the points around the modes and see what their density is and how far out they
are. Hard to judge, but you can eyeball it. Very coarsely, you can look at a histogram
of lp__ on each chain. If one's very small compared to the others, you are probably stuck
in a bad local location.

Also look at step sizes. One way to get stuck is to adapt to a very
small step size and then continually hit the tree depth boundary.

> * Good to hear you confirm that the parameter support is probably not the problem (although I don't understand how you

> can tell from the treedepth, this is beyond my grasp, although I'd love to understand it).

When models have parameters that go out of support, it usually shows
up early and often with rejections. When NUTS wanders out of support,
it stops moving. Tree depth = log2(number of steps moved), so when
tree depth is large, you moved many steps. Now it could be that you moved
many steps because of very small step sizes.

> * The quoted doc on initializing from R is exactly the one I'm using, so I'm not quite sure what I'm doing wrong. My
> init-generating function (see previous mail) basicly has the format that you describe ( /init_fun <- function(x) {
> return(init_vals); }/ ). Just to check, I don't have to initialize the transformed parameters, right?

>
> In response to Michaels' q's:
>

> * I can get the step sizes later today when I'm home. Is mean stepsize useful information, or would you like something
> more detailed?

Unless you added variable stepsize, it'll be fixed after adaptation.
See the next response.

> I am sorry but I don't know how to extract the mass matrix from a stan object in R. Note on the side:
> it would be helpful if the R help-file on stan and/or the manual would list what kind of values are in an object
> resulting from stan(fit = ..., data = ..., ...) in R, and how this information can be pulled.

It's in the doc for stanfit-class. You want get_adaptation_info.

But we definitely need better, more tutorial doc. It's on my
short-term to-do list. We keep thinking we're going to write a book
about it, but getting some better doc out would be nice.

The problem is that it's grown organically, and the doc is in the usual
R-doc format with everything documented where it lives in the code, not
necessarily where it lives conceptually for users.

> Thanks for thinking with me!

We really want to help people fit hard models.

- Bob

[Luc Coffeng]

Ok here are some of the adaptation characteristics of the models from Aug 2013 and Feb 2014:

treedepth (not sure what the negative depth is):

                       -1      0      1     2    3       4      5    6      7     8    9

Feb 2014:    403  643  581 215  63     23    38   20   10     3    1
Aug 2013:        0      0      1     7  21 1545  483   12   10   20    1

stepsize:
Feb 2014:   0.0169881250045025
Aug 2013:   0.0630848314452992

mass matrix:
See attached txt file for the output of get_adaption_info() (rstan version 2.2.0) applied to my stan objects from Aug 2013 and Feb 2014. The Aug 2013 stuff seems to be far more than the diagonal extracted from the mass matrix (as is the case for Feb 2014).

posterior density:
The most recent version of stan samples from a relatively low density area of the posterior: log-probability mostly between -480 and -420 (chain 1) and between -360 and -320 (chain 2). All four chains of the model from Aug 2013, the log-probability of most samples was between -400 and -200.

fit to data:
To illustrate the implausibility of the results from the latest analyses, see attached pdf's comparing posterior (median and 95% credible intervals) and data for the analyses from Aug 2013 and Feb 2014. I realize now that maybe the prior for the period between parasite waves may be too uninformative (the wider this interval, the less pronounced the wave-pattern). Will try to fiddle with this.

initialization:
Still not sure what's going wrong here. I'm following the instructions in the R help (and those given by Bob himself). If anyone is interested to reproduce, use the script attached to my previous post.

another thing:
R sometimes freezes after completion of a chain. So far, this has only happened when I set up multiple chains (sequentially); this may be a coincidence but still peculiar.

Thanks again,

Luc

[Luc Coffeng]

Just tested the model again with simulated data, using double the number of data points per individual, lower measurement error, and super-constrained parameters, practically telling the model what the values of the hyperparameters are. Oddly enough, the model severely underestimates all log-normally distributed child-parameters (t[], m[], and r[]) relative to the estimated mean of their parent-hyperparameter (T1, M1, and R). Even more oddly enough, the parent-hyperparameter for the means themselves do not seem to be pulled towards the lower bound of their support, and the standard deviations S_T1, S_M1, and S_R are not pushed towards their upper bound to allow the crazy values of their child-parameters. This discrepancy is most obvious for parameter M1 and its children m[]. See excerpt from output below and attached pdf with density plots. Could it be that something is off with the log-normal density in rstan 2.2.0?

Parameter constraints (see attached R script to reproduce results):
log_T1 ~ uniform(log(6), log(8))              #  ~ uniform(1.79, 2.08)
  S_T1 ~ uniform(0.001,1)
         t ~ lognormal(log_T1, S_T1)

log_M1 ~ uniform(log(100), log(400))     # ~ uniform(4.61, 5.99)
  S_M1 ~ uniform(0.001, 3)
        m ~ lognormal(log_M1, S_M1)

log_R ~ uniform(log(12), log(16))           # ~ uniform(2.48, 2.77)
  S_R ~ uniform(0.001,1)
        r ~ lognormal(log_R, S_R)


Parameter estimates
Inference for Stan model: model_sim_data.
1 chains, each with iter=4000; warmup=2000; thin=1;
post-warmup draws per chain=2000, total post-warmup draws=2000.

            mean se_mean     sd      2.5%       25%       50%       75%     97.5% n_eff Rhat
log_T1      1.89    0.00   0.00      1.89      1.89      1.89      1.89      1.89     2 2.95
T1          6.62    0.01   0.01      6.60      6.60      6.62      6.63      6.63     2 2.95
t1[1]       2.39    0.01   0.03      2.34      2.37      2.38      2.41      2.45     5 1.81
t1[2]       1.51    0.02   0.03      1.47      1.49      1.51      1.54      1.59     4 1.40
t1[3]       5.68    0.09   0.23      4.96      5.73      5.77      5.79      5.81     6 1.20
t1[4]       1.62    0.03   0.05      1.54      1.56      1.63      1.66      1.70     2 2.81
t1[5]       1.57    0.04   0.06      1.48      1.51      1.57      1.62      1.65     2 2.69
t1[6]       1.26    0.02   0.03      1.21      1.23      1.26      1.29      1.31     3 2.91
t1[7]       1.01    0.01   0.02      0.98      1.00      1.01      1.02      1.04     4 1.62
.
.
.
S_T1        0.91    0.00   0.00      0.90      0.91      0.91      0.91      0.91     5 1.00
t2          0.36    0.00   0.01      0.35      0.35      0.36      0.36      0.36     3 2.63
t3          1.11    0.01   0.02      1.08      1.09      1.11      1.13      1.14     3 2.68
s1          0.28    0.01   0.02      0.23      0.28      0.29      0.30      0.30     4 1.55
log_M1      4.85    0.00   0.00      4.84      4.85      4.85      4.85      4.86     4 1.64
M1        127.84    0.31   0.59    126.86    127.35    127.86    128.28    128.95     4 1.64
m1[1]       0.37    0.00   0.01      0.36      0.37      0.37      0.38      0.38     7 1.10
m1[2]       1.80    0.02   0.04      1.75      1.76      1.80      1.84      1.86     3 2.71
m1[3]       1.38    0.04   0.06      1.29      1.34      1.38      1.41      1.53     3 1.76
m1[4]       2.85    0.03   0.06      2.71      2.81      2.88      2.90      2.93     4 1.56
m1[5]       5.20    0.04   0.10      5.02      5.14      5.18      5.27      5.38     6 1.15
m1[6]       1.54    0.01   0.02      1.49      1.52      1.54      1.55      1.57     7 1.05
m1[7]       2.12    0.09   0.16      1.91      1.97      2.11      2.27      2.38     3 3.23
.
.
.
S_M1        2.78    0.00   0.01      2.77      2.78      2.78      2.79      2.79     3 2.30
log_R       2.56    0.00   0.00      2.56      2.56      2.56      2.56      2.56     3 1.87
R          12.98    0.01   0.01     12.96     12.97     12.97     12.98     12.99     3 1.87
r[1]        2.18    0.01   0.02      2.14      2.17      2.18      2.20      2.22    12 1.03
r[2]        1.22    0.02   0.03      1.17      1.19      1.21      1.26      1.27     3 2.91
r[3]        0.63    0.01   0.01      0.61      0.62      0.63      0.64      0.65     5 1.13
r[4]        0.96    0.01   0.01      0.93      0.95      0.96      0.97      0.98     4 1.69
r[5]        3.35    0.05   0.09      3.20      3.29      3.36      3.42      3.50     3 1.98
r[6]        1.38    0.02   0.04      1.30      1.35      1.37      1.41      1.43     3 2.16
r[7]        1.26    0.01   0.02      1.22      1.25      1.26      1.28      1.29     4 1.34
.
.
.
S_R         0.99    0.00   0.00      0.99      0.99      0.99      0.99      0.99     4 1.60
sigma       3.00    0.00   0.00      3.00      3.00      3.00      3.00      3.00     4 1.05
lp__   -13588.22  162.71 269.36 -14176.20 -13774.29 -13519.97 -13371.56 -13209.67     3 2.23

[Luc Coffeng]

Just woke up realizing that there was a mistake in my model regarding the integrating out of censored observations. Before, I was integrating out the logarithm of the detection limit using the lognormal_cdf_log function, whereas I should have used the value of the detection limit as is. After correcting this however, similarly odd results still turn up (as described in my previous post) In case you want to test the corrected model, the script is attached. I tested the lognormal_cdf_log function with a simple censored regression example, and everything looked fine there (attached script).

[Michael Betancourt]

You have to compare apples to apples — diag to diag or dense to dense.  The

numbers quoted here don’t say much.

[Luc Coffeng]

Dear Michael,

Comparing appels and oranges (or pears as we say in Dutch) is no good, no :). I'm not sure though what apples & apples you would like to see, or how to get them, or what the oranges are in this particular case. I just checked and confirmed that I've extracted adaptation info from models that had the option non_diag_mass = TRUE.

Thanks,

Luc

[Michael Betancourt]

You gave a dense mass matrix in one case and a diagonal one in the other.  Pick diag or dense for both before and after.

[Bob Carpenter]

Also, the older versions of Stan are still linked from
the mc-stan.org page for CmdStan.

Maybe you could verify that the model’s still the same
and still getting the answers you think you should get in
an older version of Stan.

- Bob

[Luc Coffeng]

Ok, I'd be happy to provide. Could you please tell how to get the dense mass matrix from a stan-fit object that was generated with option non_diag_mass = TRUE? Apparently, get_adaptation_info() doesn't do the job.

[Luc Coffeng]

I'd love to Bob. Unfortunately, I currently really don't have the time to get familiar with CmdStan, so unless rstan v1.3.0 is available, it will have to wait for a conseridable while. This model is really a pet-project (that has totally gotten out of hand btw - tends to happen with fun things).

[Ben Goodrich]

You can get old rstan releases, including v1.3.0, from

https://github.com/stan-dev/rstan/releases

Ben

[Luc Coffeng]

Thanks for the pointer Ben!

I managed to install and run the model based on simulated data (exactly the same dataset!) in both rstan 1.3.0 and 2.2.0. Turns out I'm not going crazy after all, phew! Rstan 1.3.0 neatly fit the model to the data and comes up with good parameter estimates and reasonable effective sample sizes. Rstan 2.2.0 however gets stuck somewhere in the posterior, while the only difference in the model code is the following (N.B. the simulated data are the same):

rstan 1.3.0 (vectorization and logarithmic version of log-normal cdf not yet supported):
    for (i in 1:N_cens1) {
        lp__ <- lp__ + log(lognormal_cdf(dlimit_cens1[i], log_RESP_CENS1[i], sigma));
    }
    for (j in 1:N_cens2) {
        lp__ <- lp__ + log(lognormal_cdf(dlimit_cens2[j], log_RESP_CENS2[j], sigma));
    }


rstan 2.2.0:
    increment_log_prob(lognormal_cdf_log(dlimit_cens1, log_RESP_CENS1, sigma));
    increment_log_prob(lognormal_cdf_log(dlimit_cens2, log_RESP_CENS2, sigma));


Here is the adaptation info (diagonal matrix for both versions of rstan):
Rstan v1.3.0: "# (mcmc::nuts_diag) adaptation finished\n# step size=0.158287\n# parameter step size multipliers:\n# 0.863175,0.509413,0.0207546,0.032371,0.0348721,0.030029,0.018976,0.0342302,0.0273426,0.0262984,0.021153,0.03516,0.0263005,0.0271261,0.0324769,0.0335648,0.0216606,0.0297665,0.0244665,0.0230258,0.0281096,0.0271063,0.0207699,0.0249988,0.0219923,0.0278222,0.0231806,3.43759,1.14202,0.698956,0.877757,0.849057,0.774022,0.692552,0.820241,0.834916,0.766557,0.722882,0.994702,0.756277,0.856919,0.770622,0.83905,0.765962,0.761526,0.758154,0.699609,0.830871,0.82294,0.701233,0.740589,0.705349,0.771324,0.706271,4.19726,1.32636,0.610505,0.783775,0.732568,0.683288,0.568206,0.890825,0.761874,0.597916,0.511704,0.848955,0.638742,0.759727,0.571279,0.803412,0.596836,0.665227,0.620462,0.641927,0.706293,0.575611,0.73828,0.517816,0.567233,0.619472,0.647546,0.25342,3.79881,3.20188,0.146344\n"

Rstan v2.2.0: "# Adaptation terminated\n# Step size = 0.0333261\n# Diagonal elements of inverse mass matrix:\n# 9.95611e-005, 9.78388e-005, 0.0004322, 4.77209e-005, 5.51738e-005, 0.00012751, 3.61822e-005, 3.69474e-005, 0.00014984, 7.07592e-005, 7.4052e-005, 4.02197e-005, 0.000197432, 0.000146699, 0.00019049, 2.34273e-005, 2.75238e-005, 0.000137562, 4.58576e-005, 5.7502e-005, 0.000274897, 0.000130277, 1.87729e-005, 5.92964e-005, 0.000146668, 0.000140817, 3.26834e-005, 6.49846e-005, 0.000155731, 5.47277e-005, 0.000166372, 3.63364e-005, 0.000137372, 5.00938e-005, 4.76363e-005, 0.000197817, 0.00020613, 9.2058e-005, 8.02352e-005, 2.73686e-005, 6.59215e-005, 4.17211e-005, 0.00044265, 3.68804e-005, 9.29263e-005, 0.00022285, 0.000430258, 6.62486e-005, 0.000148667, 6.33311e-005, 9.21591e-005, 7.26858e-005, 0.000102454, 7.18484e-005, 0.000102984, 0.000115399, 4.07442e-005, 5.57282e-005, 3.69454e-005, 9.09086e-005, 3.34457e-005, 2.45772e-005, 0.00016978, 2.69151e-005, 7.47625e-005, 3.4565e-005, 0.000323954, 7.71429e-005, 2.57016e-005, 7.18276e-005, 4.19812e-005, 0.00010343, 5.77683e-005, 0.000223856, 0.000275108, 0.000156034, 9.97131e-005, 8.9452e-005, 2.44022e-005, 8.95875e-005, 0.000102667, 6.17894e-005, 0.000152968, 0.000345607, 3.58389e-005\n"


The "Adaptation terminated" thrown by rstan v2.2.0 doesn't good... Any idea why this would happen? Could it have to do with the lognormal_cdf_log function? I tried replacing it with log(lognormal_cdf), but this didn't change anything.


If anyone is interested, here are links to the RData-files and scripts (google drive):
Rstan v1.3.0: R script, RData file containing simulated data and stan-fit object
Rstan v2.2.0: R script, RData file containing simulated data and stan-fit object

Best,

Luc

[Luc Coffeng]

Interestingly, rstan v2.2.0 does seem to manage when I give it the stan-fit S4 object from v1.3.0 to work with (rather than compiling the model in v2.2.0). In this case, adaptation finishes successfully, and the posterior is sampled similarly well as by rstan v1.3.0. So either something is going wrong when compiling the model in rstan v2.2.0, and/or the adaptation breaks for some reason.

[Jiqiang Guo]

On Tue, Feb 18, 2014 at 8:47 PM, Luc Coffeng <lucco...@gmail.com> wrote:

Interestingly, rstan v2.2.0 does seem to manage when I give it the stan-fit S4 object from v1.3.0 to work with (rather than compiling the model in v2.2.0).

Here I am afraid that you are using the compiled model from 1.3 in 2.2.  So in fact, v1.3 is used.  The stanfit object includes the binary shared object (dynamic link library) created from compiling the c++ code for the model.  When a saved stanfit object is reused, the shared object would be loaded and get executed for sampling.  So in v2.2, it would recompile the c++ code when using a stanfit object created in 1.3.  (I am not sure it's a bug or not, but this scenario is rare.)

Jiqiang 

[Jiqiang Guo]

On Tue, Feb 18, 2014 at 9:12 PM, Jiqiang Guo <guo...@gmail.com> wrote:

So in v2.2, it would recompile

I meant to say it would NOT ...

[Luc Coffeng]

Hi Jiqiang, that's exactly what I meant to do :). The model compiled by v1.3 seems to work well (but maybe I'm biased here, so I'm eager to hear what you stanfolks think), while the version compiled by v2.2.0 systematically terminates adaptation early and fails to sample efficiently.

On Tuesday, 18 February 2014 18:12:59 UTC-8, Jiqiang Guo wrote:

On Tue, Feb 18, 2014 at 9:12 PM, Jiqiang Guo <guo...@gmail.com> wrote:

So in v2.2, it would recompile

I meant to say it would NOT ...

[Michael Betancourt]

The problem is evident in the metric elements, which fall to very small values in 2.2.0.  This

happens when your chain gets stuck during adaptation such that the covariance becomes

very small.  In 1.3.0 the covariance was normalized so that the global scale canceled out, but

that really masked the issue (not the large number of iterations that failed after a single leapfrog

step — this is indicative of some region of parameter space behaving very poorly).

My best guess is that your model has funky tails that are causing problems, but it would take

some serious investigation to narrow down the exact source of the problem.  That said, the

fact that the log of CDFs can very quickly hit infinity is something of which to be aware.

[Luc Coffeng]

Thanks Michael!

Took me a while to understand what you said, but it makes sense to me now. I assume that iterations that fail after a single leapfrog step are reflected by treedepth = 0, or is it -1?

The funky tail business can be explained by the cyclical pattern in the regression line. You could theoretically shift the regression line by one period and scale it by the reproduction factor so that it still fits the data. I tried to anchor the model to prevent this (by making sure that the log CDF would be close to -Inf for unlikely values), but apparently this now backfires with the mass matrix no longer being normalized. I guess another way to do it would be to truncate some of the distributions; there is reasonable prior knowledge of what the parameter values should(n't) be and maybe I just wasn't telling it very well to the model :).

Luc

[Michael Betancourt]

Took me a while to understand what you said, but it makes sense to me now. I assume that iterations that fail after a single leapfrog step are reflected by treedepth = 0, or is it -1?

It’s the tree_depth =  -1 iterations that are indicative of an issue — although be careful because we changed

the definition of tree depth in 2.2.0 to be consistent with the typical computer science definition, in which

case it will be tree_depth = 0.

The funky tail business can be explained by the cyclical pattern in the regression line. You could theoretically shift the regression line by one period and scale it by the reproduction factor so that it still fits the data. I tried to anchor the model to prevent this (by making sure that the log CDF would be close to -Inf for unlikely values), but apparently this now backfires with the mass matrix no longer being normalized. I guess another way to do it would be to truncate some of the distributions; there is reasonable prior knowledge of what the parameter values should(n't) be and maybe I just wasn't telling it very well to the model :).

Here’s my guess as to what is happening: the anchor introduces a region of very high curvature when the log CDF rapidly jumps to -infinty.  In order to

explore that region the integrator requires a very small step size or it will diverge and any samples rejected.  If you naively set the target acceptance

rate the step size will be too high and this region will essentially be ignored.  This will induce some bias in your samples, although the size of the bias

will depend on the specifics of your model.  This is one of the reasons why we increased the default target acceptance probability to 0.8.

Now this also manifests in adaptation.  When you enter this region the samples are continuously rejected until the step size becomes small enough

that the integrator is reasonably accurate.  These samples, however, bias the covariance estimates low which leads to the very small metric elements.

So try running with a larger target acceptance probability and see if that helps at all.

[Luc Coffeng]

Dear Bob,

Thanks for your input! Apologies for the lack of model introduction on my part.

I'm working on a hierarchical model for blood parasite levels in experimentally infected subjects, assuming a log-normal distribution for blood parasite levels. The blood parasites appear in cycles of generations (called waves in

the model, parameter T representing the timing of the center of the first wave), each cycle being modeled by two scaled normal cdf's: one for the build-up of parasites, another for the disappearance of parasites. Both CDF's are scaled by parameter M, which represents the amplitude of the first parasite wave. Due to parasite reproduction (this happens when parasites are not visible in the blood), consecutive parasite waves have increasing amplitude (multiplication factor represented by parameter R), resulting in a typical pattern (see attached pdf). Timing (T), amplitude (M), and multiplication rate (R) are all assumed to vary between persons. To anchor the model, I'm explicitly modeling censored observations, and I assume that all observations up to the first observation above the detection limit are samples of first generation of blood parasites.

Mass matrix

The full mass matrix doesn't seem to be the problem, apart from taking longer to complete. With a diagonal mass matrix, N_eff is still low (<20 for 1000 samples in 60 seconds). Applying brute force with 20000 samples (and thin = 10, warmup = 2000) doesn't help.

Acceptance and treedepth

With diagonal mass matrix, 2000 warmup, 4000 samples, init = 0, and 0.65 target acceptance rate:

Tree depth:

-1    0     1     2      3    4    5    6    7     8   9

403 643 581 215  63  23  38  20  10   3   1

Average acceptance rate was ~65%, as targeted.

Multimodal posterior
Indeed, for some parameters the posterior seems to be multimodal. However, the multimodality of the posterior results in this case seems to be due to chains getting stuck in different areas of the posterior. Somehow, I need to get the sampler to behave more flexibly, and my impression is that it was able to explore the posterior much more efficiently before; I used to get nice caterpillar plots and chains used to converge to the same area of the posterior. Now, the traceplots are very jagged and chains end up in different places :(.

Initializing parameter values

So far, I tried initializing with random inits or init = 0 (before rstan 2.x, I let rstan initialize itself). Today, I tried feeding the model a list of inits (and lists of lists of inits for multiple chains), using the actual parameter values that I used to simulate the data (also added this to the attached r script):

  gen.init <- function() {
    list(
        log_T1  = log_T1_mu,
        S_T1    = log_T1_sd,
        t1      = rlnorm(Nind, log_T1_mu, log_T1_sd),
       
        log_M1  = log_M_mu,
        S_M1    = log_M_sd,
        m1      = rlnorm(Nind, log_M_mu, log_M_sd),
       
        log_R   = log_R_mu,
        S_R     = log_R_sd,
        r       = rlnorm(Nind, log_R_mu, log_R_sd),
       
        t23     = t2 + t3,
        t2of3   = t2 / (t2 + t3),
        s1      = 0.12,
        sigma   = 0.5
    )
  }

Unfortunately, stan threw an error when I tried to feed it these inits (trying to run a single chain):

"Error in config_argss(chains = chains, iter = iter, warmup = warmup, thin = thin,  :
  initial value list mismatchs number of chains
error in specifying arguments; sampling not done"

I'm not sure what I'm doing wrong; I used the template in the stan-helpfile available in R to generate a list of initial values for each parameter in the parameter block. I vaguely remember that in the past, inits would have to be transformed to the unconstrained scale; is this still the case? Couldn't find this information in the latest docs. Also, I assume I don't have to provide inits for the parameters in the transformed parameters block, or do I?

Parameter constraints

I'm fairly certain that the parameter support and priors align well. What may be a problem though (I'm not sure), is that I highly constrained hyperpriors for R, M, and T, while allowing their children r, m, t vary freely (given some standard deviation).

Thanks again for your advice,

Luc

On Monday, 17 February 2014 01:50:44 UTC-8, Bob Carpenter wrote: