# Why do chains get stuck, and how can I diagnose this behaviour?

**URL:** <https://discourse.mc-stan.org/t/why-do-chains-get-stuck-and-how-can-i-diagnose-this-behaviour/37505>\
**Category:** Modeling\
**Tags:** diagnostics\
**Created:** [December 4, 2024, 9:30pm UTC](https://discourse.mc-stan.org/t/why-do-chains-get-stuck-and-how-can-i-diagnose-this-behaviour/37505 "2024-12-04T21:30:29Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![PhDemetri](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/phdemetri/32/3273_2.png) [@PhDemetri](https://discourse.mc-stan.org/u/PhDemetri)\
**Post date:** [December 4, 2024, 9:30pm UTC](https://discourse.mc-stan.org/t/why-do-chains-get-stuck-and-how-can-i-diagnose-this-behaviour/37505/1 "2024-12-04T21:30:29Z")

</div>

Consider the following Beta Binomial Model

```stan
data{
  int n;
  array[n] int N;
  array[n] int y;
}
parameters{
  real<lower=0, upper=1> mu;
  real<lower=0> kappa;
}
model{
  kappa ~ cauchy(0, 1);
  mu ~ beta(1, 1);
  y ~ beta_binomial(N, mu/kappa, (1-mu)/kappa);
}
generated quantities{
  array[n] int ppc = beta_binomial_rng(N, mu/kappa, (1-mu)/kappa);
}

```

Suppose I had 100 seperate binomial experiments, each of which had the same number of trials. Unbeknownst to me, each trial had an identical probability of success

```r

N <- rpois(100, 50)
y <- rbinom(100, N, 0.1)

```

Given these data, I know that `kappa` should be 0, and so I expect there to be divergences from this model.

```r
library(cmdstanr)
library(cmdstanr)
library(tidyverse)
library(marginaleffects)
library(bayesplot)
library(tidybayes)

N <- rpois(100, 50)
y <- rbinom(100, N, 0.1)

code <- "
data{
  int n;
  array[n] int N;
  array[n] int y;
}
parameters{
  real<lower=0, upper=1> mu;
  real<lower=0> kappa;
}
model{
  kappa ~ cauchy(0, 1);
  mu ~ beta(1, 1);
  y ~ beta_binomial(N, mu/kappa, (1-mu)/kappa);
}
generated quantities{
  array[n] int ppc = beta_binomial_rng(N, mu/kappa, (1-mu)/kappa);
}
"

data <- list(N=N, y=y, n=length(N))
model <- code %>% write_stan_file() %>% cmdstan_model()
fit <- model$sample(
  data=data, 
  parallel_chains = parallel::detectCores(),
  seed = 0)

fit$diagnostic_summary()
Warning: 1123 of 4000 (28.0%) transitions ended with a divergence.
See https://mc-stan.org/misc/warnings for details.

Warning: 2877 of 4000 (72.0%) transitions hit the maximum treedepth limit of 10.
See https://mc-stan.org/misc/warnings for details.

$num_divergent
[1] 131 98 710 184

$num_max_treedepth
[1] 869 902 290 816

$ebfmi
[1] 2.009023 1.993664 1.970364 2.044654

```

Not surprised. However, when I visualize the chains, its clear that some of them get “stuck”

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/3X/a/1/a1cac5136f6134faa6e3bc2d3d0faa991d84fd08.png)

How would I diagnose why these chains get stuck? The beta binomial model above is simple but just an example. I know why the model is pathological, but I might not know why in some other model in which I experience this behaviour.

When I see chains do this, how can I diagnose what is happening and where in my model I should consider reparameterizing?

EDIT: The same behaviour occurs even when I specify the beta binomial parameters as `mu*kappa, (1-mu)*kappa`

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [December 4, 2024, 10:05pm UTC](https://discourse.mc-stan.org/t/why-do-chains-get-stuck-and-how-can-i-diagnose-this-behaviour/37505/2 "2024-12-04T22:05:12Z")

</div>

> [@PhDemetri](#):
>
> How would I diagnose why these chains get stuck?

They almost always get stuck when they enter a region of high curvature and the step size is too large to escape easily. Usually when they’re stuck, you get divergent transitions.

I’m used to seeing this done as `beta_binomial(N, mu * kappa2, (1 - mu) * kappa2)` with a prior that penalizes large values of `kappa2`. You’ve done it the other way around, which I understand is more consistent with your setting of expecting `kappa` = 0. I suspect that’ll be more stable than trying to run `kappa` down to zero, which is likely to introduce high curvature.

This is a general problem that @andrewgelman has been trying to solve since before we started working on Stan—he calls it the “unfolding flower”. The idea is that you want your model to expand as the data expands—from say, a binomial to a beta-binomial. But you want the more general model to fall back to the simpler model when the data’s simpler. For that, I like the penalized complexity priors from the INLA gang, but they’re done one-off for particular problems in that setting.

---

<div class="post-metadata">

**Author:** ![PhDemetri](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/phdemetri/32/3273_2.png) [@PhDemetri](https://discourse.mc-stan.org/u/PhDemetri)\
**Post date:** [December 4, 2024, 10:11pm UTC](https://discourse.mc-stan.org/t/why-do-chains-get-stuck-and-how-can-i-diagnose-this-behaviour/37505/3 "2024-12-04T22:11:51Z")

</div>

@Bob_Carpenter

> I’m used to seeing this done as `beta_binomial(N, mu * kappa2, (1 - mu) * kappa2)` with a prior that penalizes large values of `kappa2`

Right. I wrote this model in this way specifically to create this type of behaviour. However, I can replicate this behaviour even under the parameterization you provide. The problem now is the prior, and if I switch to something less heavy tailed than a Cauchy, the model samples well.

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/3X/d/c/dcabf5f0bdbddc47330da0ca60dfe35f4f8a91e4.png)

So the problem is curvature and stepsize. When I encounter this, do you recommend then I try decreasing step size (via adapt\_delta I believe?) or is the parameterization a better approach? I anticipate [the folk theorem of statistical computing](https://statmodeling.stat.columbia.edu/2008/05/13/the_folk_theore/) applies.

Is there a way to diagnose what part fo the model might be a target for reparameterization? Somewhat selfishly, I’m experiencing this behaviour in another model (not in Stan) but could only replicating the behaviour in a minimal working example using Stan.

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [December 4, 2024, 10:15pm UTC](https://discourse.mc-stan.org/t/why-do-chains-get-stuck-and-how-can-i-diagnose-this-behaviour/37505/4 "2024-12-04T22:15:34Z")

</div>

Any time the data want the parameters to lie on the boundary of parameter space, you’re going to have problems. Here you see in my parameterization that `kappa2` tries to go to +infinity, whereas in yours, `kappa` tries to go to zero. Both are big warning flags as that’s the boundary of parameter space. It usually means you want a simpler model that just assumes the boundary and that’s what we usually do to mitigate this kind of problem.

> [@PhDemetri](#):
>
> Is there a way to diagnose what part fo the model might be a target for reparameterization?

Not directly that I know of, but you can plot in Stan where the divergences arise. Those are often a clue as to where problems arise.

---

<div class="post-metadata">

**Author:** ![PhDemetri](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/phdemetri/32/3273_2.png) [@PhDemetri](https://discourse.mc-stan.org/u/PhDemetri)\
**Post date:** [December 4, 2024, 10:16pm UTC](https://discourse.mc-stan.org/t/why-do-chains-get-stuck-and-how-can-i-diagnose-this-behaviour/37505/5 "2024-12-04T22:16:34Z")

</div>

Appreciated!