# Divergent transitions when sampling from Dirichlet distribution

**URL:** <https://discourse.mc-stan.org/t/divergent-transitions-when-sampling-from-dirichlet-distribution/1525>\
**Category:** Modeling\
**Created:** [August 8, 2017, 5:18pm UTC](https://discourse.mc-stan.org/t/divergent-transitions-when-sampling-from-dirichlet-distribution/1525 "2017-08-08T17:18:55Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![jeffeaton](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jeffeaton/32/4994_2.png) [@jeffeaton](https://discourse.mc-stan.org/u/jeffeaton)\
**Post date:** [August 8, 2017, 5:18pm UTC](https://discourse.mc-stan.org/t/divergent-transitions-when-sampling-from-dirichlet-distribution/1525/1 "2017-08-08T17:18:55Z")

</div>

Hello all,

I have been trying to debug a model with lots of divergent transitions. The model includes a Dirichlet prior with most of the mass concentrated on one component. In the process of trying to isolate the source of divergent transitions in my model, I found that simply sampling from a Dirichlet distribution resulted in a relatively large number of divergent transitions, but also that this depended on how I ordered the prior.

If I put the parameter with the large concentration as the first parameter in the vector, I get 2-3% divergent transitions. On the other hand, if I put the parameter with the large concentration last, there are no divergent transitions:

```
> ## Sample from Dirichlet with alpha = c(20, 0.1, 0.1, ..., 0.1)
> t(sapply(get_sampler_params(fit1, inc_warmup=FALSE), colMeans))
     accept_stat __stepsize__ treedepth __n_leapfrog__ divergent __energy__
[1,] 0.7990522 0.2981183 3.380 12.230 0.018 17.24487
[2,] 0.8768769 0.2446349 3.560 13.575 0.026 17.23890
[3,] 0.8386421 0.2954739 3.388 12.306 0.034 17.60014
[4,] 0.8900672 0.2370209 3.590 13.564 0.041 17.18930
>
> ## Sample from Dirichlet with alpha = c(0.1, 0.1, ..., 0.1, 20)
> t(sapply(get_sampler_params(fit2, inc_warmup=FALSE), colMeans))
     accept_stat __stepsize__ treedepth __n_leapfrog__ divergent __energy__
[1,] 0.8415861 0.09923680 4.133 23.700 0 17.94246
[2,] 0.9099294 0.08892594 4.387 28.132 0 17.33387
[3,] 0.9150783 0.07699482 4.487 30.936 0 18.04562
[4,] 0.8939827 0.09872973 4.347 27.772 0 17.88692

```

I am wondering if this is ‘expected’ behaviour? I imagine that it may be related to the ‘stick breaking’ algorithm used for the simplex parameterization. Is an intuitive explanation for why it occurs, which might support guidance for most efficient parameterization for more complex models? I also note the caution on page 40 of the manual that high dimensional simplex parameters might require smaller step sizes. Any guidance on the most efficient parameterizations would be useful.

Full stan code for the above example pasted below.

Many thanks,  
Jeff

```
library(rstan)

mod1 <- '
transformed data {
  vector[10] alpha;
  alpha[1] = 20;
  for(i in 2:10)
    alpha[i] = 0.1;
}
parameters {
  simplex[10] x;
}
model {
  x ~ dirichlet(alpha);
}
'

mod2 <- '
transformed data {
  vector[10] alpha;
  for(i in 1:9)
    alpha[i] = 0.1;
  alpha[10] = 20;
}
parameters {
  simplex[10] x;
}
model {
  x ~ dirichlet(alpha);
}
'

fit1 <- stan(model_code = mod1)
fit2 <- stan(model_code = mod2)

## Sampling from Dirichlet with alpha = c(20, 0.1, 0.1, ..., 0.1)
t(sapply(get_sampler_params(fit2, inc_warmup=FALSE), colMeans))

## Sampling from Dirichlet with alpha = c(0.1, 0.1, ..., 0.1, 20)
t(sapply(get_sampler_params(fit1, inc_warmup=FALSE), colMeans))

```

---

<div class="post-metadata">

**Author:** ![bgoodri](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bgoodri/32/4451_2.png) [@bgoodri](https://discourse.mc-stan.org/u/bgoodri)\
**Post date:** [August 8, 2017, 6:00pm UTC](https://discourse.mc-stan.org/t/divergent-transitions-when-sampling-from-dirichlet-distribution/1525/2 "2017-08-08T18:00:07Z")

</div>

> [@jeffeaton](#):
>
> Any guidance on the most efficient parameterizations would be useful.

For this case, the parameterization of the simplex that rescales `gamma(alpha, 1)` primitives avoids the divergences problem and the ordering arbitrariness. See

> [@Ragged array of simplexes](http://discourse.mc-stan.org/t/ragged-array-of-simplexes/1382):
>
> Hello gurus, I would like to implement a model with a ragged array of simplexes. In response to this query about a year ago, Bob mentions (I think) that there was an intention to expose the simplex\_constrain function and log determinant of the Jacobian, which would facilitate this: [https://groups.google.com/forum/#!msg/stan-users/9WTI41kBf4I/sb\_WAt6yBQAJ](https://groups.google.com/forum/#!msg/stan-users/9WTI41kBf4I/sb_WAt6yBQAJ) I was wondering if this was implemented? (Or perhaps I misunderstood.) Many thanks, Jeff

For the case you are probably actually referring to (i.e. with data), it is hard to tell in advance which parameterization will work best.

---

<div class="post-metadata">

**Author:** ![jeffeaton](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jeffeaton/32/4994_2.png) [@jeffeaton](https://discourse.mc-stan.org/u/jeffeaton)\
**Post date:** [August 11, 2017, 2:33pm UTC](https://discourse.mc-stan.org/t/divergent-transitions-when-sampling-from-dirichlet-distribution/1525/3 "2017-08-11T14:33:05Z")

</div>

Hi Ben,

Many thanks for the guidance again. Yep, for the full model with data both approaches – the rescaled gamma or re-ordering the simplex such that the parameter with most concentration is last – reduce the number of divergent transitions to very few. But both also result in adaptation to very small step size and slow sampling.

I guess no free lunch here. I’ll play around a bit more…

Thanks,  
Jeff

---

<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:** [August 14, 2017, 3:23pm UTC](https://discourse.mc-stan.org/t/divergent-transitions-when-sampling-from-dirichlet-distribution/1525/4 "2017-08-14T15:23:10Z")

</div>

> [@bgoodri](#):
>
> Bob mentions (I think) that there was an intention to expose the simplex\_constrain function and log determinant of the Jacobian,

Not yet—still on the to-do list, along with vectorizing our derivatives for the transforms and allowing varying lower and upper bounds in an array.

That’s a quote of a quote, but the original quote doesn’t seem to be there any more!