# Bayesian updates for Dirichlet-multinomial with Gamma prior

Let \begin{aligned} X_i &\sim \text{Dir-multinom}(X\mid\lambda)\\ \lambda_{j} &\sim \text{Gamma}(\lambda_j\mid\alpha,\beta)\\ \end{aligned} where $$i$$ iterates over observations, $$j$$ iterates over categories.

I'm trying to perform inference on $$\lambda$$. My thought was to calculate the latent $$\pi$$, such that the model becomes: \begin{aligned} X_i &\sim \text{Multinomial}(X_i\mid\pi_i)\\ \pi_i &\sim \text{Dirichlet}(\pi\mid\lambda)\\ \lambda_j &\sim \text{Gamma}(\lambda_j\mid\alpha,\beta) \end{aligned} Then posterior updating for $$\pi$$ in the Gibbs sampler becomes: $$\pi_i \mid X_i, \lambda \sim \text{Dirichlet}(\pi_i\mid \lambda + X_i)$$ Which is frequently sampled as $$\rho_{ij}\mid X_{ij}, \lambda_j \sim \text{Gamma}(\rho_{ij}\mid \lambda_{j} + X_{ij}, 1)$$ then Gamma $$\rho_i$$ is reduced to Dirichlet $$\pi_i$$ by projecting onto the simplex (divide by the sum). With these $$\rho_i$$'s, I should be able to produce estimates for $$\lambda$$ as owing to the second model formulation, the full conditional for $$\lambda$$ is proportional to $$f(\lambda_j\mid\rho) \propto \prod_{i = 1}^{n}\text{Gamma}(\rho_{ij}\mid\lambda_j, 1)\times \text{Gamma}(\lambda_j\mid\alpha,\beta)$$ then sampling $$\lambda$$ using Metropolis Hastings.

However, when I implement this sampler, my estimates for $$\lambda$$ and subsequently $$\rho$$ just grow larger the longer I run the sampler.

I suppose I should mention that I can use the first model formulation; updating is slower as it can't be done marginally. Is there any valid reason the second model can not be implemented? I have reason for wanting the latent $$\pi$$'s/$$\rho$$'s.

• There is no apparent reason why Gibbs sampling shouldn't work in this case. Have you considered the possibility that there is a mistake in your implementation ? (did you test it with different starting points?) May 10, 2022 at 18:12
• There's always the possibility of a bug in the code. But, this model is really simple, the multinomial dirichlet updating is well documented, and the Gamma Gamma inference steps use the same functions I wrote for previous Gamma Gamma inference models that worked. I don't see where an issue could have inserted itself. May 10, 2022 at 18:58
• Looking at the update equations it does look like $\lambda$ should grow indefinitely because the $\rho$'s are sampled from a Gamma with mean $\lambda+X > \lambda$ which in turn increases the mean of $f(\lambda|\rho)$. So I wonder if using $\text{Gamma}(\rho|\lambda,1)$ instead of $\text{Dir}(\pi|\lambda)$ is actually correct here May 10, 2022 at 19:30

Let's look carefully at the conditional distribution of $$\lambda$$ : Using the original parametrization

$$P(\lambda | \pi ) \propto \prod_i \text{Dirichlet}(\pi_i | \lambda) \times \prod_j \text{Gamma}(\lambda_j|\alpha,\beta) \\ = \left(\frac{1}{B(\lambda)}\right)^n \prod_{i,j} \pi_{ij}^{\lambda_j-1} \times \prod_j \text{Gamma}(\lambda_j|\alpha,\beta)$$

We can see immediately that the $$\lambda_j$$'s are not independent conditional on $$\pi_{ij}$$ (due to to the $$B(\lambda)$$ factor).

Alternatively using the $$\rho$$ parametrization we can express the model as :

$$\rho_{ij} \sim \text{Gamma}(\lambda_j,1)$$ $$X_i \sim \text{Multinomial}\left(\frac{\rho_i}{\sum_j \rho_{ij}} \right)$$

However now the posterior probability of $$\rho$$ is

$$P(\rho|X,\lambda) \propto \prod_{ij} \left(\frac{\rho_{ij}}{\sum_j \rho_{ij}}\right)^{X_{ij}} \times \text{Gamma}(\rho_{ij} | \lambda_j,1) \\ = \prod_i \frac{1}{(\sum_j \rho_{ij})^{\sum_j X_{ij}}} \prod_j \rho_{ij}^{X_{ij}} \times \text{Gamma}(\rho_{ij}|\lambda_j,1)$$

which is not a product of independent Gamma distributions (due to the $$\sum_j \rho_{ij}$$ term).

So although now the $$\lambda_j$$'s are independent conditional on $$\rho_{ij}$$, the $$\rho_{ij}$$'s themselves are not independent conditional on the data. Therefore in either way we can't fully decouple the dependencies between the categories as you attempted to do.

To facilitate the analysis, let $$\mathbf{x} = (\mathbf{x}_1,...,\mathbf{x}_n)$$ denote the matrix of observed sample vectors and let $$\mathbf{x}_i = (x_{i,1},...,x_{i,K})$$ denote the individual sample vectors. We will also let $$\dot{\lambda} \equiv \sum_{k=1}^K \lambda_k$$ denote the parameter sum. The likelihood function in this case is:

\begin{align} L_\mathbf{x}(\boldsymbol{\lambda}) &= \prod_{i=1}^n \text{Dir-Mu}(\mathbf{x}_i | \boldsymbol{\lambda}) \\[6pt] &\propto \frac{\Gamma(\dot{\lambda})}{\Gamma(n + \dot{\lambda})} \prod_{i=1}^n \prod_{k=1}^K \frac{\Gamma(x_{i,k} + \lambda_k)}{\Gamma(\lambda_k)}, \\[6pt] &\propto \frac{\Gamma(\dot{\lambda})}{\Gamma(n + \dot{\lambda})} \prod_{k=1}^K \frac{1}{\Gamma(\lambda_k)} \prod_{i=1}^n \Gamma(x_{i,k} + \lambda_k), \\[6pt] \end{align}

and the prior kernel is:

\begin{align} p(\boldsymbol{\lambda}) &= \prod_{k=1}^K \text{Gamma}(\lambda_k | \alpha, \beta) \quad \quad \quad \quad \quad \quad \ \ \\[6pt] &\propto \prod_{k=1}^K \lambda_k^{\alpha-1} \exp(-\beta \lambda_k). \\[6pt] \end{align}

Consequently, the posterior kernel is:

\begin{align} p(\boldsymbol{\lambda} | \mathbf{x}_1,...,\mathbf{x}_1) &\propto L_\mathbf{x}(\boldsymbol{\lambda}) \cdot p(\boldsymbol{\lambda}) \\[6pt] &\propto \bigg[ \frac{\Gamma(\dot{\lambda})}{\Gamma(n + \dot{\lambda})} \prod_{k=1}^K \frac{1}{\Gamma(\lambda_k)} \prod_{i=1}^n \Gamma(x_{i,k} + \lambda_k) \bigg] \bigg[ \prod_{k=1}^K \lambda_k^{\alpha-1} \exp(-\beta \lambda_k) \bigg] \\[6pt] &\propto \frac{\Gamma(\dot{\lambda})}{\Gamma(n + \dot{\lambda})} \prod_{k=1}^K \frac{\lambda_k^{\alpha-1} \exp(-\beta \lambda_k)}{\Gamma(\lambda_k)} \prod_{i=1}^n \Gamma(x_{i,k} + \lambda_k). \\[6pt] \end{align}

This is a nasty-looking posterior kernel, and the corresponding scaling value leading to the posterior density does not have a closed form. There are a number of ways you could estimate the scaling value to get the posterior density (e.g., importance sampling, numerical integration methods, etc.), or create algorithms to sample directly from the posterior (e.g., MCMC methods).