Compound Poisson Distribution with sum of exponential random variables

Question

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.$

Answer 1

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.

Answer 2

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.

Answer 3

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)

Answer 4

This comes a bit late, but I've just spent some time on this question, so here's a summary of the way the reasoning goes, based on the previously proposed answers:

• the sum of $n$ i.i.d. exponentially distributed variables is distributed as $f_0 = \delta(x)$, $f_n = \frac{\alpha^n x^{n-1} e^{-\alpha x}}{(n-1)!}$ (in which case it is a gamma distribution).
• then the probability of there being $n$ contributions in the total sum is given by the Poisson distribution of $n$, $w_n = e^{-\lambda} \frac{\lambda^n}{n!}$.

So, making a mix of the original post and parts of the different previous answers, we can make up an expression for the total pdf as: $$ f(x) = e^{-\lambda} \left[ \delta(x) + \sum_{n=1}^\infty \frac{\lambda^n}{n!} \frac{\alpha^n x^{n-1} e^{-\alpha x}}{(n-1)!}\right], $$ or, shifting the indices around: $$ f(x) = e^{-\lambda} \left[ \delta(x) + e^{-\alpha x} \lambda \alpha \sum_{n=0}^\infty \frac{(\lambda \alpha x)^n}{(n+1)! n!}\right]. $$ Where I contribute is by referencing Eq. 9.6.10 of Abramowitz & Stegun, which gives an expression for the ``hard series'' mentioned by the OP in terms of modified Bessel functions as $$ I_\nu(z) = (z/2)^\nu \sum_{k=0}^\infty \frac{(z^2/4)^k}{k! \Gamma(\nu+k+1)}, $$ from which we finally get $$ f(x) = e^{-\lambda} \left[ \delta(x) + e^{-\alpha x} \sqrt{\frac{\alpha \lambda}{x}} \mathrm{I_1} (2 \sqrt{\alpha \lambda x})\right]. $$

Answer 5

Here is an alternative approach by considering the process as a 2d random walk with drift:

$$X(t) = \sum_{k=1}^t X_k \\ Y(t) = \sum_{k=1}^t Y_k$$

where each $X_k$ and $Y_k$ are exponentially distributed with rate $\lambda$. Below is a simulation for $\lambda = 15$.

The question is then similar to the distribution of the position of $Y(t-1)$ for the lowest $t$ where $X(t) > 1$. Or approximately similar to the distribution of the position of $Y(t)$ for the lowest $t$ where $X(t) > 1$.

Or in simple words: what was/is the position of the distribution of $Y(t)$ when $X(t)$ passes the barrier $x=1$?

If the rate $\lambda$ is high, which means that there are many small steps, then the above process becomes approximately a 2D Wiener process with a drift $\nu = \lambda^{-1}$ and diffusion $\sigma = \lambda^{-1}$.

We could consider the

• Time $t$ that the $x$ coordinate hits the barrier. This is inverse Gaussian distributed.

  $$f(t) = \frac{1}{\sqrt{2\pi s^{2} t^3}} \exp \left( \frac{-(1-st)^2}{2s^2 t} \right)$$ where $s = \lambda^{-1}$

• Position $y(t)$ conditional on the barrier hit time. This is Gaussian distributed. $$g(y;t) = \frac{1}{\sqrt{2\pi s^2 t}} \exp \left( \frac{-(y-s t)^2}{2 s^2 t}\right)$$

And the distribution of $y$ when the coordinate $x$ hits the barrier is distributed as the compound distribution

$$h(y) = \int_0^\infty g(y;t)f(t) dt$$

when we work the integral out then we get

$$\begin{array}{} h(y) &=& \int_0^\infty \frac{1}{\sqrt{2\pi s^2 t^3}} \exp \left( \frac{-(1-s t)^2}{2s^2 t} \right)\frac{1}{\sqrt{2\pi s^2 t}} \exp \left( \frac{-(y-s t)^2}{2 s^2 t}\right) dt\\ &=&\frac{1}{2\pi s^2 } \int_0^\infty \frac{1}{t^2} \exp \left( \frac{-(1-s t)^2-(y-s t)^2}{2 s^2 t} \right) dt\\ &=& \frac{\sqrt{2}}{\pi s} \exp\left(\frac{y+1}{s}\right)\sqrt{\frac{1}{y^2+1}}K_1\left(\sqrt{\frac{2y^2+2}{s^2}}\right) \end{array}$$

where $K_1$ is the first order modified Bessel function of the second kind.

The graph below, with an emperical density of the simulated values along with the derived function, shows that this approach gives a reasonable approximation.

Non central $\chi^2$-distribution.

The other answers, using an exact approach, ended up with an expression that uses a slightly different Bessel function. Interestingly, those expressions relate to a non-central chi-squared distribution, which has the density

$$f_Y(y;n,\lambda) = \sum_{i=1}^\infty \frac{e^{-\lambda/2} (\lambda/2)^i}{i!} f_{Z_{k+2i}}(y)$$

where $f_{Z_{k+2i}}(y)$ is the pdf for the chi-squared distribution with $k+2i$, which for even values is equal to the Erlang distribution with rate $\lambda = 1$.

The case with $k = 0$ is equivalent to the problem here. This problem has been discussed before. Siegel, A. F. (1979), "The noncentral chi-squared distribution with zero degrees of freedom and testing for uniformity"

This relationship with the chi-squared distribution I still have to investigate further. To be continued. Intuition behind occurence of non central chi squared distribution in conditional coordinates of a random walk