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.

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.
