# Numerically/approximately integrating over independent gamma variables

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

• But doesn't your model too (like the Dirichlet-Multinomial) only accomodate overdispersion? If your data is underdispersed as you say, then it is not surprising that you problems with convergence? – Jarle Tufto Mar 6 '19 at 13:15
• Sorry, I was a bit unclear - the DM model is underdispersed with regards to the data, the data are overdispersed with regards to the DM model. Will clarify in the text. – Martin Modrák Mar 6 '19 at 13:18
• Are all of the $\{\alpha_{1..n}, \beta_{1..n}\}$ parameters known and fixed, or is your goal to optimize $log(P)$ w.r.t. these parameters? That is, will your optimal solution have to be efficient in modifying these parameters, or is it enough for the solution to be a simple black-box with those parameters as inputs? – Peter Leopold Mar 8 '19 at 15:56
• No, the $\alpha_i$ and $\beta_i$ are not fixed - will clarify in the question. – Martin Modrák Mar 8 '19 at 15:58
• What is a ballpark number for $n$ in your application? – Peter Leopold Mar 9 '19 at 0:11

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{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}

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

• Thanks for the answer - I see a lot of effort went into it, but note that the main sampling scheme you proposed is what I mention in the "Things I've tried" section (albeit in less detail) and it turns out it is very slow (I need roughly 1e6 samples to reliably converge when $n = 2$) . The importance sampling step is a neat addition, I will check if it can rescue the inefficiency of the main sampling) – Martin Modrák Mar 9 '19 at 14:18

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?

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.

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
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)
$$$$
`
• I like this approach - it is something completely different than what I would have come up with, which is great. I am not sure it will be practical for my case, but it's definitely inspiring. Could you share your Python sript as a Github gist or something, so that I can play with it? – Martin Modrák Mar 9 '19 at 16:21
• OK, that was the encouragement I needed to keep going with it. I will certainly share my end point, which will be done before the bounty expires. When I picked dynamic programming, I moved from R to python as R isn't a good recursion platform, in my experience. Also, this problem is only solvable if it focuses most of its probability narrowly on small subset of $\{\theta_i\}$ so $n_{effective}<<200$. In fact, if there is any way to break the problem into pieces, that would help tremendously. Is it truly multinomial over all $n$ states at the same time? Or can we do better than that? – Peter Leopold Mar 9 '19 at 16:34
• Unfortunately, I can't really break the problem into pieces, it is multinomial over everything (and if I'm happy with ignoring the multinomial, I can just use neg. binomial). Not sure it is really sensible for you to spend a lot of time with it - I think I might need to break the assumptions that make the approach competitive. Will have more time to think about this properly tomorrow and will let you know. Also, if you do develop further, note that you can use the Dirichlet-Multinomial special case to check your accuracy against analytic solution. – Martin Modrák Mar 9 '19 at 17:55
• So on second thought, your approach is still interesting, but I don't think I can work with it - at least not directly. Some problems I see: selecting $M$ is not easy as I somehow need to be sensitive to differences between small values while also handling quite large values. Further, the $\epsilon$ pruning is questionable, as in my tests, you sometimes have a huge region of $\Theta$ where each point has small density, but the overall contribution of this region to $P(X, \alpha, \beta)$ is important, beacause it is large.. :-( – Martin Modrák Mar 10 '19 at 15:13
• I see in your comment above that $E(x_i)=100$ and $var(x_i) \sim 100$, too, probably. Also, you have $n$ values of $x_i$ and you are trying to infer $2n$ values $\alpha_{1..n},\beta_{1..n}$ so you have to use strong regularization, probably in the form of a prior. What would that be? – Peter Leopold Mar 10 '19 at 15:26