# Compound Binomial-Exponential: Closed form for the PDF?

Setup:

Consider the random variable $$Y_N$$ derived from $$Y_N = \sum_{i=0}^N X_i$$ where $$N$$ is a random variable with distribution $$p_n = {M \choose n} p^n q^{M-n}$$ (binomial i.e. $$M$$ Bernoulli trials) and the $$X_i$$ are random variables with distribution $$f(x) = \lambda e^{-\lambda x}$$ (exponential). I'd like to get the pdf of $$y$$. Previously I asked this question which led me to expect the distribution is the sum $$p(y) = \sum_{i=1}^M {M \choose i} p^i q^{M-i} \frac{\lambda^iy^{i-1}e^{-\lambda y}}{(i-1)!}$$ or the integral (inverse Fourier transform of the characteristic function) $$p(y) = \int_{-\infty}^{\infty} \frac{dk}{2\pi} e^{-iky}(\frac{p\lambda}{\lambda-ik}+q)^M = \int_{-\infty}^{\infty} \frac{dk}{2\pi}e^{-iky}(\frac{i\lambda + kq}{i\lambda + k})^M$$ The second equality follows from $$p+q=1$$. I'd like to solve either one of these problems in closed form if possible. The moment generating function these derive from is $$p(s) = (\frac{p\lambda}{\lambda-s}+q)^M.$$ The characteristic function is $$p(s=ik)$$.

Question:

What is the distribution $$p(y)$$ of $$Y_N$$?

Attempts:

I tried binomial expanding $$p(s)$$ in the fourier transform form, which leads to a sum over integrals like $$\int_{-\infty}^{\infty} du u^{-m}e^{-a u}$$ where $$m$$ is a positive integer -- almost a sum of gamma functions but with negative parameters and over the wrong limits.

I tried manipulaing the sum form into something I could use the binomial formula on, but I can't obtain binomial coefficients from the factorial prefactors in the sum.

Any thoughts? I had thought a compound binomial-exponential process would have been well studied with a closed form solution but I couldn't find anything immediately relevant.

Other:

This paper explores a very similar problem but I'll have to read carefully to see how it relates (link).

• Your definition of the random number is a little confusing is it N or (as it appears in usage) n or i? Dec 14, 2018 at 22:24
• $M$ is a parameter of the binomial distribution (it was $N$ before), i.e. the number of Bernoulli trials of probabilities $p$ and $q$ which contribute to the distribution. $N$ is a random variable, and $n$ is a particular realization of that random variable. Edited thanks @MikeEarnest Dec 14, 2018 at 22:31
• The question has no answer without any independence assumptions. Dec 14, 2018 at 23:28
• The $X_i$ are all independent and identically distributed. $X_i$ and $N$ are independent. These assumptions are embedded into the two expressions I gave for $p(y)$. Dec 14, 2018 at 23:35
• Your first expression for $p(y)$ can be written $e^{-\lambda y}\frac{q^M}y\sum_{i=0}^M \binom{M}i \frac{t^i}{(i-1)!}$, where $t=(py\lambda/q)$. I am pretty sure this has no closed form. Dec 16, 2018 at 19:24

I was able to compute the integral you presented using contour integration in the complex plane.

For the case of $$\lambda > 0$$, $$y <0$$, it is easy to show via contour integration that

$$\int_{-\infty}^\infty \dfrac{dk}{2\pi}e^{-iyk}\left(\dfrac{i\lambda p}{k - (-i\lambda)}+q\right)^M = 0$$

For the more interesting case of $$\lambda > 0$$, $$y > 0$$, one can show via contour integration, and algebraic manipulations to generate a Laurent Series, that

\begin{align*}\int_{-\infty}^\infty \dfrac{dk}{2\pi}e^{-iyk}\left(\dfrac{i\lambda p}{k - (-i\lambda)}+q\right)^M &= -2\pi i \space \mathrm{Res}_{z=-i\lambda}\left[\dfrac{1}{2\pi}e^{-iyz}\left(\dfrac{i\lambda p}{z - (-i\lambda)}+q\right)^M\right]\\ \\ &= -2\pi i \space \mathrm{Res}_{z=-i\lambda}\left[\dfrac{1}{2\pi}e^{-iyz}\sum_{j=0}^{M}\binom{M}{j}\left(\dfrac{i\lambda p}{z - (-i\lambda)}\right)^jq^{M-j}\right]\\ \\ &= -2\pi i \space \mathrm{Res}_{z=-i\lambda}\left[\dfrac{1}{2\pi}\left(\sum_{n=0}^\infty \dfrac{\left(-iy\right)^n e^{-\lambda y}}{n!}\left(z-\left(-i\lambda\right)\right)^n\right)\left(\sum_{j=0}^{M}\binom{M}{j}\left(\dfrac{i\lambda p}{z - (-i\lambda)}\right)^jq^{M-j}\right)\right]\\ \\ &= -2\pi i \dfrac{1}{2\pi}\sum_{j=1}^M \dfrac{\left(-iy\right)^{j-1}e^{-\lambda y}}{\left(j-1\right)!}\binom{M}{j}\left(i\lambda p\right)^jq^{M-j} \\ \\ &= \sum_{j=1}^M \binom{M}{j} p^j q^{M-j}\dfrac{\lambda^j y^{j-1} e^{-\lambda y}}{\left(j-1\right)!}\\ \\ &\mbox{(confirmation of exactly what you already thought to be true)}\\ \\ &= \dfrac{q^M e^{-\lambda y}}{y}\sum_{j=1}^M \binom{M}{j} \left(\dfrac{\lambda y p}{q}\right)^j \dfrac{1}{\left(j-1\right)!}\\ \end{align*}

For the case of $$\lambda > 0$$, $$y =0$$, the integral diverges due to the "$$+q$$". The integral diverges at this single value of $$y$$ with a weight of $$q^M$$, given some hindsight with the convolution theorem and the definition of a CDF.

The integral really does not provide any new information regarding the series. You need to manipulate the series form into something you find palatable.

The above development of the Laurent Series, to find the residue, does provide some insights into how portions of each term in the final sum arise. The binomial coefficients come from the binomial expansion, as you would expect. The factorial factors come from the Taylor Series expansion of the exponential function. The final summation comes from multiplying series terms together and adding to find the coefficient of the $$(z-(-i\lambda))^{-1}$$ term of the Laurent Series.

Consolidating the results for all real values of $$y$$, the expression should really be written as:

\begin{align*}p(y) &= q^M\delta(y)+H(y)\dfrac{q^M e^{-\lambda y}}{y}\sum_{j=1}^M \binom{M}{j} \left(\dfrac{\lambda y p}{q}\right)^j \dfrac{1}{\left(j-1\right)!}\\ \\ &= q^M\delta(y)+ H(y) \dfrac{p}{q} q^{M} \lambda e^{-\lambda y}\sum_{j=0}^{M-1} \binom{M}{j+1} \left(\dfrac{p}{q}\lambda y\right)^{j} \dfrac{1}{j!}\\ \end{align*}

Where $$H(y)$$ is the Heaviside unit step and $$\delta(y)$$ is the Dirac Delta function.

Using this expression, convolving the $$M=1$$ case with itself $$M-1$$ times, yields the correct expression for the $$M^{th}$$ case, as expected given the convolution theorem of the Fourier Transform.

Using the above expression for $$p(y)$$, we find

$$\int_{-\infty}^\infty p(y) dy = (q+p)^M =1$$

which is what one would expect.

To replace the series, Wolfram Alpha comes up with this "closed" form in terms of a hypergeometric function:

\begin{align*}p(y) &= q^M\delta(y)+H(y) \dfrac{p}{q} q^{M} \lambda e^{-\lambda y}\sum_{j=0}^{M-1} \binom{M}{j+1} \left(\dfrac{p}{q}\lambda y \right)^{j} \dfrac{1}{j!}\\ \\ &= q^M\delta(y)+ H(y) \dfrac{p}{q} q^{M} \lambda e^{-\lambda y} M \space {}_1F_1 \left(1-M;2;-\dfrac{p}{q}\lambda y\right) \end{align*}

Update to derive the CDF from the PDF

To derive the CDF, let's rewrite the PDF a little bit. I'll be using the following expression involving the Incomplete Upper Gamma Function in the process

$$\dfrac{\partial}{\partial y}\Gamma(j,\lambda y) =\dfrac{\partial}{\partial (\lambda y)}\Gamma(j,\lambda y) \cdot \dfrac{\partial}{\partial y}\lambda y = -(\lambda y)^{j-1}e^{-\lambda y}\lambda$$

So working the PDF a bit

\begin{align*}p(y) &= q^M\delta(y)+H(y) \dfrac{p}{q} q^{M} \lambda e^{-\lambda y}\sum_{j=0}^{M-1} \binom{M}{j+1} \left(\dfrac{p}{q}\lambda y \right)^{j} \dfrac{1}{j!}\\ \\ &= \delta(y)q^M+H(y) \lambda e^{-\lambda y}\sum_{j=0}^{M-1} \binom{M}{j+1} q^{M-(j+1)}p^{j+1}\left(\lambda y \right)^{j} \dfrac{1}{j!}\\ \\ &= \delta(y)q^M+H(y) \lambda e^{-\lambda y}\sum_{j=1}^{M} \binom{M}{j} q^{M-j}p^{j}\dfrac{1}{(j-1)!}\left(\lambda y \right)^{j-1} \\ \\ &= \delta(y)q^M-H(y) \sum_{j=1}^{M} \binom{M}{j} q^{M-j}p^{j}\dfrac{1}{(j-1)!}(-1)\left(\lambda y \right)^{j-1} e^{-\lambda y}\lambda \\ \\ &= \delta(y)q^M-H(y) \sum_{j=1}^{M} \binom{M}{j} q^{M-j}p^{j}\dfrac{\partial}{\partial y}\dfrac{\Gamma(j,\lambda y)}{\Gamma(j)} \\ \end{align*}

To integrate the PDF to get the CDF, I'll be using the following relations:

$$\int_{-\infty}^{x} H(t)f'(t) dt = H(x)\int_0^x f'(t) dt = H(x)\left[f(x) - f(0)\right]$$

$$\int_{-\infty}^x \delta(t) f(t)dt = H(x)f(0)$$

So, integrating the PDF to find the CDF:

\begin{align*}P(y) &= \int_{-\infty}^y p(t) dt\\ \\ &= \int_{-\infty}^y \delta(t)q^M dt - \sum_{j=1}^{M} \binom{M}{j} q^{M-j}p^{j}\int_{-\infty}^y H(t)\dfrac{\partial}{\partial t}\dfrac{\Gamma(j,\lambda t)}{\Gamma(j)} dt\\ \\ &= H(y)q^M -\sum_{j=1}^{M} \binom{M}{j} q^{M-j}p^{j} H(y)\left[\dfrac{\Gamma(j,\lambda y)}{\Gamma(j)}-\dfrac{\Gamma(j,0)}{\Gamma(j)}\right]\\ \\ &= H(y)\left[q^M -\sum_{j=1}^{M} \binom{M}{j} q^{M-j}p^{j} \left(\dfrac{\Gamma(j,\lambda y)}{\Gamma(j)}-1\right)\right]\\ \\ &= H(y)\left[q^M +\sum_{j=1}^{M} \binom{M}{j} q^{M-j}p^{j} -\sum_{j=1}^{M} \binom{M}{j} q^{M-j}p^{j}\dfrac{\Gamma(j,\lambda y)}{\Gamma(j)}\right]\\ \\ &= H(y)\left[(q+p)^M -\sum_{j=1}^{M} \binom{M}{j} q^{M-j}p^{j}\dfrac{\Gamma(j,\lambda y)}{\Gamma(j)}\right]\\ \\ P(y) &= H(y)\left[1 -\sum_{j=1}^{M} \binom{M}{j} q^{M-j}p^{j}\dfrac{\Gamma(j,\lambda y)}{\Gamma(j)}\right]\\ \end{align*}

• Solid @Andy. I especially appreciate the contour integration to show my two expressions are equivalent. That's satisfying. I guess the hypergeometric solution is as good as it gets. I was hoping for something simpler for the application which brought this problem up, but it is what it is. Dec 22, 2018 at 2:36
• Thanks. I just added an update to rework the PDF a little, so I could also derive the CDF. Dec 23, 2018 at 17:55
• Yeah, there seems to be no easy way to simplify those $\Gamma(,)$ related weights from the otherwise binomial expansion. Dec 23, 2018 at 18:01