# Dirichlet\_log: probabilities is not a valid simplex

**URL:** <https://discourse.mc-stan.org/t/dirichlet-log-probabilities-is-not-a-valid-simplex/1805>\
**Category:** Modeling\
**Created:** [September 6, 2017, 5:59am UTC](https://discourse.mc-stan.org/t/dirichlet-log-probabilities-is-not-a-valid-simplex/1805 "2017-09-06T05:59:12Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![kbchoi](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/kbchoi/32/861_2.png) [@kbchoi](https://discourse.mc-stan.org/u/kbchoi)\
**Post date:** [September 6, 2017, 5:59am UTC](https://discourse.mc-stan.org/t/dirichlet-log-probabilities-is-not-a-valid-simplex/1805/1 "2017-09-06T05:59:12Z")

</div>

Hi,

I am trying to fit a 4-component mixture model where the probability of component membership of each sample ( **pi** ) is estimated in a separate stan fit. I am using dirichlet to introduce uncertainty but getting many of “Exception thrown at line 20: dirichlet\_log: probabilities is not a valid simplex. sum(probabilities) = nan, but should be 1”. I am using the first component to capture samples whose _n_ is zero. Where do you think I am getting nan? Any comment would be greatly appreciated.

```
data {
    int<lower=1> N;
    int<lower=0> n[N];
    vector<lower=0,upper=1>[4] pi[N];
}
parameters {
    real<lower=0,upper=1> a4;
    real<upper=max(n)> a5;
    real<lower=1> b5;
    real<upper=max(n)> a6;
    real<lower=1> b6;
    real<upper=max(n)> a7;
    real<lower=1> b7;
}
model {
    for (i in 1:N) {
        vector[4] lps;
        vector[4] tau;
        tau ~ dirichlet(pi[i]);
        lps[1] = tau[1] * exp(neg_binomial_2_lpmf(n[i]|a4, 1));
        lps[2] = tau[2] * exp(neg_binomial_2_lpmf(n[i]|a5, b5));
        lps[3] = tau[3] * exp(neg_binomial_2_lpmf(n[i]|a6, b6));
        lps[4] = tau[4] * exp(neg_binomial_2_lpmf(n[i]|a7, b7));
        target += log(sum(lps));
    }
}

```

---

<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:** [September 6, 2017, 6:12am UTC](https://discourse.mc-stan.org/t/dirichlet-log-probabilities-is-not-a-valid-simplex/1805/2 "2017-09-06T06:12:02Z")

</div>

> [@kbchoi](#):
>
> Where do you think I am getting nan?

from the entire `model` block:

```no-highlight
    for (i in 1:N) {
        vector[4] lps;
        vector[4] tau;
        tau ~ dirichlet(pi[i]);
        lps[1] = tau[1] * exp(neg_binomial_2_lpmf(n[i]|a4, 1));
        lps[2] = tau[2] * exp(neg_binomial_2_lpmf(n[i]|a5, b5));
        lps[3] = tau[3] * exp(neg_binomial_2_lpmf(n[i]|a6, b6));
        lps[4] = tau[4] * exp(neg_binomial_2_lpmf(n[i]|a7, b7));
        target += log(sum(lps));
    }

```

First of all, `tau` is undefined in this loop. You need something like

```no-highlight
simplex[4] tau[N];

```

in the `parameters` block.

Once you fix that, you have to use numerically stable constructions like:

```no-highlight
        lps[1] = log(tau[i,1]) + neg_binomial_2_lpmf(n[i]|a4, 1);
        lps[2] = log(tau[i,2]) + neg_binomial_2_lpmf(n[i]|a5, b5);
        lps[3] = log(tau[i,3]) + neg_binomial_2_lpmf(n[i]|a6, b6);
        lps[4] = log(tau[i,4]) + neg_binomial_2_lpmf(n[i]|a7, b7);
        target += log_sum_exp(lps);

```

---

<div class="post-metadata">

**Author:** ![kbchoi](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/kbchoi/32/861_2.png) [@kbchoi](https://discourse.mc-stan.org/u/kbchoi)\
**Post date:** [September 6, 2017, 6:17am UTC](https://discourse.mc-stan.org/t/dirichlet-log-probabilities-is-not-a-valid-simplex/1805/3 "2017-09-06T06:17:07Z")

</div>

I appreciate your response. One quick related general question. Do I need to declare _tau_ in the `parameters` block even if I am not trying to estimate it?

---

<div class="post-metadata">

**Author:** ![kbchoi](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/kbchoi/32/861_2.png) [@kbchoi](https://discourse.mc-stan.org/u/kbchoi)\
**Post date:** [September 6, 2017, 7:06am UTC](https://discourse.mc-stan.org/t/dirichlet-log-probabilities-is-not-a-valid-simplex/1805/4 "2017-09-06T07:06:17Z")

</div>

Now I am getting “dirichlet\_log: prior sample sizes[1] is 0, but must be \> 0!” with the following fix.

```
data {
    int<lower=1> N;
    int<lower=0> n[N];
    vector<lower=0,upper=1>[4] pi[N]; // data shape should be (N, 4)
}
parameters {
    real<lower=0,upper=1> a4;
    real<upper=max(n)> a5;
    real<lower=1> b5;
    real<upper=max(n)> a6;
    real<lower=1> b6;
    real<upper=max(n)> a7;
    real<lower=1> b7;
    simplex[4] tau[N];
}
model {
    for (i in 1:N) {
        tau[i] ~ dirichlet(pi[i]);
    }
    for (i in 1:N) {
        vector[4] lps;        
        lps[1] = log(tau[i,1]) + neg_binomial_2_lpmf(n[i]|a4, 1);
        lps[2] = log(tau[i,2]) + neg_binomial_2_lpmf(n[i]|a5, b5);
        lps[3] = log(tau[i,3]) + neg_binomial_2_lpmf(n[i]|a6, b6);
        lps[4] = log(tau[i,4]) + neg_binomial_2_lpmf(n[i]|a7, b7);
        target += log_sum_exp(lps);
    }
}

```

---

<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:** [September 6, 2017, 1:36pm UTC](https://discourse.mc-stan.org/t/dirichlet-log-probabilities-is-not-a-valid-simplex/1805/5 "2017-09-06T13:36:03Z")

</div>

> [@kbchoi](#):
>
> dirichlet\_log: prior sample sizes[1] is 0, but must be \> 0!

How many times? Underflowing `tau` to 0 is not as big a deal.

---

<div class="post-metadata">

**Author:** ![kbchoi](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/kbchoi/32/861_2.png) [@kbchoi](https://discourse.mc-stan.org/u/kbchoi)\
**Post date:** [September 6, 2017, 1:54pm UTC](https://discourse.mc-stan.org/t/dirichlet-log-probabilities-is-not-a-valid-simplex/1805/6 "2017-09-06T13:54:53Z")

</div>

Many many times. If the message implies underflowing, I think it is all right since my pi’s contains lots of zeros although they sum to one. The fitting did not finish with the following message at the end.

```
Initialization between (-2, 2) failed after 100 attempts. 
 Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.

```

I guess having zeros in dirichlet parameter is not good?

---

<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:** [September 6, 2017, 2:14pm UTC](https://discourse.mc-stan.org/t/dirichlet-log-probabilities-is-not-a-valid-simplex/1805/7 "2017-09-06T14:14:18Z")

</div>

Yes, the concentration parameters should be positive although they do not have to sum to 1.

---

<div class="post-metadata">

**Author:** ![kbchoi](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/kbchoi/32/861_2.png) [@kbchoi](https://discourse.mc-stan.org/u/kbchoi)\
**Post date:** [September 6, 2017, 5:21pm UTC](https://discourse.mc-stan.org/t/dirichlet-log-probabilities-is-not-a-valid-simplex/1805/8 "2017-09-06T17:21:39Z")

</div>

The fitting finishes successfully with the following change although it may have introduced some bias.

```
tau[i] ~ dirichlet(pi[i]+0.000001);

```

Thanks Ben!

---

<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:** [September 10, 2017, 12:27pm UTC](https://discourse.mc-stan.org/t/dirichlet-log-probabilities-is-not-a-valid-simplex/1805/9 "2017-09-10T12:27:28Z")

</div>

Adding small epsilons is rarely what you want to do. How were those `pi` set in the first place? You should just set them in the data rather than setting them zero then adding epsilon. You should also run sensitivity analyses to see how much they affect inference—these things can sometimes have a suprrisingly strong effect.

When you have Dirichlet parameters less than one, the prior is putting most of its mass on solutions with zeros (technically, very very small non-zero values that round to zero because 1e-350 is roughly the smallest number that can be represented).

Negative binomials are already hard to fit—I imagine making a mixture of them is going to be extra difficult unless the clusters are well separated. Michael Betancourt wrote a case study on mixtures (see web site \>\> users \>\> documentation \>\> case studies).

You can simplify this a bit:

```
vector[4] lps;        
lps[1] = log(tau[i,1]) + neg_binomial_2_lpmf(n[i]|a4, 1);
lps[2] = log(tau[i,2]) + neg_binomial_2_lpmf(n[i]|a5, b5);
lps[3] = log(tau[i,3]) + neg_binomial_2_lpmf(n[i]|a6, b6);
lps[4] = log(tau[i,4]) + neg_binomial_2_lpmf(n[i]|a7, b7);

```

to

```
vector[4] lps = log(tau[i]);        
lps[1] += neg_binomial_2_lpmf(n[i]|a4, 1);
lps[2] += neg_binomial_2_lpmf(n[i]|a5, b5);
lps[3] += neg_binomial_2_lpmf(n[i]|a6, b6);
lps[4] += neg_binomial_2_lpmf(n[i]|a7, b7);

```

I don’t understand the naming of the `a` and `b`. Why not follow the ordering and make those

```
vector<lower = 1> b[3];
vector<upper = max(n)> a[4];

```

You should also give the `a` and `b` priors—like othe precision-like parameters, they can tend to drift upwards because ever higher values correspond to ever lower variance.

---

<div class="post-metadata">

**Author:** ![betanalpha](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/betanalpha/32/4_2.png) [@betanalpha](https://discourse.mc-stan.org/u/betanalpha)\
**Post date:** [September 10, 2017, 1:55pm UTC](https://discourse.mc-stan.org/t/dirichlet-log-probabilities-is-not-a-valid-simplex/1805/10 "2017-09-10T13:55:22Z")

</div>

This is ultimately a pathology of the empirical Bayes strategy of plugging a point estimate from some previous fit into the priors. It’s just made worse by the fact that the point estimates are themselves lie on the corners of parameter space. Adding epsilon to all of the pi regularizes the point estimates and avoids the immediate pathology, but one should still be skeptical of the validity of the resulting fit.

---

<div class="post-metadata">

**Author:** ![kbchoi](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/kbchoi/32/861_2.png) [@kbchoi](https://discourse.mc-stan.org/u/kbchoi)\
**Post date:** [September 11, 2017, 2:13am UTC](https://discourse.mc-stan.org/t/dirichlet-log-probabilities-is-not-a-valid-simplex/1805/11 "2017-09-11T02:13:51Z")

</div>

I am dealing with a set of binomial (n\_i, p\_i) samples where there are four different classes according to p\_i’s. Let’s assume they cluster pretty well. After I estimate the expectation on 4-dim latent variable ( **pi** ), I wanted to get the shrinkage estimators for n\_i’s based upon the probabilistic class membership of each sample. In order to implement this idea, I had to separate classification and shrinkage estimation. Do you think there is a better way than using dirichlet? Do you think there’s a way to combine the classification and the shrinkage in one stan fitting? I would really appreciate any comment.

---

<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:** [September 11, 2017, 10:57am UTC](https://discourse.mc-stan.org/t/dirichlet-log-probabilities-is-not-a-valid-simplex/1805/12 "2017-09-11T10:57:09Z")

</div>

> [@kbchoi](#):
>
> Do you think there’s a way to combine the classification and the shrinkage in one stan fitting?

Not only is the answer “yes”, it’s the standard approach in Bayesian hierarchical modeling.

The Dirichlet is a very weak prior in that it only gives you a way to model location and a kind of precision, not anything like covariance. A multivariate logistic normal is more flexible. You can see that in the “Correlated topic model” section of the manual (in the latent Dirichlet allocation section of the “Clustering Models” chapter).

---

<div class="post-metadata">

**Author:** ![kbchoi](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/kbchoi/32/861_2.png) [@kbchoi](https://discourse.mc-stan.org/u/kbchoi)\
**Post date:** [September 12, 2017, 12:58pm UTC](https://discourse.mc-stan.org/t/dirichlet-log-probabilities-is-not-a-valid-simplex/1805/13 "2017-09-12T12:58:16Z")

</div>

@Bob_Carpenter Here’s my original problem. Because I am trying to “reuse” **pi** , I don’t feel this is a standard situation for stan. Any advice would be great.

```
data {
    int<lower=1> N;
    int<lower=0> n[N];
    int<lower=0> x[N];
}
parameters {
    simplex[4] pi;
    real<lower=1> a[4];
    real<lower=1> b[4];
    real<upper=max(n)> c[4];
    real<lower=1> d[4];
}
model {
    vector[4] log_pi = log(pi);
    for (i in 1:N) {
        vector[4] lps = log_pi;
        for (j in 1:4)
            lps[j] += beta_binomial_lpmf(x[i]|n[i], a[j], b[j]);
        target += log_sum_exp(lps);
    }

    //
    // Now I also want to fit n's to estimate c's and d's using the same **pi** ,
    // but I have already used "target" to fit x above.
    //
    // for (i in 1:N) {
    // vector[4] lps2 = log_pi;
    // for (j in 1:4)
    // lps2[j] += neg_binomial_2_lpmf(n[i]|c[j], d[j]);
    // target += log_sum_exp(lps2);
    //

}

```

---

<div class="post-metadata">

**Author:** ![stijn](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/stijn/32/105_2.png) [@stijn](https://discourse.mc-stan.org/u/stijn)\
**Post date:** [September 12, 2017, 1:25pm UTC](https://discourse.mc-stan.org/t/dirichlet-log-probabilities-is-not-a-valid-simplex/1805/14 "2017-09-12T13:25:01Z")

</div>

It’s not a problem to use target more than once. You would just be adding more components to the log likelihood. It sounds like you have two equations per observation which is fine. It can help to identify the pi’s but it can also lead to a more complex geometry. I would lookout for divergences.

---

<div class="post-metadata">

**Author:** ![kbchoi](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/kbchoi/32/861_2.png) [@kbchoi](https://discourse.mc-stan.org/u/kbchoi)\
**Post date:** [September 12, 2017, 2:12pm UTC](https://discourse.mc-stan.org/t/dirichlet-log-probabilities-is-not-a-valid-simplex/1805/15 "2017-09-12T14:12:34Z")

</div>

@stijn It takes a long time and did not converge. In fact, I don’t want the neg\_binomial components to influence the estimation of **pi** and that is the cause of my headache. I guess it is also related to a more general question: **What are my options when n is not fixed in my binomial mixture data?**

---

<div class="post-metadata">

**Author:** ![stijn](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/stijn/32/105_2.png) [@stijn](https://discourse.mc-stan.org/u/stijn)\
**Post date:** [September 13, 2017, 6:40am UTC](https://discourse.mc-stan.org/t/dirichlet-log-probabilities-is-not-a-valid-simplex/1805/16 "2017-09-13T06:40:40Z")

</div>

I am just trying to think generatively here. This may or may not be helpful.

Your data generation process is something like the following. You have 4 different groups. In each group n ~ neg\_binomial(c\_j, d\_j) and x ~ beta\_binomial(n, a\_j, b\_j) with a\_j, b\_j, d\_j \> 0 and c\_j \< max(n).

You can already use that information to simulate x’s and n’s outside of Stan. First simulate groups of 4 a\_j, b\_j, c\_j, d\_j, then simulate n, and next x. You should then check whether these simulations are reasonable for your data. If not, you can start thinking about what is wrong with how you have written the model so far. For instance, you probably need more informative priors on a, b, c, and d. You have a mixture model with discrete distribution. Those are red flags for identifiability problems.

> [@kbchoi](#):
>
> In fact, I don’t want the neg\_binomial components to influence the estimation of pi and that is the cause of my headache.

That made me think that maybe you expect the 4 groups to have different a\_j, b\_j but roughly similar c\_j and d\_j. If that is the case you do not need the mixture model for the neg\_binomial. You can just work with a singular c\_0 and d\_0. When the neg\_binomial component does not help to identify the groups (my previous assumption) and you are not interested in the c’s and d’s, you may not even need to estimate the neg\_binomial components. (I am probably wrong in my assumptions but I am trying to understand why you want the neg\_binomial in there).

> [@kbchoi](#):
>
> What are my options when n is not fixed in my binomial mixture data?

What do you mean by ‘not fixed’? n varies per observation or n has measurement error or n varies per group or something else?

---

<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:** [September 13, 2017, 11:09am UTC](https://discourse.mc-stan.org/t/dirichlet-log-probabilities-is-not-a-valid-simplex/1805/17 "2017-09-13T11:09:53Z")

</div>

> [@kbchoi](#):
>
> I don’t want the neg\_binomial components to influence the estimation of pi and that is the cause of my headache.

This is just how Bayesian inference works. BUGS includes a hack called “cut” that breaks Bayesian inference by preventing information flow between random variables. We intentionally avoided that in Stan.

> [@stijn](#):
>
> I am just trying to think generatively here.

This is the way to think about things!

---

<div class="post-metadata">

**Author:** ![kbchoi](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/kbchoi/32/861_2.png) [@kbchoi](https://discourse.mc-stan.org/u/kbchoi)\
**Post date:** [September 13, 2017, 11:42am UTC](https://discourse.mc-stan.org/t/dirichlet-log-probabilities-is-not-a-valid-simplex/1805/18 "2017-09-13T11:42:56Z")

</div>

I was trying to test whether n’s are different between clusters by comparing their confidence intervals, along the same sampled **pi**. I still needed to estimate c’s and d’s for that purpose. But I also found I was not thinking generatively. I appreciate your help @stijn and @Bob_Carpenter.

---

<div class="post-metadata">

**Author:** ![stijn](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/stijn/32/105_2.png) [@stijn](https://discourse.mc-stan.org/u/stijn)\
**Post date:** [September 14, 2017, 2:52am UTC](https://discourse.mc-stan.org/t/dirichlet-log-probabilities-is-not-a-valid-simplex/1805/19 "2017-09-14T02:52:46Z")

</div>

> [@kbchoi](#):
>
> test whether n’s are different between clusters by comparing their confidence intervals, along the same sampled pi. I still needed to estimate c’s and d’s for that purpose

I think you can do that in the generated quantities block or outside of stan. You don’t need the c’s and d’s or the negative binomial in your model.

1. Simulate group membership ~ categorial\_rng(pi) for each observation in each draw
2. Calculate the mean or median of n or 80% interval or whatever floats your boat for each simulated group
3. For each simulation, check whether the mean of group 1 \> group 2 (or any other comparison you want to make).
4. Interpret the distribution of the variable you calculate in 3. E.g. if in 50% of the draws group 1 \> group 2 and in 50% group 2 \> group 1, the groups have roughly the same mean n.

---

<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:** [September 15, 2017, 4:53pm UTC](https://discourse.mc-stan.org/t/dirichlet-log-probabilities-is-not-a-valid-simplex/1805/20 "2017-09-15T16:53:59Z")

</div>

> [@kbchoi](#):
>
> comparing their confidence intervals,

Stan’s generating posterior intervals rather than confidence intervals.

> [@stijn](#):
>
> Simulate group membership ~ categorial\_rng(pi) for each observation in each draw  
> Calculate the mean or median of n or 80% interval or whatever floats your boat for each simulated group

You can usually just do these categorical estimates in expectation directly rather than simulating. That is, you get the weight of `pi` in each iteration for the expectation calculation. Same expected result as simulating, but more efficient due to Rao-Blackwell.

You can’t do this easily for continuous examples because the simple sums turn into integrals—so for continuous cases the RNG approach will usually be the best way to code it.

[Next page](https://discourse.mc-stan.org/t/dirichlet-log-probabilities-is-not-a-valid-simplex/1805.md?page=2)