# A Bernoulli mixture model with a Dirichlet prior on the parameters

Let's $$x_1,...,x_N$$ be a set of observation coming from the following generative process:

$$\boldsymbol{\theta} \sim \text{Dirichlet}(\boldsymbol{\alpha})\qquad\boldsymbol{\theta},\boldsymbol\alpha\in\mathbb R_+^K\\ \mathbf{z} \sim \prod_{n=1}^N\text{Multinomial}(1, \boldsymbol{\pi})\qquad\mathbf{z}\in\{1,\ldots,K\}^N\\ \mathbf x \sim \prod_{n=1}^N \text{Bernoulli}(\theta_{z_{n}})\qquad\mathbf{x}\in\{0,1\}^N$$

That is, each observation $$x_k$$ is drawn from one among $$K$$ Bernoulli components. The Bernoulli we draw from is decided by the selector $$\mathbf{z}_{n}$$ that we represent here as a one-hot encoded vector. And the probabilities of the selector pointing to one or another component are encoded in the vector $$\boldsymbol{\pi}$$. The Bernoulli parameters are jointly drawn from a Dirichlet.

My question is: can $$\mathbf{\theta}$$ be marginalized out?

If so, we might be able to get a collapsed Gibbs sampler that only samples the indicator variables $$\mathbf{z}_{n}$$. The problem I see is that the Dirichlet is not a conjugate prior for the Bernoulli. If we put a Beta prior over each Bernoulli parameter, it works out since we have K independent integrals of something nice thanks to conjugacy. But in the Dirichlet case, we can not decompose it in K easy integrals since the parameters $$\theta_1,...,\theta_K$$ are coupled, i.e., they need to sum up to 1.

We have been told that this model can ben marginalized, but we don't find the trick to do this, if any, and it would be useful to have second opinions on whether we are hitting a wall here.

This is the expression we get after writing down the joint probability, ignoring the terms that do not depend on $$\mathbf{\theta}$$ and rearranging the terms:

$$\frac{\Gamma(\sum_k \alpha_k)}{\prod_k\Gamma(\alpha_k)} \int_{\mathfrak S_K}\prod_{k=1}^K \theta_k^{\alpha_k + \sum_{1\le n\le N} x_n\mathbb I_{z_{n}=k}-1} (1-\theta_k)^{\sum_{1\le n\le N} (1-x_n) \mathbb I_{z_{n}=k}} \text{d}\boldsymbol{\theta}$$ where $$\mathfrak S_K$$ denotes the simplex of $$\mathbb R^K$$ and which looks like the integral over a product of Beta distributions where their parameters need to sum up to 1 (because of the Dirichlet prior).

If we denote$$m_k=\sum_{1\le n\le N} x_n\mathbb I_{z_{n}=k}\quad\text{and}\quad r_k=\sum_{1\le n\le N} (1-x_n) \mathbb I_{z_{n}=k}$$ then $$\prod_{k=1}^K \theta_k^{\alpha_k + \sum_{1\le n\le N} x_n\mathbb I_{z_{n}=k}-1} (1-\theta_k)^{\sum_{1\le n\le N} (1-x_n) \mathbb I_{z_{n}=k}}\\=\prod_{k=1}^K \theta_k^{\alpha_k + m_k-1} (1-\theta_k)^{r_k}\\=\prod_{k=1}^K \theta_k^{\alpha_k + m_k-1}\sum_{j_k=0}^{r_k}{r_k\choose j_k}\theta_k^{j_k}\\ =\sum_{j_1=0}^{r_1}\cdots\sum_{j_K=0}^{r_K}\prod_{k=1}^K {r_k\choose j_k}\theta_k^{\alpha_k + m_k-1+j_k}$$ from which the integral can be derived: $$\int_{\mathfrak S_K}\prod_{k=1}^K \theta_k^{\alpha_k + \sum_{1\le n\le N} x_n\mathbb I_{z_{n}=k}-1} (1-\theta_k)^{\sum_{1\le n\le N} (1-x_n) \mathbb I_{z_{n}=k}} \text{d}\boldsymbol{\theta}\\ =\int_{\mathfrak S_K}\sum_{j_1=0}^{r_1}\cdots\sum_{j_K=0}^{r_K}\prod_{k=1}^K {r_k\choose j_k}\theta_k^{\alpha_k + m_k-1+j_k} \text{d}\boldsymbol{\theta}\\ =\sum_{j_1=0}^{r_1}\cdots\sum_{j_K=0}^{r_K}\prod_{k=1}^K {r_k\choose j_k}\int_{\mathfrak S_K}\prod_{k=1}^K\theta_k^{\alpha_k + m_k-1+j_k} \text{d}\boldsymbol{\theta}\\ =\sum_{j_1=0}^{r_1}\cdots\sum_{j_K=0}^{r_K}\prod_{k=1}^K {r_k\choose j_k}\Gamma(\alpha_k + m_k+j_k)\Big/\Gamma\left(\sum_{k=1}^K\alpha_k + m_k+j_k\right)$$ While this is a closed-form expression, it involves $$K\prod_{n=1}^K r_k$$ terms and requires heavy book-keeping. I would rather suggest integrating $$\mathbf z$$ out since the marginal distribution of $$\mathbf z$$ is a standard Bernoulli mixture with known weights. The likelihood is then $$\prod_{n=1}^n \sum_{i=1}^K \pi_i \theta_i^{x_n} (1-\theta_i)^{1-x_n}$$ which can be easily computed as a function of $$\boldsymbol \theta$$ and hence opens the door to a manageable Metropolis-Hastings implementation on the simplex $$\mathfrak S_K$$.
• Thanks @Xi'an. So the key was in the Multinomial theorem! As for your suggestion, I've been following this lead for some hours and I always end up in the same intractable cases. Marginalising out $z$ gives a mixture model where the Bernoulli parameters sum up to one, which differs from a standard one (K Beta priors) and I cannot do Gibbs sampling. The only solution I see is to augment the model by drawing another indicator variable, $c$, with probabilities $\theta$, and then $x_n$ is 1 if $z_n=c_n$. Then you can marginalise out $\theta$. But the resulting Gibbs sampler is not very elegant. – alberto Oct 13 '19 at 13:44
• Concerning the "intractability" of the mixture representation, I beg to disagree. The components of $\boldsymbol\theta$ need not be treated separately in a Gibbs sampler but together in a Metropolis-Hastings version. – Xi'an Oct 13 '19 at 17:02
• Is an indicator function missing in the definition of the $m_k$? – jochen Oct 30 '19 at 8:55