Bayesian updates for Dirichlet-multinomial with Gamma prior

Question

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.

Answer 1

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.

Answer 2

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