Poor mixing of chains in model with randomly sampled errors - pointers?

[pcoll...@berkeley.edu]

I'm running models in STAN to try to estimate parameters in a bi-exponential function describing pathogen die-off:

ln(C(t)) = ln(C(0)) + ln(a*e^(-k1*t)+(1-a)*e^(-k2*t))

Where C(t) is a function of pathogen concentration over time t, C(0) is the initial concentration, a is a proportion of pathogens undergoing rapid decay specified by rate parameter k1, and 1-a is the complementary proportion undergoing slower die-off. In the model, a is estimated based on the location of the zero of the third derivative of the function, which is constrained by timepoints at which measurements were taken (can't be earlier than point midway between first two observations, can't be later than second to last data point).

Some experimental data for which I'm attempting to estimate parameters include measures of uncertainty and others don't. To attempt to integrate both data types into the same framework, I've created a model (attached) that imputes errors at each observed concentration timepoint according to a gamma distribution fit to the reported errors from other experiments at the same time as attempting to estimate the decay parameters. The model seems to run reasonably well much of the time, though chains will occasionally get stuck in particular regions of the sample space (see goodrun_trace.pdf, goodrun_dens.pdf, goodrun_data.pdf). Occasionally I run into a set of experimental data for which the effective sample size is very low, even for extremely high (>1,000,000) numbers of iterations (see badrun_trace.pdf, badrun_dens.pdf, badrun_data.pdf). In these cases, I've observed that the mixing of chains is quite poor, and chains spend more time stuck in narrow ranges of the probability space.

Are there any adjustments to model specification or data transformation that I should consider in order to eke out better performance? I have a relatively simplistic understanding of the workings of STAN at this point, so if there is material I should familiarize myself with prior to soliciting specific input, I would welcome direction towards that as well. 

[Andrew Gelman]

Hi, I don’t see your Stan model in your attachments.  But my quick questions are:  (a) do you have informative priors for the tricky-to-estimate parameters, and (b) if you do have informative priors, are you using good initial values for your Stan runs?  I’ve often found that a combination of informative priors and good inits can solve this sort of problem of chains getting stuck in a bad neighborhood.  The priors can be weakly informative—you don’t have to “cheat” and stick your answer into your prior; the key is just to rule out various pathological zones.

A

--
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.
<badrun_data.pdf><badrun_data.rdump><badrun_dens.pdf><badrun_trace.pdf><goodrun_data.pdf><goodrun_data.rdump><goodrun_dens.pdf><goodrun_trace.pdf><Biphasic_Val_Stan_Imputes.rda>

[pcoll...@berkeley.edu]

Hi Dr. Gelman,

   Thanks for your reply. The STAN model is the last attachment 'Biphasic_Val_Stan_Imputes.stan'. Priors are uninformative, though the parameters are constrained within the declaration block. I could probably supply better initial parameter values. I'm not sure of the best way to rule out pathological zones, since the parameter values that seem to result in poor mixing are not particularly unreasonable within the context of the model. I suppose I could put more aggressive constraints on Chgpt, though I'm a little wary of doing so, as I would like the model formulation to be as generalizable as possible. I'll play around with supplying informative initial values first and see if that helps.

  Thanks,

   Philip

[Andrew Gelman]

Hi Philip.  There was no .stan attachment in that email of yours.  I think it could make sense to put a bit of prior in your model.  Just a bit of curvature keeping you away from iffy zones could make a difference.  I understand your desire to make the formulation as general as possible, but sometimes you just have to think of different regions of parameter space as representing different models!

A

[pcoll...@berkeley.edu]

Hi Dr. Gelman,

  Strange... I can see the attachment in the original post on google groups. At any rate, I was able to improve things a bit by creating dynamic constraints for the parameters of the model that vary according to properties of each dataset I want to run the model for. I also was able to get a bit better performance by making the prior distributions for variance and mean parameters a bit less flat. Examining trace plots of problematic runs, I noticed that periods where sampling efficiency of a chain was very poor seemed to correspond to times when a couple of key variance parameters were simultaneously very close to 0. Placing constraints on these parameters so that they could not take on values below 1E-6 seemed to prevent the worst of the inefficient mixing behavior. Model code posted below:

// Biphasic Decay Model With Imputed Concentration Variability 

data{

   int<lower=0> Nobs; 

   vector[Nobs] ln_Conc;

   vector[Nobs] Time;

   real K1upp;

   real K2upp;

   real diffhi;

   real difflo;

   real K1lo;

   real K2lo;

   real alphaupp;

   real alphalo;

   real ln_Conc0;

   real<lower=0> Errs2rate;

   real<lower=0> Errs2shape;

   }

parameters {

   real<lower=0> var_ln_Conc[Nobs];

   real<lower=0,upper=25> mu_ln_Conc0;

   real<lower=K1lo,upper=K1upp> K1;

   real<lower=difflo,upper=diffhi> diff;

   real<lower=difflo,upper=diffhi> mu_diff;

   real<lower=K1lo,upper=K1upp> mu_K1;

   real<lower=1E-6> var_diff;

   real<lower=1E-6> var_K1;

   real<lower=1E-6> var_ln_Conc0;

   real<lower=alphalo, upper=alphaupp> alpha;

   }

transformed parameters {

   real<lower=0> sigma_ln_Conc[Nobs];

   real<lower=K2lo,upper=K2upp> K2;

   real<lower=0> sigma_diff;

   real<lower=0> sigma_K1;

   real<lower=0> sigma_ln_Conc0;

   K2<-K1-exp(diff);

   sigma_diff<-sqrt(var_diff);

   sigma_K1<-sqrt(var_K1);

   for (i in 1:Nobs){

      sigma_ln_Conc[i]<-sqrt(var_ln_Conc[i]);

      }

   sigma_ln_Conc0<-sqrt(var_ln_Conc0);

   }

model{

   mu_diff~normal(0,8);

   mu_K1~normal(0,100);

   var_diff~cauchy(0,1);

   var_K1~cauchy(0,15);

   diff~normal(mu_diff, sigma_diff);

   K1~normal(mu_K1,sigma_K1);

   var_ln_Conc0~gamma(Errs2shape,Errs2rate);

   mu_ln_Conc0~normal(ln_Conc0, sigma_ln_Conc0);

   var_ln_Conc~gamma(Errs2shape,Errs2rate);

   ln_Conc~normal((mu_ln_Conc0+log(alpha*exp(-Time*K1)+(1-alpha)*exp(-Time*K2))), sigma_ln_Conc); 

   }

Thanks for the help and suggestions!

   Phil 

[Bob Carpenter]

> On May 16, 2016, at 4:10 PM, pcoll...@berkeley.edu wrote:
>
> Hi Dr. Gelman,
>

> Strange... I can see the attachment in the original post on google groups.

It didn't come through on my mailer, either. I had 9 attachments.

You really need to consider the non-centered parameterizations for
efficiency if there are small data sizes int he groups.

And you don't want super-diffuse gamma priors (I didn't try to
figure out how to inspect your data). And you can just mode scales
directly --- Stan doesn't use any conjugacy information.

Is this right:

> mu_ln_Conc0~normal(ln_Conc0, sigma_ln_Conc0);

sigma and ln_Conc0 are vectors, but mu_ln_Conc0 isn't, so you're
modeling mu_ln_Conc0 Nobs times.

- Bob

> At any rate, I was able to improve things a bit by creating dynamic constraints for the parameters of the model that vary according to properties of each dataset I want to run the model for. I also was able to get a bit better performance by making the prior distributions for variance and mean parameters a bit less flat. Examining trace plots of problematic runs, I noticed that periods where sampling efficiency of a chain was very poor seemed to correspond to times when a couple of key variance parameters were simultaneously very close to 0. Placing constraints on these parameters so that they could not take on values below 1E-6 seemed to prevent the worst of the inefficient mixing behavior.

If lower=1e-6 makes a difference, then you have a problem with estimates
near the boundary and adding that lower bound is just going to mask that.

What do the posteriors look like?

- Bob

[pcoll...@berkeley.edu]

Hi Bob,

   Thanks for the reply. Didn't see it until just now. 

On Tuesday, May 17, 2016 at 12:03:35 PM UTC-7, Bob Carpenter wrote:

> On May 16, 2016, at 4:10 PM, pcoll...@berkeley.edu wrote:
>
> Hi Dr. Gelman,
>
>   Strange... I can see the attachment in the original post on google groups.

It didn't come through on my mailer, either.  I had 9 attachments.

The 9th attachment is the STAN model. What showed up in the email? Sorry for all the confusion around this, I'm still not sure what's going on...

You really need to consider the non-centered parameterizations for
efficiency if there are small data sizes int he groups.

Thank you for suggesting this. I'm struggling somewhat to understand how centered vs. non-centered parameterizations work. Basically the sampling of the estimated parameters is independent from the likelihood estimation in the NCP? In the case of the above model, if for example I wanted to have a non-centered parameterization for K1, would I declare it like so?

Centered:

model{

K1~normal(mu_K1, sigma_K1)

mu_K1~normal(0,100)

sigma_K1~cauchy(0,15)

}

Non-centered:

transformed parameters{

K1<-mu_K1+sigma_K1*z_K1

}

model{

mu_K1~normal(0,100)

sigma_K1~cauchy(0,15)

z_K1~normal(0,1)

}

 

And you don't want super-diffuse gamma priors (I didn't try to
figure out how to inspect your data).  

The gamma distribution estimated from the known variances is not very diffuse, I think. Parameters are .4 and .54, tails don't get out beyond ~25. 

 

And you can just mode scales
directly --- Stan doesn't use any conjugacy information.

Meaning I should go ahead and draw samples of the sigma parameters directly, rather than worrying about doing it through the variances?

 

Is this right:

>    mu_ln_Conc0~normal(ln_Conc0, sigma_ln_Conc0);

sigma and ln_Conc0 are vectors, but mu_ln_Conc0 isn't, so you're
modeling mu_ln_Conc0 Nobs times.

- Bob

ln_Conc0 and sigma_ln_Conc0 are single values in this model, whereas ln_Conc and sigma_ln_Conc are vectors. mu_ln_Conc0 is only sampled once per iteration.

 

> At any rate, I was able to improve things a bit by creating dynamic constraints for the parameters of the model that vary according to properties of each dataset I want to run the model for. I also was able to get a bit better performance by making the prior distributions for variance and mean parameters a bit less flat. Examining trace plots of problematic runs, I noticed that periods where sampling efficiency of a chain was very poor seemed to correspond to times when a couple of key variance parameters were simultaneously very close to 0. Placing constraints on these parameters so that they could not take on values below 1E-6 seemed to prevent the worst of the inefficient mixing behavior.

If lower=1e-6 makes a difference, then you have a problem with estimates
near the boundary and adding that lower bound is just going to mask that.

What do the posteriors look like?

- Bob

Ah, I was afraid that might be a bit cheaty. I'll rerun the models overnight with all parameters saved and post some posterior distributions (I assume you were interested in the scale parameters?) tomorrow.

Thanks again,

   Philip

[Bob Carpenter]

> On May 31, 2016, at 10:16 PM, pcoll...@berkeley.edu wrote:
>
> Hi Bob,
>
> Thanks for the reply. Didn't see it until just now.
>
> On Tuesday, May 17, 2016 at 12:03:35 PM UTC-7, Bob Carpenter wrote:
>
> > On May 16, 2016, at 4:10 PM, pcoll...@berkeley.edu wrote:
> >
> > Hi Dr. Gelman,
> >
> > Strange... I can see the attachment in the original post on google groups.
>
> It didn't come through on my mailer, either. I had 9 attachments.
>
> The 9th attachment is the STAN model. What showed up in the email? Sorry for all the confusion around this, I'm still not sure what's going on...

I really don't want to type the name of 9 files out.
There were 6 .pdf, 2 .rdump, and 1 .rda file. I was
looking for something ending in .stan.

But really 9 attachments are more than we realistically
have time to inspect.

> You really need to consider the non-centered parameterizations for
> efficiency if there are small data sizes int he groups.
>
> Thank you for suggesting this. I'm struggling somewhat to understand how centered vs. non-centered parameterizations work.

See: https://arxiv.org/pdf/1312.0906.pdf

> Basically the sampling of the estimated parameters is more independent from the likelihood estimation? In the case of the above model, if for example I wanted to have a non-centered parameterization for K1, would I declare it like so?

>
> Centered:
> model{
> K1~normal(mu_K1, sigma_K1)
> mu_K1~normal(0,100)
> sigma_K1~cauchy(0,15)
> }
>
> Non-centered:
> transformed parameters{
> K1<-mu_K1+sigma_K1*z_K1
> }
>
> model{
> mu_K1~normal(0,100)
> sigma_K1~cauchy(0,15)
> z_K1~normal(0,1)
> }

Exactly. I'd recommend spaces and indentation for readability.

Do you really expect mu_K1 to take absolute values > 100
and sigma_K1 to take arbitrarily large values? If not, you can
tighten these priors, too.

> And you don't want super-diffuse gamma priors (I didn't try to
> figure out how to inspect your data).
>
> The gamma distribution estimated from the known variances is not very diffuse, I think. Parameters are .4 and .54, tails don't get out beyond ~25.
>
>
>
>
>

> And you can just mode scales
> directly --- Stan doesn't use any conjugacy information.
>
>
> Meaning I should go ahead and draw samples of the sigma parameters directly, rather than worrying about doing it through the variances?

We usually just model sigma directly (with a prior).
The reason it's easier than variance or precision is that
it's on the natural scale of the data, so easier to formulate
reasonable priors for.

Stan samples everything jointly, though --- so you never sample
sigma directly all by itself.

> > At any rate, I was able to improve things a bit by creating dynamic constraints for the parameters of the model that vary according to properties of each dataset I want to run the model for. I also was able to get a bit better performance by making the prior distributions for variance and mean parameters a bit less flat. Examining trace plots of problematic runs, I noticed that periods where sampling efficiency of a chain was very poor seemed to correspond to times when a couple of key variance parameters were simultaneously very close to 0. Placing constraints on these parameters so that they could not take on values below 1E-6 seemed to prevent the worst of the inefficient mixing behavior.
>
> If lower=1e-6 makes a difference, then you have a problem with estimates
> near the boundary and adding that lower bound is just going to mask that.
>
> What do the posteriors look like?
>
> - Bob
>
> Ah, I was afraid that might be a bit cheaty. I'll rerun the models overnight with all parameters saved and post some posterior distributions (I assume you were interested in the scale parameters?) tomorrow.

The problem is that if the properly formulated models wants to draw values
near 1E-6, then you don't want to put an artificial wall there.

- Bob

[pcoll...@berkeley.edu]

Thanks Bob,

   Apologies about the excessive attachments. The .rda file is the STAN file, so I guess I will have to change the file type beyond editing the file name in the future. As it turns out, the typical posterior distributions for scale parameters in the model don't have a lot of density near 0. All the same, I'll proceed without those lower bounds and see if I can get better performance using non-centered parameterizations.

  Philip