Numerically/approximately integrating over independent gamma variables

Question

Problem Statement

For a problem in biology, I am testing out a joint distribution of the form:

$$ X \sim Multinomial(\frac{\theta_1}{\sum \theta_i}, ...,\frac{\theta_n}{\sum{\theta_i}}) \\ \theta_i \sim Gamma(\alpha_i, \beta_i) $$ (where I use the shape-rate parametrization of the Gamma distribution).

I am interested in the (logarithm of) marginal distribution $$ P(X | \alpha_{1..n}, \beta_{1..n}) = \int Multinomial(X|\Theta) \prod_i Gamma(\theta_i | \alpha_i, \beta_i) d \Theta $$

I know that in the special case of $\beta_1 = \beta_2 = ... = \bar\beta$, the marginal distribution is Dirichlet Multinomial (DM). The DM distribution however does not fit the data I am encountering, the data are more dispersed than what the DM distribution can accomodate. How could I evaluate/approximate the general integral at least slightly efficiently? I would need to this within a larger inference, so for many combinations of the parameters.

The $\alpha_i$ and $\beta_i$ come from a deeper level of the model and they are constrained so that $\sum E(log(\theta_i)) = 0$.

It is possible this is actually a very hard problem and pointer to why this is hard would be a sufficient answer.

Things I've tried

Letting $\beta_{1..n}$ to vary and treating $\theta_{1..n}$ as explicit latent variables, this model fits the data reasonably well, but is very slow even for small $n$ due to the large number of latent variables (I am using Stan).

I don't believe there is an analytic solution for $P(X | \alpha_{1..n}, \beta_{1..n})$ in the general case, but would it be possible to write an approximation or have a numerical integration scheme that would let me to compute the (log of) this density efficiently? It feels like the structure of independent Gammas could be somehow exploitable (in fact, the integral may be restated to be over i.i.d Gammas taking advantage of the Gamma's scaling properties).

I've tried a naive Monte-carlo scheme (sample from the gammas, compute the multinomial density, average over samples) which was hopelessly slow to converge even with few dimensions and a simple importance sampling scheme, using samples from the Dirichlet distribution as proposals for $\theta_1 / \sum \theta_i$ (to sample near the maximum likelihood area of the multinomial), but that only made things worse.

I've tried to find the Laplace approximation to $P(log(\Theta)|X, \alpha_i, \beta_i)$ which looks reasonably Gaussian-like for many parameter combinations, but turns out to be problematic, as when $x_i = 0$ and $\alpha_i \leq 1$, the mode is not defined. $P(\Theta)|X, \alpha_i, \beta_i)$ does not resemble Gaussian at all.

From approximation point of view the question What is the expected value of modified Dirichlet distribution? (integration problem) is related. But all my attempts ended as very poor approximations.

---

Slightly philosphical note: I believe DM does not fit my data because DM arises also from multinomial sampling of negative binomial variables with a fixed Fano factor (mean/variance ratio) - see e.g. Analytically solving sampling with or without replacement after Poisson/Negative binomial. Neg. binomial is Gamma-Poisson, where the $\beta$ parameter determines the Fano factor. But my data correspond to multinomial sampling of neg. binom variables where the Fano factor varies

Answer 1

Here is a sketch of a solution. The general strategy is to introduce the scalar sum $\theta_T=\sum\theta_i, $ then marginalize over it, then use it to turn the complicated multidimenional sum over the $\theta_{1..n}$ into a dynamic programming problem with aggressive cutoffs to prune solution space as the dynamic program runs.

We know that $\theta_T$ is a sum of gamma RVs with heterogeneous parameters. Since $\beta_i>0$ for all $i$, the variances are finite, and if $n$ is large enough, this will succumb to the CLT, so $\theta_T \rightarrow_d ~~ N(\mu_T,\sigma^2_T)$ where $\mu_T =\sum \alpha_i/\beta_i$ and $\sigma^2_T=\sum \alpha_i/\beta_i^2.$ This gives us a credible subset of $\theta_T$ over which to focus our attention on $\theta \in [\theta_{T_0},\theta_{T_1}]$. Even if $n$ isn't particularly large, we are able to approximate $p(\theta_T)$ through numerical methods and other methods, which are contingent on the data.

For each $\theta$ selected from this interval, we choose an integer $M$ so that $\theta/M$ is smallest quantity of $\theta$ to be allocated to $\theta_i$ in the partitioning of $\theta$ among the $\theta_{1..n}$. This turns the problem into a simple dynamic programming problem in which $\sum \theta_{1..n}=\theta_T$ is equivalent to $\sum m_i=M$ where $0\le m_i \le M, m_i \in \mathbb{N}. $ The number of solutions is $\binom{M+n-1}{n-1}$. For numbers like $M=20, n=10, \binom{29}{9} \sim 10^8$. Fortunately, we may implement the dynamic programming with an eye to selecting the solutions for which $P(\theta_i|\alpha_i,\beta_i)>\epsilon$ which will reduce the number tremendously. This is equivalent to $\sum m_i=M$ requiring that,e.g., $1 \le m_i \le M$, and now the number of solutions is $\binom{(M-n)+n-1}{n-1}$. For $M=20, n=10, \binom{19}{9}\sim 10^5,$ or 3 orders of magnitude smaller with a simple, obviously-correct requirement.

The cleverness enters when you memo-ize the dynamic program in order to reuse the solutions, then use a rule for varying $\alpha_{1..n},\beta_{1..n}$ perhaps to maintain constant values for $\mu_T$ and $\sigma^2_T$, and which will also maximize reuse of the dynamic programming solution. This constancy of $\mu_T$ and $\sigma^2_T$ may well be a within the realm of your needs.

Other tricks to making the problem manageable is what I would put in a more-than-just-a-sketch answer to this question.

In summary, $$ \begin{align} P(X|\alpha_{1..n},\beta_{1..n}) &= \sum_{\theta_T} P(X|\theta_T,\alpha_{1..n},\beta_{1..n})P(\theta_T| \alpha_{1..n},\beta_{1..n} )\\ &= \sum_{\theta_T} \sum_{\{\theta_{1..n}|\sum \theta_i=\theta_T\}} P(X|\alpha_{1..n},\beta_{1..n})N(\theta_T|\mu_T,\sigma^2_T). \end{align} $$

Then, $P(X|\alpha_{1..n},\beta_{1..n})= \binom{X}{x_1,x_2..x_n}\theta_1^{x_1} \cdots \theta_n^{x_1}$ need only be evaluated only for those $\theta_{1..n}$ values that are high probability in the first place. Note that the $\theta_{1..n}$ are implicitly functions of $\alpha_{1..n},\beta_{1..n}$.

Some observations

1. Finding a dynamic programming problem in the wild is always a lot of fun. It would make for a fairly high profile research result, I should think.

2. There are numbers $M$ and $n$ here where this approach becomes computationally very hard, even with dynamic programming. If $n$ has to be large, you may compensate by using smaller $M$ or larger $\epsilon.$ Ultimately, the problem may only be solvable if you are willing to a) use some restrictive priors, and/or b) accept a low resolution solution.

3. The choice to marginalize over $\theta_T$ means that your use of log-likelihood will have to be term-by-term in the sum, rather than all at once, which is not so inconvenient in the grand scheme of things.

Is this close to what you are looking for?

---

Addendum

A 20 line dynamic programming python script is able to enumerate $\binom{M+n-1}{n-1} M=100, ~n=20$ with a solution $4.910371215196107e+21$ in $0.18s$. This involves counting solutions by visiting them in the dynamic programming sort of way in which redundant solutions do not need to be revisited, just accounted for. In this enumeration, $\epsilon=0$. Visit-every-state enumeration is much slower, of course. For M=20 and n=10, the full enumeration that visits every one of the $10^8$ solutions currently takes 42s, which is where $\epsilon$ is valuable, and requirements that permit low resolution reporting of results.

---

The python code is included here as requested.

# -*- coding: utf-8 -*-
"""
Created on Fri Mar  8 21:38:57 2019

@author: Peter
"""
import numpy as np
import scipy.special
from datetime import datetime as dt

Solutions=[] # is a list of lists of stars grouped by bars

"""
This is the recursive routine that uses optimal substructure
-- the key design element of dynamic programming --
to *count* all of the distinct ways in which M identical stars
can be located in n pidgeon holes. This is an exact enumeration
of the solution to the stars-and-bars problem.

This algorithm can be adapted to *visit* explicitly instead of
leveraging dynamic programming and Memoization to merely *count* 
every solution, It is significantly slower if it visits the 
solutions. The advantage of visiting a solution is to eliminate 
solutions with a kpredictable low probability.

"""

"""
countMe 
"""
##########################################################
def countMe(m #The index of the pidgeon hole to start placing stars
            ,i #The number of stars that have to be placed
            ,V #
            ,visitEverySolution
            ):

    global Memo
    global Called;  Called += 1 # this counts the number of times CountMe was called
    if i<0: return # base case
    solns =0 # this is a local count of solutions for 
             # this instance, with M and i specified 
             # in fxn signature

    for x in reversed(range(0,m+1)): #counting down from m to 0
        VV=V # VV is mutable; V is not
        VV[i]=x
        for j in range(0,i):
            VV[j]=0 # clears the cubby holes to the left of i
                    # needed for recursion
        if sum(VV)==M: # we have placed all of the stars
            global Solns; Solns += 1
            solns += 1 # local solution
            global Solutions; Solutions.append(list(VV))
        else:
          if visitEverySolution or Memo[m-x,i-1]==0:
                countMe(m-x,i-1,VV,visitEverySolution)
          solns += Memo[m-x,i-1]
    Memo[m,i]=solns
    return solns                
##############################################################

"""
The beginning of execution is here.
"""

(M # number of stars
 ,n # number of bars+1 (or number of pidgeon holes)
 #)=(10,7) # 8008 solutions. use for testing
 )=(20,10) # 10015005 solutions
           # visitEverySolution=True takes 66 seconds
           # visitEverySolution=False takes 0.003 secodns



"""
Memo is a matrix containing the number of solutions with 
M stars up to and including the ith pidgeon hole.
Memo is the hallmark of dynamic programming.
"""
Memo=np.zeros((M+1,n)) 

# V this the running occupancy of the pidgeon holes
V=[0 for i in range(n)] #start with all pidgeon holes empty
Solns=0 # global number of solutions
Called=0 #global number of calls to countMe

print("An exact solution can computed to be",scipy.special.binom(M+n-1,n-1))

visitEverySolution=False #count only, not visit
t0=dt.now()
countMe(M,n-1,V,visitEverySolution) # enter the counting routine here
t1=dt.now()
# We want to locate M stars in all pidgeon holes up to 
# and including n-1 (the last one)

print("Elapsed time was",t1-t0)
print("The number of solutions with M=",M,"stars and n=",n,"pidgeon holes is",Memo[M,n-1])
print("Solns=",Solns)
print("CountMe was called",Called,"times.")
#print(Solutions)
```

Answer 2

Your question does not specify the trials parameter in the multinomial distribution, so I'm going to use the standard generic case where you have $n$ trial and $m$ categories. I will let $\tilde{\boldsymbol{\theta}} = (\tilde{\theta}_1,...,\tilde{\theta}_m)$ denote the normalised values $\tilde{\theta}_k = \theta_k / \sum_i \theta_i$, so the model of interest is:

$$\mathbf{X} | \boldsymbol{\theta} \sim \text{Mu}(n,\tilde{\boldsymbol{\theta}}) \quad \quad \quad \theta_k \sim \text{Ga}(\alpha_k, \beta_k).$$

---

The required integral: The marginal mass function for $\mathbf{X}$ is given by:

$$\begin{equation} \begin{aligned} p(\mathbf{x}|\boldsymbol{\alpha}, \boldsymbol{\beta}) &= \int \limits_\Theta \text{Mu}(\mathbf{x} | \tilde{\boldsymbol{\theta}}) \prod_{i=1}^m \text{Ga}(\theta_i | \alpha_i, \beta_i) \ d \boldsymbol{\theta} \\[6pt] &\overset{\mathbf{x}}{\propto} {n \choose \mathbf{x}} \int \limits_\Theta \frac{\prod_{i=1}^m \theta^{x_i + \alpha_i - 1} \exp( - \beta_i \theta_i )}{(\sum_{i=1}^m \theta_i)^n} \ d \boldsymbol{\theta}. \\[6pt] \end{aligned} \end{equation}$$

If $\beta_1 = ... = \beta_m = \beta$ then $\tilde{\boldsymbol{\theta}}$ has a Dirichlet distribution which means that $\mathbf{X}$ then has a multinomial-Dirichlet distribution. (In your question you incorrectly state that $\mathbf{X}$ has a Dirichlet distribution in this case. That can't be right because the Dirichlet distribution generates non-integer vectors of values.) There are various ways to approximate the above integral, and I will show you one of these methods below. (If time permits, I may come back to this answer and add other methods.)

---

Approximation by direct sampling (useful if hyperparameters are fixed): Since the object of interest is a mass function obtained via a known conditional distribution, the simplest way to approximate the integral is simply to generate a large number of conditioning values and approximate the integral by an analogous sum. To do this, choose some large number of simulations $H \in \mathbb{N}$ and generate independent values:

$$\theta_k^{(h)} \sim \text{Ga}(\alpha_k, \beta_k) \quad \quad \quad \text{for all } k=1,...,m \text{ and } h =1,...,H.$$

You can then use the approximation:

$$\hat{p}_\text{DS}(\mathbf{x}|\boldsymbol{\alpha}, \boldsymbol{\beta}) = \frac{1}{H} \sum_{h=1}^H {n \choose \mathbf{x}} \frac{\prod_{k=1}^m \theta_k^{(h) \ x_k}}{(\sum_{k=1}^m \theta_k^{(h)})^n} .$$

This method requires you to generate $H \times m$ scalar parameter values and then evaluate the mass function as an average of $H$ multinomial mass functions. Note that this method evaluates the mass function for a fixed set of parameter values $\boldsymbol{\alpha}$ and $\boldsymbol{\beta}$, so it is not useful if you would like your approximation to serve as a function of those parameter values.

---

Approximation by importance sampling (useful if hyperparameters are variables): The method of direct sampling builds the parameters $\boldsymbol{\alpha}$ and $\boldsymbol{\beta}$ into the simulation, so it is useful when these are fixed (i.e., when you don't need the function to vary over these values). If you would like your approximating function to also be able to vary over these parameters then you can use importance sampling instead. To do this, choose some mid-ranged hyperparameter values $\bar{\alpha}$ and $\bar{\beta}$, choose some large number of simulations $H \in \mathbb{N}$ and generate independent values:

$$\theta_k^{(h)} \sim \text{Ga}(\bar{\alpha}, \bar{\beta}) \quad \quad \quad \text{for all } k=1,...,m \text{ and } h =1,...,H.$$

You can then use the approximation:

$$\hat{p}_\text{IS}(\mathbf{x}|\boldsymbol{\alpha}, \boldsymbol{\beta}) = \frac{1}{H} \sum_{h=1}^H {n \choose \mathbf{x}} \frac{\prod_{k=1}^m \theta_k^{(h) \ x_k}}{(\sum_{k=1}^m \theta_k^{(h)})^n} \cdot \prod_{k=1}^m \frac{\text{Ga}(\theta_k^{(h)} | \alpha_k, \beta_k)}{\text{Ga}(\theta_k^{(h)} | \bar{\alpha}, \bar{\beta})}.$$

This method requires you to generate $H \times m$ scalar parameter values and then evaluate the mass function as an average of $H$ weighted multinomial mass functions, where the weighting is used to adjust for the hyperparameters. Note that this method treats the hyperparameters as variables that can be adjusted in the calculation of the approximating function (without generating new simulated values). Here is some R code to implement this latter method:

#Load required libraries and set seed
library(stats);
library(matrixStats);
set.seed(1);

#Set parameter values for simulation
m <- 3;
A <- 2;
B <- 3;
H <- 10^4;

#Simulate values of theta
THETA <- array(rgamma(H*m, shape = A, rate = B), dim = c(H, m));

#Create function to generate approximation using importance sampling 
LOGMASS  <- function(x, alpha, beta) { 

    LOGWEIGHT <- array(0, dim = c(H, m));
    for (h in 1:H) {
    for (k in 1:m) {
        T1 <- dgamma(THETA[h,k], shape = alpha[k], rate = beta[k], log = TRUE);
        T2 <- dgamma(THETA[h,k], shape = A,        rate = B,       log = TRUE);
        LOGWEIGHT[h,k] <- (T1 - T2); } }
    LOGW <- rowSums(LOGWEIGHT);

    LOGTERMS <- rep(0, H);
    for (h in 1:H) {
    LOGTERMS[h] <- (stats::dmultinom(x, size = sum(x), prob = THETA[h, ],
                                     log = TRUE) + LOGW[h]); }

    matrixStats::logSumExp(LOGTERMS) - log(H); } 

#Find approximate mass function at an example point
x     <- c(3,3,4);
alpha <- c(4,3,2);
beta  <- c(2,1,2);
LOGMASS(x, alpha, beta);

[1] -4.132541

This code will allow you to generate the log-mass for your distribution at any input value for the observable vector and the hyperparameters. By setting H reasonably large you should get a good approximation to the true log-mass function. Some statistical theory and accuracy bounds for this kind of simulation can be found in O'Neill (2009).

---