Solution (source code)

= Solution

Use the shape-rate convention for the <gamma distribution>. A shape $r$ and rate $\beta$ give <expected value> $r/\beta$ and <variance> $r/\beta^2$. The two moments here force $r=2$ and $\beta=2/\mu$. Thus $M_X(t)=\beta^2/(\beta-t)^2$, and the <geometric-sum moment-generating function> becomes
$$
M_S(t)=\frac{p\beta^2}{(\beta-t)^2-(1-p)\beta^2}.
$$
Put $r_0=\sqrt{1-p}$, $a=\beta(1-r_0)$ and $b=\beta(1+r_0)$. Then $0<a<b$, $ab=p\beta^2$, and
$$
M_S(t)=\frac{a}{a-t}\frac{b}{b-t}\qquad(t<a).
$$
Consequently the <geometric sum of shape-two gamma variables> has the same <probability distribution> as the sum of two independent <exponential distributions> with rates $a,b$. Their <convolution of independent random variables> yields
$$
\boxed{f_S(s)=\frac{ab}{b-a}(e^{-as}-e^{-bs})
=\frac{p}{\mu\sqrt{1-p}}
\left(e^{-2(1-\sqrt{1-p})s/\mu}-e^{-2(1+\sqrt{1-p})s/\mu}\right),\quad s>0.}
$$
The <probability density function> is zero for $s<0$. It is nonnegative because $a<b$, and its integral is
$$
\frac{ab}{b-a}\left(\frac1a-\frac1b\right)=1.
$$
This verifies normalization directly. An alternative check uses the conditional <gamma distribution> with shape $2n$:
$$
\sum_{n\geq1}p(1-p)^{n-1}\frac{\beta^{2n}s^{2n-1}e^{-\beta s}}{(2n-1)!}
=\frac{p\beta}{r_0}e^{-\beta s}\sinh(\beta r_0s),
$$
which is exactly the same density. The <expected value> $1/a+1/b=\mu/p$ supplies a further check.