general questions regarding divergent transitions after warmup

[Jaslene Lin]

Dear Stan developers:

I am relatively new to Stan and I am currently running a hierarchical model using Rstan and my model complied however it returns me with the error message that there are around 3 (sometimes upto 5) divergent transitions after warmup. I have tried increasing (adapt_delta=0.9999999999,max_treedepth=18,stepsize=0.001) however, it still behaves this way.

Firstly, if there is any divergent transitions then my posterior samples cannot be trusted right? So my top priority is to work out what causes the divergent transitions before using the posteriors distribution samples to make any inferences right? I have tried give a more informative priors but it does not seem to work. Can someone please provide me with some ideas? Any thoughts would be greatly appreciated.

Kindest regards,

Jaslene

[Jaslene Lin]

I want to estimate 7 observed species categories along with 5 unobserved species categories as a bivariate model transferring the data into c(midpoints, log(range)) as data come in as intervals. some data only have midpoints(centre estimates recorded and thus their corresponding log(range)) are coded as parameters to be estimated in the model. There are potentially 32 combinations of models (2^5, depending on which UNOBSERVED species categories are in the model and rather than writing 32 separate files for each combination of model, I would write a single model using indicator such (1,1,1,1,1) denoting which UNOBSERVED species categories are in that particular model. In the model, I record DIC and the aim is to compare each combination of model and choose the one with the smallest DIC to be able to say that particular unobserved species are indeed in the model/ be able to estimated by the model given the data.

 

Also there is intrinsic hierarchical structure within the species categories. Such as, insects+arthropods=terrestrial etc.  

I have attached my Rstan code in the file. With the current construction, the model ran for (1,1,1,1,1), in other words, all other unobserved species are in the model but when I ran for another extreme case (0,0,0,0,0) the model ran extremely slowly. Also when the model did run, there is always error saying that there are divergent transitions and I do not want to proceed comparing DIC without fixing this problem.

One odd thing is that when estimating missing data, (midpoints observed, log(range)=NA), the model returns all observed midpoints are reported with Rhat=NaN, I think that is because they are data and the model does not need to estimate them rather using them to estimate their corresponding log(range), which makes sense. However, for some of the observations, the model seems to report n_eff=very low number with Rhat roughly around 1. Does this indicate that the model tries to estimate the data I imported?

Can someone please provide me with some insights into my problems and any suggestions would be greatly appreciated.

Thanks in advance,

Jaslene

[Jaslene Lin]

I also set constraints for midpoints parameters to ensure that none of them are negative and also constraints for log(range) which is a function of midpoints to ensure that lower bounds of species intervals are non-negative.  

[Bob Carpenter]

I'm afraid I don't have time to dig through long
models. Is there some way to use loops to cut down on
the code size? Or to summarize what you're trying to do?
Often, problems arise from hierarchical models with centered
parameterizations and the non-centered approach will fit
much better.

The divergent transitions are just places where Stan has
numerical problems simulating the Hamiltonian. The low
n_eff is more worrying.

I didn't follow the missing data question. Stan doesn't
allow NA data from R, so it all has to be modeled. If
it's a parameter (being estimated), then it's not data.
R-hat will be NaN in cases where all of the draws have the
exact same value.

- 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.
> <combo model real data.stan>

[Jaslene Lin]

Hi Bob:

Thank you so much for your reply and I will look into using non-centred parameterisation for my model.

For the missing data problems, I have 16 observations where their log(range) are missing and I code them as parameters to be estimated in the model and the following the are the results. Where [5,1] indicating the 5th missing observations and [5,1] is the observed midpoints and [5,2] is the log(range) estimated by the model.

y1_complete[5,1]          0.30    0.00     0.00       0.30       0.30       0.30       0.30       0.30     2 0.99

y1_complete[5,2]          0.97    0.12     3.37      -4.96      -0.89       0.87       2.65       8.11   733 1.00

y1_complete[6,1]          8.00    0.00     0.00       8.00       8.00       8.00       8.00       8.00   800  NaN

y1_complete[6,2]          1.99    0.03     0.81       0.55       1.46       1.93       2.44       3.69   747 1.00

y1_complete[7,1]          2.00    0.00     0.00       2.00       2.00       2.00       2.00       2.00   800  NaN

y1_complete[7,2]         -0.28    0.03     0.83      -2.04      -0.81      -0.23       0.29       1.12   760 1.00

y1_complete[8,1]          4.00    0.00     0.00       4.00       4.00       4.00       4.00       4.00   800  NaN

y1_complete[8,2]          0.47    0.03     0.75      -1.07       0.01       0.51       0.93       1.84   800 1.00

y1_complete[9,1]          8.00    0.00     0.00       8.00       8.00       8.00       8.00       8.00   800  NaN

y1_complete[9,2]          2.04    0.03     0.83       0.56       1.50       1.98       2.49       3.89   800 1.00

y1_complete[10,1]        30.00    0.00     0.00      30.00      30.00      30.00      30.00      30.00   800  NaN

y1_complete[10,2]         3.15    0.04     1.04       1.24       2.56       3.11       3.73       5.36   800 1.00

y1_complete[11,1]         6.60    0.00     0.00       6.60       6.60       6.60       6.60       6.60     2 1.01

y1_complete[11,2]         1.47    0.03     0.82      -0.12       0.94       1.47       1.98       3.14   800 1.00

 

I highlighted the part where I found a bit odd. These rows are data that I imported via the data block and the model does not need to estimate them, so I expect for all [,1] observations they should have n_eff=800 and Rhat=NaN but for 5 of them the output shows n_eff=2 and some values for Rhat. Is it correct? or there is something wrong with my model?

Thanks in advance,

Jaslene

[Jaslene Lin]

Also, may I ask again about the divergent transitions after warmup, even if I set my control = list(adapt_delta=0.9999999999,max_treedepth=12,stepsize=0.0001) and sometimes I got divergent transitions varying between 1-7, should I be worried about it? Can I trust my posterior distributions for further analyses? I read one of the papers written by stan developers, there it claimed that ideally, there should be no divergent transitions after warmup, however, in reality I did get few, can I get away with it or it is considered to be quite bad and indicating my model is not well constructed?

Thanks in advance,  

Jaslene

[Michael Betancourt]

Jaslene,

The immediate answer is that the Rhat and ESS codes have trouble with constants, so

you shouldn’t worry about those results.  Ultimately this is because there is a very thin

line between a parameter not moving at all (constant) and a parameter moving very, 

very slowly (large autocorrelations).

Note, however, that when you have missing data trying to simulate the missing values 

doesn’t actually do much.  The problem is that you can never recover the missing 

information, so simulating those values doesn’t improve the fit at all.  It is much more

preferable to marginalize the missing values out of the likelihood, for example

p(d_obs | theta) = \int d(d_missing) p(d_obs, d_missing | theta),

and then fit using the marginalized likelihood.

[Michael Betancourt]

I would not advise trusting the posterior.

A divergence indicates that there is some pathological neighborhood in parameter

space that Stan is not able to explore.  That neighborhood might be extremely small,

in which case you could possibly ignore it, or it might be really big, in which case

you absolutely cannot without biasing your fit significantly.

If the pathological neighborhood is small then you should be able to remove the

divergences by decreasing the step size, which is accomplished by increasing

the adapt_depth.  Hence the recommendation in RStan.

But if increasing adapt_delta doesn’t remove the divergences then you have a

large pathological region.  In particular, small increases to adapt_delta let you

explore the neighborhood a bit more, but that just lets you see more pathologies

and hence more divergences pop up.  So the fact that the divergences persist

even as you increase adapt_delta is a sign that there is something seriously

wrong.

If your model is hierarchical then you’ll definitely need the non-centered

parameterization if you’re seeing divergences.

[Bob Carpenter]

> On Aug 8, 2016, at 3:58 AM, Michael Betancourt <betan...@gmail.com> wrote:
>
> Jaslene,
>
> The immediate answer is that the Rhat and ESS codes have trouble with constants, so
> you shouldn’t worry about those results. Ultimately this is because there is a very thin
> line between a parameter not moving at all (constant) and a parameter moving very,
> very slowly (large autocorrelations).
>
> Note, however, that when you have missing data trying to simulate the missing values
> doesn’t actually do much.

In many models, it does absolutely nothing.

> The problem is that you can never recover the missing
> information, so simulating those values doesn’t improve the fit at all.

That depends on the problem and whether you have predictors
or not. Semi-supervised learning's hugely helpful in natural
language processing!

> It is much more
> preferable to marginalize the missing values out of the likelihood, for example
>
> p(d_obs | theta) = \int d(d_missing) p(d_obs, d_missing | theta),
>
> and then fit using the marginalized likelihood.

If they don't contribute to the posterior, they'll drop out as constants;
if they do contribute, this will be the most efficient way to implement
them. But sometimes the integral/sums are too tricky.

- Bob

[Bob Carpenter]

> On Aug 8, 2016, at 4:03 AM, Michael Betancourt <betan...@gmail.com> wrote:
>
> I would not advise trusting the posterior.
>
> A divergence indicates that there is some pathological neighborhood in parameter
> space that Stan is not able to explore. That neighborhood might be extremely small,
> in which case you could possibly ignore it, or it might be really big, in which case
> you absolutely cannot without biasing your fit significantly.
>
> If the pathological neighborhood is small then you should be able to remove the
> divergences by decreasing the step size, which is accomplished by increasing
> the adapt_depth.

adapt_delta parameter (closer to 1) --- it controls target acceptance
rate.

- Bob

[Jaslene Lin]

Thank you all for your suggestions and explanations, I am modifying my model and trying to reduce and keep what is essence in my model. Hopefully, the model would work.

Kindest regards,

Jaslene 

[Jaslene Lin]

Dear Stan developers:

I have made changed to my model which made a huge improvement in terms of the time of simulating and reducing the number of divergent transitions after warmup. However, occasionally, my model can still get divergent transitions after warmup between 1-3, however, if I re-run the simulation again, it might not have divergent transitions. Hence, to get a non-divergent simulation results, I might need to run multiple times, I am wondering is this normal? In all cases, my adapt-acceptance is nearly 1, 0.99999999999. I am not sure whether setting an extremely high adapt_acceptance rate and occasional divergent runs would indicate there is still some problems in my model. Thanks for your insights and explanations.

Kindest regards,

Jaslene 

[Bob Carpenter]

It can lead to biases, but if you only ave a couple occasionally
you should be OK (Michael may say differently).

You should also simulate data from the model and make sure you
can recover the parameters with the expected coverage (e.g., 50% of
parameters fall in the 50% intervals [assuming binomial distribution]).

- Bob

[Michael Betancourt]

Again, the concern here is that the divergences persist even for very small step sizes,

which suggests that the pathological neighborhoods are relevant.  Note also that once

you get past adapt_delta = 0.999 or so it becomes harder and harder to get smaller

step sizes by increasing adapt_delta, so there’s not much else you can do with tuning.

Why would these pathologies be relevant even if you don’t always see divergences?

Intuitively, the number of divergences are proportional to the _boundary_ of the pathological

neighborhood, not it’s volume or probability, which determine how strongly ignoring those

neighborhoods might bias your MCMC estimators.

If you see lots of divergences then you have a pathological neighborhood with a large

boundary that suggests large volume/probability and large biases.  This is clearly bad.

If you see occasional divergences then the pathological neighborhood has a small boundary,

and this could mean two things.  First, the neighborhood could just be small and ignoring

it might be okay.  But it could also just be that the neighborhood is elongated and exposes

just a small surface to the sampler even though it has appreciable probability and could

be important for accurate estimators.

The fact that the divergences persist even for small step sizes would suggest the latter

case, which is why I would still be suspicious and continue to try to improve your model

implementation.  In particular, try to identify where in parameter space the divergences

arise to see if you can figure out what correlations might be causing them.

[Bob Carpenter]

> On Aug 10, 2016, at 6:59 PM, Michael Betancourt <betan...@gmail.com> wrote:
>
> Again, the concern here is that the divergences persist even for very small step sizes,
> which suggests that the pathological neighborhoods are relevant. Note also that once
> you get past adapt_delta = 0.999 or so it becomes harder and harder to get smaller
> step sizes by increasing adapt_delta, so there’s not much else you can do with tuning.
>
> Why would these pathologies be relevant even if you don’t always see divergences?
> Intuitively, the number of divergences are proportional to the _boundary_ of the pathological
> neighborhood, not it’s volume or probability, which determine how strongly ignoring those
> neighborhoods might bias your MCMC estimators.
>
> If you see lots of divergences then you have a pathological neighborhood with a large
> boundary that suggests large volume/probability and large biases. This is clearly bad.
>
> If you see occasional divergences then the pathological neighborhood has a small boundary,
> and this could mean two things. First, the neighborhood could just be small and ignoring
> it might be okay. But it could also just be that the neighborhood is elongated and exposes
> just a small surface to the sampler even though it has appreciable probability and could
> be important for accurate estimators.
>
> The fact that the divergences persist even for small step sizes would suggest the latter
> case, which is why I would still be suspicious and continue to try to improve your model
> implementation. In particular, try to identify where in parameter space the divergences
> arise to see if you can figure out what correlations might be causing them.

Is there a way to figure out where the divergences happen? I only
see displays in RStan of where the initial point is, but the diverging
neighborhood is only indirectly related to these locations.

Can running adaptation longer help? The reason I ask is that if step
sizes / mass matrix is too big, it'll also cause divergences.

- Bob

[Michael Betancourt]

Is there a way to figure out where the divergences happen?   I only
see displays in RStan of where the initial point is, but the diverging
neighborhood is only indirectly related to these locations.

It’s not the initial point that is shown, but the point sampled from the

trajectory prior to proposing the new subtree with the divergence.

It’s still not exactly where the divergence occurred, but it can be closer 

than the initial point.  Also, a divergence is technically a property of

the trajectory and not a single point, so it can be difficult to track back

and blame a single point even if the entire trajectory was saved.

Can running adaptation longer help?  The reason I ask is that if step
sizes / mass matrix is too big, it'll also cause divergences.  

Yes, but it would cause constant divergences, not spurious ones.

Riemannian could help, but in general the mass matrix would have

to be especially poor to cause serious divergence problems.

[Bob Carpenter]

I think the doc may be wrong there! Or at least word on the
street.

But you're right either way --- what we want to see is the whole
path causing the divergence. I'd just like to be sure that they're
not problems where we fully cycle around the level set, then start
spiraling out of control. But I guess the U-turn condition would kick
in if that were the case.

- Bob

[Andrew Gelman]

Soon we'll be able to recommend that people try RHMC on such models, right?

[Bob Carpenter]

If by "soon" you mean roughly Q1 or Q2 2017, then yes.

It's waiting on the refactor, which is getting near to being done.
Then we need to figure out how to do the builds. It's going to
take at least three times as long to compile programs for RHMC
due to higher-order autodiff required and we don't want to have
to pay that cost when we do EHMC, so all the builds and interfaces
are going to need to be enhanced to control what gets compiled.

- Bob