#### membership probabilities

So in principle this should give me the posterior membership probabilities, right?

```
generated quantities {
real pi_post[N];
for (i in 1:N) {
pi_post[i] = exp(neg_binomial_2_lpmf(y[i,] | mu[2], phi)) * pi / (exp(neg_binomial_2_lpmf(y[i,] | mu[2], phi)) * pi + exp(neg_binomial_2_lpmf(y[i,] | mu[1], phi)) * (1 - pi));
}
}
```

In your derivation, is \theta = (\mu, \phi)? In contrast to what you have written, I have no term similar to p(\mu, \theta). Is my solution still correct?

#### sampling

The sampling with just one \pi seems to be more problematic than I thought previously. I am not sure whether I changed anything, but now I regularly get results like (20,8), (19, 11) for the \mu s (instead of (20,5)) if I fix \pi = 0.75. So now the smaller \mu seems affected.

If I donâ€™t fix \pi but have it estimated as a parameter, than I get tons of divergent transitions. Even in the scenario with the fixed \pi I have the feeling that the effective sample size is rather low for all parameters.

Here again the full code:

```
generatemixedcounts <- function(n = 6, N = 50, phi = 1, mu1 = 20, mu2 = 5){
## select some of the entries of the matrix
m <- matrix(NA, nrow = N, ncol = n)
m[] <- rnbinom(n = n*N, mu = mu1, size = phi)
select <- sample(c(T,F), size = N, replace = T, prob = c(0.25, 0.75))
## assign rnbinom with different mu
m[select, ] <- rnbinom(n = sum(select) * n, mu = mu2, size = phi)
return(m)
}
toy <- generatemixedcounts(N = 60, mu1 = 20, mu2 = 5, n = 6, phi = 1)
standata <- list(y = toy, N = nrow(toy), n = ncol(toy))
write("
data {
int<lower=0> N;
int<lower=0> n;
int<lower=0> y[N,n];
}
parameters {
real <lower=0, upper=1> pi;
positive_ordered[2] mu;
real<lower=0> phi;
}
model {
// real pi = 0.75;
vector[2] contributions;
for (j in 1:n){
for (i in 1:N) {
contributions[1] = log(pi) + neg_binomial_2_lpmf(y[i,j] | mu[2], phi);
contributions[2] = log(1 - pi) + neg_binomial_2_lpmf(y[i,j] | mu[1], phi);
target += log_sum_exp(contributions);
}
}
}
// generated quantities {
// real pi = 0.25;
// real pi_post[N];
// for (i in 1:N) {
// pi_post[i] = exp(neg_binomial_2_lpmf(y[i,] | mu[2], phi)) * pi / (exp(neg_binomial_2_lpmf(y[i,] | mu[2], phi)) * pi + exp(neg_binomial_2_lpmf(y[i,] | mu[1], phi)) * (1 - pi));
// }
// }
", file = "./mixturenegbin.stan")
stanfit <- rstan::stan(file = "./mixturenegbin.stan",
data = standata,
iter = 1000, warmup = 500, thin = 1)
stanfit
```

Thanks a lot!