A Bernoulli mixture model with a Dirichlet prior on the parameters

Question

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).

Answer 1

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$.