Compound Poisson Distribution with sum of exponential random variables I'm trying to find the distribution and parameters in a Compound Poisson 
$S=\displaystyle\sum_{j=1}^{N}Y_{j},$ where $Y_{j}$ are exponential random variables independent and distributed identically with parameter $\alpha$ and $N$ is Poisson random variable with parameter $\lambda.$
I'm stuck because my attempts are based in computing the density probability function and I find a hard series: $$e^{-\lambda}e^{-\alpha x}\displaystyle\sum_{n=1}^{\infty}\frac{(\alpha\lambda)^{n}}{n!}\frac{x^{n-1}}{(n-1)!}.$$
The expression above was obteined using theorem of conditional probability in the random variable $N.$
If I use generating function I get $M_{S}(t)=e^{\lambda((\frac{\alpha}{\alpha-t})-1)},$ but I don't know the distribution of $S.$ This expression was got using conditional expectation again in random variable $N.$
Any kind of help is thanked in advanced.
 A: The Question
Let  $X \sim \text{Exponential}(\alpha)$, and let $\{X_1, X_2,\dots, X_n\}$ denote an iid sample of size $n$, where the sample size $n$ (instead of being fixed) is itself a Poisson random variable $N=n$. The OP seeks the distribution of the sample sum:
$$S = X_1 + X_2 +  \dots + X_n \quad \quad \text{where} \quad N \sim \text{Poisson}(\lambda)$$
As $N$ is Poisson, and the domain of support of a Poisson includes 0, it follows that the sample size $N$ can be 0, in which case $S = 0$. This is important, because it means that $P(S = 0)$ will have discrete mass. 
Solution
To proceed, first note that the sum of $n$ independent identical $\text{Exponential}(\alpha)$ variables has a $\text{Gamma}(n,\alpha)$ distribution i.e. $S$ has pdf $f(s \; \big| \; N = n)$:

where parameter $N \sim \text{Poisson}(\lambda)$ with pmf $g(n)$:

We seek the parameter mixture distribution of $S$ and $N$.
Unconditional pdf of $S$


*

*Discrete Part:  $S = 0$ iff $n = 0$. This occurs with probability $P(N=0)$:





*

*Continuous Part: 
The parameter-mix distribution, for $S>0$, is given by $\mathbb{E}_g[f]$:



where:


*

*I am using the Expect function from the mathStatica package for Mathematica

*Hypergeometric0F1Regularized denotes the confluent hypergeometric function
In summary, the unconditional pdf of $S$ is:
$$\text{pdf}(S) = \left\{
\begin{array}{cc}
e^{-\lambda} & \text{ if } s = 0 \\
\text{sol} & \text{ if } s > 0 \\
\end{array}\right.$$
which is a mixed discrete-continuous distribution.
A: Let $\{Y_i\}_{i\geq 1}$ be a sequence of IID $\mathrm{Exp}(\alpha)$ random variables, and let 
$$
  S_N=\sum_{i=1}^N Y_i,
$$
in which $N\sim\mathrm{Poisson}(\lambda)$. We know that $S_N\mid N=n\sim\mathrm{Gamma}(n,\alpha)$, for $n
\geq 1$. Hence,
\begin{align*}
  \Pr\{S_N\in B\} &= \mathrm{E}\left[\Pr\{S_N\in B\mid N\}\right] \\ 
&=e^{-\lambda}I_B(0)+\sum_{n=1}^\infty \left(\frac{e^{-\lambda}\lambda^n}{n!} \int_B \frac{\alpha^n}{(n-1)!}\,u^{n-1}e^{-\alpha u}\,I_{(0,\infty)}(u)\,du\right).
\end{align*}
TOL <- 0.01
lambda <- 4
alpha <- 2

B <- c(0, 3.5)

prob <- 0
if (B[1] == 0) prob <- exp(-lambda)

n <- 1
repeat {
  next_term <- exp(-lambda+n*log(lambda)-lfactorial(n)) * 
                 (pgamma(B[2], shape = n, rate = alpha) - 
                  pgamma(B[1], shape = n, rate = alpha))
  if (next_term < TOL) break
  prob <- prob + next_term
  n <- n + 1
}

print(prob)

A: The cumulative distribution does not have a simple closed form
expression, nor does the density.  Note that there is an atom at $S =
0$ with mass $\mathrm{Pr}\{N = 0\} = e^{-\lambda}$, so the density is
for $S \vert S > 0$.
The series in the density can be related to the Bessel functions
$I_0(x)$ and $I_1(x)$. But since this is a special case of the compound
  Poisson-Gamma distribution which itself is a special case of the
Tweedie distribution, usable computing tools can be found
under this name.
EDIT To derive an expression of the density, consider the following series 
$$
R(y) := \sum_{n=0}^{\infty} \frac{1}{n!\,n!} \, y^{n},  \qquad
r(y) := \sum_{n=1}^{\infty} \frac{1}{n!\,(n-1)!} \, y^{n-1} = R'(y). 
$$
Note that $R(y) = I_0(2 \sqrt{y})$ where $I_0(z)$ is the usual
modified Bessel function,
with derivative $I_0'(z) = I_1(z)$. So, using the expression given in the
question and some simple algebra we get the density: $f(x) = p
\,\delta(x) + (1 - p) \, f_1(x)$ where $\delta(x)$ (abusively) stands
for a Dirac density, $p:= e^{-\lambda}$ and
$$
   f_1(x) =  \frac{p}{1 -p} \, e^{ - \alpha x } \alpha \lambda \,  
   \frac{I_1(2 \sqrt{\alpha \lambda x})}{\sqrt{\alpha \lambda x}}
   \qquad \text{for } x > 0,
$$
which is the density of $S$ conditional on $S > 0$. This must be the
same solution as that of the answer by @wolfies, up to the $\alpha
\leftrightarrow 1/ \alpha$ change of notation therein.  The Bessel
functions $I_\nu(z)$ are widely available, e.g. in R using besselI.
