Generic sum of Gamma random variables I have read that the sum of Gamma random variables with the same scale parameter is another Gamma random variable. I've also seen the paper by Moschopoulos describing a method for the summation of a general set of Gamma random variables. I have tried implementing Moschopoulos' method but have yet to have success.
What does the summation of a general set of Gamma random variables look like? To make this question concrete, what does it look like for:
$\text{Gamma}(3,1) + \text{Gamma}(4,2) + \text{Gamma}(5,1)$
If the parameters above are not particularly revealing, please suggest others.
 A: According to Ansari et al. 2012, the PDF and CDF of independent gamma random variables with different distribution can be expressed in terms of Fox's Ĥ function (H-bar function).  The paper referenced below also contains an implementation of this function in the Wolfram language.
REFERENCES:

*

*Ansari, I.S., Yilmaz, F., et al., "New Results on the Sum of Gamma Random Variates With Application to the Performance of Wireless Communication Systems over Nakagami-m Fading Channels", arXiv:1202.2576 [cs.IT], 2012.

A: First, combine any sums having the same scale factor: a $\Gamma(n, \beta)$ plus a $\Gamma(m,\beta)$ variate form a $\Gamma(n+m,\beta)$ variate.
Next, observe that the characteristic function (cf) of $\Gamma(n, \beta)$ is $(1-i \beta  t)^{-n}$, whence the cf of a sum of these distributions is the product
$$\prod_{j} \frac{1}{(1-i \beta_j  t)^{n_j}}.$$
When the $n_j$ are all integral, this product expands as a partial fraction into a linear combination of $(1-i \beta_j  t)^{-\nu}$ where the $\nu$ are integers between $1$ and $n_j$.  In the example with $\beta_1 = 1, n_1=8$ (from the sum of $\Gamma(3,1)$ and $\Gamma(5,1)$) and $\beta_2 = 2, n_2=4$ we find
$$\begin{aligned}&\frac{1}{(1-i t)^{8}}\frac{1}{(1- 2i t)^{4}} = \\
&\frac{1}{(t+i)^8}-\frac{8 i}{(t+i)^7}-\frac{40}{(t+i)^6}+\frac{160 i}{(t+i)^5}+\frac{560}{(t+i)^4}-\frac{1792 i}{(t+i)^3}\\
&-\frac{5376}{(t+i)^2}+\frac{15360 i}{t+i}+\frac{256}{(2t+i)^4}+\frac{2048 i}{(2 t+i)^3}-\frac{9216}{(2t+i)^2}-\frac{30720 i}{2t+i}.
\end{aligned}$$
The inverse of taking the cf is the inverse Fourier Transform, which is linear: that means we may apply it term by term.  Each term is recognizable as a multiple of the cf of a Gamma distribution and so is readily inverted to yield the PDF.  In the example we obtain
$$\begin{aligned}
&\frac{e^{-t} t^7}{5040}+\frac{1}{90} e^{-t} t^6+\frac{1}{3} e^{-t} t^5+\frac{20}{3} e^{-t} t^4+\frac{8}{3} e^{-\frac{t}{2}} t^3+\frac{280}{3} e^{-t} t^3\\
&-128 e^{-\frac{t}{2}} t^2+896 e^{-t} t^2+2304 e^{-\frac{t}{2}} t+5376 e^{-t} t-15360 e^{-\frac{t}{2}}+15360 e^{-t}
\end{aligned}$$
for the PDF of the sum.
This is a finite mixture of Gamma distributions having scale factors equal to those within the sum and shape factors less than or equal to those within the sum.  Except in special cases (where some cancellation might occur), the number of terms is given by the total shape parameter $n_1 + n_2 + \cdots$ (assuming all the $n_j$ are different).

As a test, here is a histogram of $10^4$ results obtained by adding independent draws from the $\Gamma(8,1)$ and $\Gamma(4,2)$ distributions.  On it is superimposed the graph of $10^4$ times the preceding function.  The fit is very good.


Moschopoulos carries this idea one step further by expanding the cf of the sum into an infinite series of Gamma characteristic functions whenever one or more of the $n_i$ is non-integral, and then terminates the infinite series at a point where it is reasonably well approximated.
A: The Welch–Satterthwaite equation could be used to give an approximate answer in the form of a gamma distribution. This has the nice property of letting us treat gamma distributions as being (approximately) closed under addition. This is the approximation in the commonly used Welch's t-test.
(The gamma distribution is can be viewed as a scaled chi-square distribution, and allowing non-integer shape parameter.)
I've adapted the approximation to the $k, \theta$ parametrization of the gamma distriubtion:
$$
k_{sum} = { (\sum_i \theta_i k_i)^2 \over \sum_i \theta_i^2 k_i }
$$
$$
\theta_{sum} = { { \sum \theta_i k_i } \over k_{sum} }
$$
Let $k=(3,4,5)$, $\theta=(1,2,1)$
So we get approximately Gamma(10.666... ,1.5)
We see the shape parameter $k$ has been more or less totalled, but slightly less because the input scale parameters $\theta_i$ differ. $\theta$ is such that the sum has the correct mean value.
A: I will show another possible solution, that is quite widely applicable, and with todays R software, quite easy to implement.  That is the saddlepoint density approximation, which ought to be wider known!
For terminology about the gamma distribution, I will follow https://en.wikipedia.org/wiki/Gamma_distribution  with the shape/scale parametrization, $k$ is shape parameter and $\theta$ is scale.  For the saddlepoint approximation I will follow  Ronald W Butler: "Saddlepoint approximations with applications" (Cambridge UP).  The saddlepoint approximation is explained here: How does saddlepoint approximation work?
here I will show how it is used in this application.
Let $X$ be a random variable with existing momentgenerating function
$$
  M(s) = E e^{sX}
$$ which must exist for $s$ in some open interval that contains zero.  Then define the cumulant generating function by
$$
   K(s) = \log M(s)
$$
It is known that $E X = K'(0), \text{Var} (X) = K''(0)$.  The saddlepoint equation is $$ K'(\hat{s}) = x$$ which implicitly defines $s$ as a function of $x$ (which must be in the range of $X$). We write this implicitely defined function as
$\hat{s}(x)$. Note that the saddlepoint equation always has exactly one solution, because the cumulant function is convex.
Then the saddlepoint approximation to the density $f$ of $X$ is given by
$$
   \hat{f}(x) = \frac1{\sqrt{2\pi K''(\hat{s})}} \exp(K(\hat{s}) - \hat{s} x)
$$
This approximate density function is not guaranteed to integrate to 1, so is the unnormalized saddlepoint approximation.  We could integrate it numerically and the renormalize to get a better approximation.  But this approximation is guaranteed to be non-negative.
Now let $X_1, X_2, \dots, X_n$ be independent gamma random variables, where $X_i$ has the distribution with parameters $(k_i, \theta_i)$. Then the cumulant generating function is
$$
  K(s) = -\sum_{i=1}^n k_i \ln(1-\theta_i s)
$$ defined for $s<1/\max(\theta_1, \theta_2, \dots, \theta_n)$.
The first derivative is
$$
  K'(s) = \sum_{i=1}^n \frac{k_i \theta_i}{1-\theta_i s}
$$
and the second derivative is
$$
  K''(s) =   \sum_{i=1}^n \frac{k_i \theta_i^2}{(1-\theta_i s)^2}.   
$$
In the following I will give some R code calculating this, and will use the parameter values $n=3$, $k=(1,2,3)$, $\theta=(1,2,3)$.  Note that the following R code uses a new argument in the uniroot function introduced in R 3.1, so will not run in older R's.
shape <- 1:3 #ki
scale <- 1:3 # thetai
    # For this case,  we get expectation=14,  variance=36
make_cumgenfun  <-  function(shape, scale) {
    # we return list(shape, scale, K, K', K'')
    n  <-  length(shape)
    m <-   length(scale)
    stopifnot( n == m, shape > 0, scale > 0 )
    return( list( shape=shape,  scale=scale, 
          Vectorize(function(s) {-sum(shape * log(1-scale * s) ) }),
          Vectorize(function(s) {sum((shape*scale)/(1-s*scale))}) ,
          Vectorize(function(s) { sum(shape*scale*scale/(1-s*scale))
      }))    )
    }
    
solve_speq  <-  function(x, cumgenfun) {
       # Returns saddle point!
       shape <- cumgenfun[[1]]
       scale <- cumgenfun[[2]]
       Kd  <-   cumgenfun[[4]]
       uniroot(function(s) Kd(s)-x,lower=-100,
                      upper = 0.3333, 
                      extendInt = "upX")$root
    }
    
make_fhat <-  function(shape,  scale) {
    cgf1  <-  make_cumgenfun(shape, scale)
    K  <-  cgf1[[3]]
    Kd <-  cgf1[[4]]
    Kdd <- cgf1[[5]]
    # Function finding fhat for one specific x:
    fhat0  <- function(x) {
    # Solve saddlepoint equation:
    s  <-  solve_speq(x, cgf1)
    # Calculating saddlepoint density value:
    (1/sqrt(2*pi*Kdd(s)))*exp(K(s)-s*x)
        }
        # Returning a vectorized version:
   return(Vectorize(fhat0))
} #end make_fhat
    
fhat  <-  make_fhat(shape, scale)
plot(fhat, from=0.01,  to=40, col="red", 
         main="unnormalized saddlepoint approximation\n
               to sum of three gamma variables")

resulting in the following plot:

I will leave the normalized saddlepoint approximation as an exercise.
A: An exact solution to the convolution (i.e., sum) of $n$ gamma distributions is given as Eq. (1) in the linked pdf by DiSalvo. As this is a bit long, it will take some time to copy it over here. For only two gamma distributions, their exact sum in closed form is specified by Eq. (2) of DiSalvo and without weights by Eq. (5) of Wesolowski et al., which also appears on the CV site as an answer to that question. That is, $$\mathrm{G}\mathrm{D}\mathrm{C}\left(\mathrm{a}\kern0.1em ,\mathrm{b}\kern0.1em ,\alpha, \beta; \tau \right)=\left\{\begin{array}{cc}\hfill \frac{{\mathrm{b}}^{\mathrm{a}}{\beta}^{\alpha }}{\Gamma \left(\mathrm{a}+\alpha \right)}{e}^{-\mathrm{b}\tau }{\tau^{\mathrm{a}+\alpha-1}}{}_1F_1\left[\alpha, \mathrm{a}+\alpha, \left(\mathrm{b}-\beta \right)\tau \right],\hfill & \hfill  \tau >0\hfill \\ {}\hfill \kern2em 0\kern6.6em ,\hfill  \kern5.4em \tau \kern0.30em \le \kern0.30em 0\hfill \end{array}\right.,$$
where the notation in the questions above; $Gamma(a,b) \rightarrow \Gamma(a,1/b)$, here. That is, $b$ and $\beta$ are rate constants here and not time scalars.
