Gibbs sampling with Poisson Gamma models I am trying to obtain a Gibbs sampler for a Poisson-Gamma topic model. Essentially, for each document $d$, the likelihood of $d$ depends on a Poisson parameter $\lambda_d = \sum_k \pi_{k,d}\phi_{k,w}$. In turn, a Gamma prior is assigned to $\pi_{k,d}$ (the strength of topic $k$ in $d$). Likewise, a Gamma prior is assigned to $\phi_{k,w}$ (the frequency of word $w$ in topic $k$).
Can you please show me how to derive a Gibbs sampler for the above model? Is it possible to conveniently collapse?
 A: Most surprisingly, this question and its resolution appear in a question you asked five years ago, in connection with a paper of Gopalan et al. (2015). As I have written my answer before realising this duplicate existed, I'll keep it posted (if only because the original post found no answer at the time).
If$$N_{w,d}\sim\mathcal P(\sum_{1\le k\le K} \pi_{k,d}\phi_{k,w})\tag{1}$$ the likelihood writes as
$$\prod_d\prod_{w \in d}  \frac{1}{n_{w,d}!} (\textstyle \sum_{1\le k\le K} \pi_{k,d}\phi_{k,w})^{n_{w,d}} \cdot e^{-\sum_{1\le k\le K} \pi_{k,d}\phi_{k,w}}$$
and the summation in each term creates a difficulty when simulating the parameters.
A way to avoid the summation is to create latent variables, using the property that a sum of independent Poisson variates is a Poisson variate. If we introduce the latent (unobserved) independent variables $$N_{w,d,k}\sim P(\pi_{k,d}\phi_{k,w})$$ and condition on the fact that $$N_{w,d}=\sum_{1\le k\le K} N_{w,d,k}$$
is their observed sum, then $N_{w,k}$ is distributed as in (1), hence the distribution of $N_{w,d}$ is the marginal of the joint distribution on the $N_{w,d,k}$'s
$$\prod_d\prod_{w \in d}\prod_{k=1}^K \frac{1}{n_{w,d,k}!} [\pi_{k,d}\phi_{k,w}]^{n_{w,d,k}} \, e^{-\pi_{k,d}\phi_{k,w}} \mathbb I_{x_{w,d,1}+\cdots+x_{w,d,K}=x_{w,d}}
$$
Since the $N_{w,d,k}$'s are not observed, they need be simulated as well in the Gibbs sampler. The simulation conditional on $N_{w,d}$ is however straightforward as the conditional distribution of $(N_{w,d,1},\ldots,N_{w,d,K})$ given their sum $N_{w,g}$ is a Multinomial distribution with probability vector
$$(\pi_{1,d}\phi_{k,w},\ldots,\pi_{K,d}\phi_{k,w})\Big/\textstyle \sum_{1\le k\le K} \pi_{k,d}\phi_{k,w}$$
The other components of the Gibbs sampler then turn out to be straightforward as well since
$$p(\pi_{k,d}) \propto \pi_{k,d}^{\textstyle\sum_{w\in d} n_{w,d,k}}\,e^{-\pi_{k,d}\textstyle\sum_{w\in d}\phi_{k,w}}\,p_0(\pi_{k,d})$$
and
$$p(\phi_{k,w})\propto \phi_{k,w}^{\textstyle\sum_{d\ni w} n_{w,d,k}}\,e^{-\phi_{k,w}\textstyle\sum_{d\ni w}\pi_{k,d}}\,p_0(\phi_{k,w})$$
if $p_0$ denotes the appropriate prior pdf and $p$ the full conditional posterior pdf.
In the case when the priors are Poisson-conjugate, i.e. when$$\pi_{k,d}\sim \mathcal G(\alpha,\beta)
\quad\phi_{w,k}\sim\mathcal G(\gamma,\delta)$$
$p(\phi_{k,w})$ corresponds to a
$$\mathcal G\left(\gamma+\sum_{d\ni w} n_{w,d,k},\delta+\sum_{d\ni w}\pi_{k,d}\right)$$density and $p(\pi_{k,d})$ to a $$\mathcal G\left(\alpha+\sum_{w\in d} n_{w,d,k},\beta+\sum_{w\in d}\phi_{k,w}\right)$$density.
