# Concentration inequality for sums of independent gamma random variables

I am dealing with the following problem:

Say $$X_1, \ldots, X_n$$ are independent Gamma random variables, each one having shape and rate parameters $$\alpha_i$$ and $$\beta_i$$, respectively. Let $$S_n = \sum_{i=1}^n X_i$$ and I want to obtain a concentration inequality for $$S_n$$, essentially an upper bound for $$\mathbb{P}( \lvert S_n - \mathbb{E}(S_n)\rvert > t)$$.

I am aware of concentration inequalities that exist out there (e.g. the Bernstein or the Chernoff bound or the Moment bound etc.) but I am a bit confused as $$S_n$$ is not bounded. The moment bound could be used here I guess, even if the sum of these Gamma random variables is not really a Gamma random variable (given that the rate parameters differ). Any suggestions on what I could use would be more than welcome!

• The Chernoff bound is applicable: is it insufficient somehow? Commented Jul 2 at 11:59
• It is not insufficient, I know it can be used here but I was wondering if there's a better alternative; essentially a tighter bound. According to the Wikipedia page for the Chernoff bound (although I reckon it is not the most trustworthy source): "For unbounded random variables the bound is nowhere tight, though it is asymptotically tight up to sub-exponential factors ("exponentially tight"). Individual moments can provide tighter bounds, at the cost of greater analytical complexity." - this is why I mentioned the moment bound as a better alternative. Is there something even tighter? Commented Jul 2 at 12:15
• So what's the goal then? To try to find some sort of best generic approximation of the sum? The Chernoff should still be a decent approximation. The "never tight" remark is accurate, but it is also kind of a high bar to reach. We are usually happy once we get the right order of magnitude, and I think that's the case here? Commented Jul 2 at 12:30
• My advice would be. 1. Figure out exactly what you want. ie concrete succes metrics for the result of your calculation 2. Check if the Chernoff bound is what you want by running experiments to compare the Chernoff bound to monte-carlo estimates of the CDF of your sums. 3. Report here once that is done. (PS: If you struggle with 1, give me the wider context of your question. Is this part of some larger research project or is this just an exercise you are doing for fun?) Commented Jul 2 at 12:33
• Just be careful that you do not lose too much time chasing after a mirage. You can always improve a bound (until you reach the true CDF, but you rarely do). In your case, you should be able to prove that there is exponential decay in the tail. If you need to go beyond Chernoff, check out if there are bounds on the tail probability of log-concave densities. This review should be a decent starting point projecteuclid.org/journals/statistics-surveys/volume-8/… Commented Jul 2 at 12:53

There are tractable integral representations of the PDF and CDF of sums of gamma variables with arbitrary parameters (see Equation 1.3 here and Equation 7 here).

$$f_S(x)=\frac{1}{\pi}\int_0^\infty\frac{\cos\big(\sum_{k=1}^n\alpha_k\arctan(\beta_kt)-xt\big)}{\prod_{k=1}^n(1+(t\beta_k)^2)^{\alpha_k/2}}dt$$

$$F_S(x)=\frac{1}{2}-\frac{1}{\pi}\int_0^\infty\frac{\sin\big(\sum_{k=1}^n\alpha_k\arctan(\beta_kt)-xt\big)}{t\prod_{k=1}^n(1+(t\beta_k)^2)^{\alpha_k/2}}dt$$

(Note that $$\beta_k$$ are the scale parameters.)

Here is an R implementation that I've used before:

library(cubature)

f <- function(z, x, alpha, beta) {
bz <- beta%*%z
sin(crossprod(alpha, atan(bz)) - x*z)/z/
exp(crossprod(alpha/2, log1p(bz^2)))
}

psumgamma <- function(x, alpha, tau) { # exact gamma-sum CDF
0.5 - hcubature(f, 0, Inf, tol = 1e-8, vectorInterface = TRUE,
x = x, alpha = alpha, beta = 1/tau)$integral/pi }  Example usage: (seed <- sample(.Machine$integer.max, 1))
#> [1] 1972923043
set.seed(seed)

n <- 1e3              # number of gamma variables
mu <- rexp(n, 1e4)    # means
tau <- rexp(n, 1e-2)  # rate parameters
alpha <- mu*tau       # shape parameters
ES <- sum(mu)         # expected value of the sum of the gamma variates
t <- sum(mu)/2        # bounding value

system.time(
p <- 1 - psumgamma(ES + t, alpha, tau) + psumgamma(ES - t, alpha, tau)
)
#>    user  system elapsed
#>    0.06    0.00    0.06


Compare the probability against an estimate from simulation:

system.time({
library(parallel)

rsumgamma <- function(N, alpha, tau) {
# function to generate random samples of sums of N gamma variables with
# shape parameter alpha and rate parameter tau
n <- max(length(alpha), length(tau))

if (mean(alpha) < 1) {
# initial samples in log space
x <- matrix(log(rgamma(N*n, alpha + 1)) - rexp(N*n, alpha) - log(tau), n)
m <- colMaxs(x, TRUE)
# convert out of log space and sum
colsums(exp(setop(x, "-", m, rowwise = TRUE)))*exp(m)
} else colSums(matrix(rgamma(n*N, alpha, tau), n))
}

nc <- detectCores() - 1L
cl <- makeCluster(nc)
invisible(clusterEvalQ(cl, {library(Rfast); library(collapse)}))
N <- 1e6 # number of gamma sums to sample
N1 <- N%/%nc
N2 <- N%%nc
S <- unlist(
parLapply(cl, rep(c(N1, N1 + 1), c(nc - N2, N2)), rsumgamma, alpha, tau)
)
stopCluster(cl)
psim <- mean(abs(S - ES) > t)
})
#>    user  system elapsed
#>    0.07    0.11   20.47

c("Exact" = p, "Simulated" = psim)
#>     Exact Simulated
#> 0.2200846 0.2202790