Distribution of sum of random variables Let $X_1, X_2, . . .$ be independent exponential random variables with mean $1/\mu$ and let $N$ be
a discrete random variable with
$P(N = k) = (1 − p)p^{k-1}$  for $k = 1, 2, . . . $
where $0 ≤ p < 1$ (i.e. $N$ is a shifted geometric random variable). Show that $S$ defined as $S =\sum_{n=1}^{N}X_n$ is again exponentially distributed with parameter $(1 − p)\mu$.
My approach:
$S = \sum_{n=1}^{N}\sum_{k=1}^{\infty}\frac{1}{\mu}e^{-\frac{n}{\mu}}(1-p)p^{k-1}$
How to solve this? Is it the right approach to solve this problem?
 A: Not quite, as $S$ will be a convolution of distribution functions. This is much easier to do with either characteristic or moment generating functions. Specifically, let
$M_X(t)=E[e^{tX}],$
be the moment generating function of $X$. For an exponential with parameter $\lambda$, it's not hard to show that:
$$M(t)=\frac{\lambda}{\lambda-t}, \ t<\lambda.$$
Then by independence:
\begin{align*}
M_S(t)=\sum_{N\geq 1} M_X(t)^NP(N)&=\sum_{N\geq 1}\left(\frac{\lambda}{\lambda-t}\right)^N(1-p)p^{N-1}\\
&=(1-p)\left(\frac{\lambda}{\lambda-t}\right)\frac{1}{1-p\lambda/(\lambda-t)}\\
&=(1-p)\frac{\lambda}{\lambda-t-p\lambda}\\
&=\frac{\nu}{\nu-t},
\end{align*}
where $\nu:=(1-p)\lambda$, which we recognize to be the MGF of an exponential with parameter $\nu$. 
A: I would calculate the characteristic function :
$$E[ e^{itS} ] = E\left[ \prod_{i=1^N} e^{itX_1} \right]$$
$$= \sum_{k=1}^\infty (1-p)p^k \prod_{i=1}^k \int_{\mathbb{R}^+} e^{itx} \mu e^{-\mu x} dx$$
$$= \sum_{k=1}^\infty (1-p)p^{k-1} \prod_{i=1}^k \frac{\mu}{\mu-it}$$
$$  =\frac{(1-p)}{p}\sum_{k=1}^\infty \left( \frac{ \mu p }{\mu - it} \right)^k$$
$$ = \frac{(1-p)}{p} \frac{ 1 }{ 1-  \frac{ \mu p }{\mu - it} } - (1-p)$$
$$= \frac{(1-p)}{p} \frac{\mu - it - (\mu - it -   \mu p)}{ \mu - it -   \mu p}$$
$$ = \frac{(1-p)\mu }{ (1-p) \mu - it}$$
And this is exactly the characteristic function pof an exponential variable with parameter $(1-p)\mu$
A: Your approach is in fact OK if you compute more carefully (although characteristic function method is always preferable).
Given $N = n$, $S$ has Gamma distribution $\Gamma(n, \mu)$ (see here), i.e., the conditional density of $S$ given $N = n$ is 
$$f_{S|N = n}(s) = \frac{\mu^n}{\Gamma(n)}s^{n - 1}e^{-\mu s}, \quad s > 0. $$
Therefore the marginal cdf of $S$ can be computed as
\begin{align*}
& P(S \leq x) = \sum_{n = 1}^\infty P(S \leq x|N = n)P(N = n) \\
= & \sum_{n = 1}^\infty\left(\int_0^x \frac{\mu^n}{\Gamma(n)}s^{n - 1}e^{-\mu s} ds\right) \times (1 - p)p^{n - 1} \\
= & \mu(1 - p) \int_0^x \left(\sum_{n = 1}^\infty\frac{(\mu p s)^{n - 1}}{(n - 1)!}\right)e^{-\mu s}ds \\
= & \mu(1 - p)\int_0^x e^{\mu p s}e^{-\mu s} ds \\
= & \mu(1 - p) \int_0^x e^{-\mu(1 - p)s} ds \\
= & 1 - e^{-\mu(1 - p)x},
\end{align*}
which shows that $S \sim \exp(\mu(1 - p))$. The interchange of $\int$ and $\sum$ follows from the Fubini's theorem.
