# Stopped-sum of Exponential random variables

Let $\xi_1, \xi_2, \ldots \xi_n, \ldots$ - independent random variables having exponential distribution $p_{\xi_i} (x) = \lambda e^{- \lambda x}, \; x \ge 0$ and $p_{\xi_i} (x) = 0, \; x < 0$. Let $\nu = \min \{n \ge 1 : \xi_n > 1\}$. Need to find the distribution function of a random variable $g = \xi_1 + \xi_2 + \ldots \xi_{\nu}$ that is, find the probability $\mathbb{P}(g < x) = \mathbb{P} (\xi_1 + \xi_2 + \ldots \xi_{\nu} < x)$.

I made the following calculations:

$\mathbb{P} (\xi_1 + \xi_2 + \ldots \xi_{\nu} < x) = \sum_{k = 1}^{\infty} \mathbb{P} (\xi_1 + \xi_2 + \ldots \xi_k < x, \nu = k) = \sum_{k = 1}^{\infty} \mathbb{P} (\xi_1 + \xi_2 + \ldots \xi_k < x, \xi_1 \le 1, \ldots \xi_{k-1} \le 1, \xi_k > 1)$.

The probability of the sum can be represented as integral:

$\mathbb{P} (\xi_1 + \xi_2 + \ldots \xi_k < x, \xi_1 \le 1, \ldots \xi_{k-1} \le 1, \xi_k > 1) = \int\limits_D \lambda^k e^{- \lambda u_1} e^{- \lambda u_2} \ldots e^{- \lambda u_k} {d}u_1 \ldots {d}u_k$, where $D = \{ u_1 + \ldots u_k < x, u_1 \le 1, \ldots u_{k-1} \le 1, u_k > 1\}$. I'm afraid that this integral cannot be calculated.

Is it somehow easier to find the distribution function $\mathbb{P} (g < x)$?

• Is seems like you are implicitly assuming that the $\xi_n$ are ordered
– user76844
Jan 7 '15 at 1:28
• I've deleted my answer because it wasn't very useful. However, I did find this paper which addresses this exact subject: ams.org/journals/tpms/2008-76-00/S0094-9000-08-00740-0/… Jan 8 '15 at 22:25
• @Math1000 Thanks for the interesting link. However, the paper you reference is actually surprisingly different because it considers the sum $\sum_{i=1}^{N-1} X_i$ where $N$ is a rv ... whereas the OP considers the same sum PLUS the final $x_n$ term. This makes a substantive difference to the cdf: (i) Excluding the $x_n$ term will result in a large number of sample sums = 0, because the sample terminates at $n=1$, which in turn means that the cdf of the sum will have a discrete mass/jump at 0: see Fig 1 of the paper. (ii) By contrast, the OP's cdf of the sum will be continuous (no jumps). Jan 14 '15 at 15:34
• Yeah, I noticed that. However, it is probably the closest reference you're going to find. I worked on the problem for several hours trying to find a solution before finding this paper and realizing the problem is a bit more involved than I thought :( Jan 14 '15 at 15:55

The fact that this involves "convolutions" and sums of i.i.d. random variables makes me think of trying to deduce the distribution from Moment generating functions. Using the independence of the $\xi_i$, we have (for $t<\lambda$),

$$\begin{eqnarray*} \mathbb{E}\left[e^{t\sum\limits_{n=1}^{\nu}\xi_{n}}\right]&{}={}&\sum\limits_{r=1}^{\infty}\int{\textbf{1}}_{\left\{\xi_{1}\leq 1,\,\ldots\,,\,\xi_{r-1}\leq 1\,,\,\xi_{r}>1 \right\}}e^{t\left(\sum\limits_{n=1}^{r}z_{n}\right)}f_{z_1}\ldots f_{z_r}dz_1\ldots dz_r\newline &{}={}&\sum\limits_{r=1}^{\infty}\left(\dfrac{\lambda}{\lambda{}-{}t}\right)^r e^{t-\lambda}\left(1{}-{}e^{t-\lambda}\right)^{r-1}\newline &{}={}&e^{t-\lambda}\dfrac{\lambda}{\lambda{}-{}t}\sum\limits_{r=1}^{\infty}\left(\dfrac{\lambda}{\lambda{}-{}t}\right)^{r-1} \left(1{}-{}e^{t-\lambda}\right)^{r-1}\newline &{}={}&\left(\dfrac{e^{t-\lambda}\dfrac{\lambda}{\lambda{}-{}t}}{1{}-{}\left(\dfrac{\lambda}{\lambda{}-{}t}\right) \left(1{}-{}e^{t-\lambda}\right)}\right)\,\newline &{}={}&\dfrac{1}{1{}-{}e^{\lambda{}-{}t}t/\lambda}\,. \end{eqnarray*}$$

This looks like $\sum\limits_{n=1}^{\nu}\xi_{n}$ is trying to be exponentially distributed with "rate" $\lambda/e^{\lambda{}-{}t}$, but I do not know this functional form by heart. Any ideas?

Edit: Not knowing the explicit inverse of the final generating function form above, I thought of examining each term in the equivalent series representation: perhaps the individual terms have nicer inverses. If this is the case, then a series representation might be sufficient. If we perform the substitution $u{}={}t/\lambda$, so that $u<1$, note that the moment generating function may be re-written as

$$\sum\limits_{r=1}^{\infty}\left(\dfrac{1}{1-u}\right)^r\left(1-e^{\lambda\left(u-1\right)}\right)^{r-1}e^{\lambda\left(u-1\right)}{}={}\sum\limits_{r=1}^{\infty}\sum\limits_{k=0}^{r-1}{r-1\choose k}\left(\dfrac{1}{1-u}\right)^r(-)^ke^{(k+1)\lambda(u-1)}\,.$$

A series representation may be obtained, therefore, if we can invert the "atomic" moment generating functions

$$\left(\dfrac{1}{1-u}\right)^r e^{(k+1)\lambda(u-1)}\,.$$

Heuristically, we wish to solve an integral of the kind

$$\int\limits_{-\infty}^{1}\left(\dfrac{1}{1-u}\right)^r e^{(k+1)\lambda(u-1)-xu}\,\,\mbox{d}u\,.$$

For $x<\lambda(k+1)$, the integral's solution has the form $$(-\lambda)^{r}e^{-x}\left(\dfrac{x^{r-1}}{(r-1)!}\log(\lambda(k+1)-x){}+{}f_{r-1}(k)\right)$$

where $f_{r-1}(k)$ is a rational function involving terms of, at most, degree "$r-1$" in $k$.

(Note: the explicit solution can be obtained by integrating the expression $\dfrac{(-\lambda)^re^{-x}}{\lambda(k+1)-x}$, "$r$"-times, w.r.t "$k$". A justification of this follows by differentiating the integral we wish to solve. Note, also, that our "$u$" substitution above was merely to make this presentation look nicer and puts this solution "off" by a factor of $\lambda$: the actual solution follows analogous operations using the $t$ variable, instead).

• Even if your derived mgf/cf is does not have a tractable symbolic inversion to pdf (and even if it is not numerically invertible/stable) ... it still yields the very neat result that the expected sum ($E[Y]$, in my notation) is $e^\lambda/\lambda$ ... which is very elegant! Jan 10 '15 at 18:34

Let:

• $$X \sim \text{Exponential}(\lambda)$$ with pdf $$p(x) = \lambda e^{-\lambda x}$$, for $$x>0$$.

• $$(X_1, X_2, \dots)$$ denote successive random draws on $$X$$, where the sample $$(X_1, X_2, \dots, X_n)$$ is terminated as soon as $$X_n > 1$$ is attained.

• $$Y = \sum_{i=1}^n X_i$$

The stopped-sum constraint

The problem has two complications: The first is that the number of draws $$N=n$$ is itself a random variable. The second is that we are not just seeking the sum of Exponentials ... but rather that $$Y$$ is the sum of:

• $$(n-1)$$ 'truncated above' Exponentials (each conditional on $$X_i <1$$), PLUS:
• the final term $$X_n$$ which is a 'truncated below' Exponential (conditional on $$X_n >1$$).

This sum appears unlikely to have a tractable closed-form solution.

From theory to approximation

However, one can obtain very good approximations for small $$\lambda$$ (e.g. $$\lambda <\frac15$$) and large $$\lambda$$ (e.g. $$\lambda > 5$$). To see why, consider a plot of the Exponential pdf given different values of parameter $$\lambda$$:

Consider the two extremes:

# Case 1: $$\lambda$$ small

As $$\lambda \rightarrow 0$$ , $$P(X_i>1) \rightarrow 1$$, so in the extreme, the model should simplify to the pdf of $$(X \big| X>1)$$, which is a 'truncated below' Exponential with pdf:

$$\phi_{n=1}(y) = \frac{\lambda e^{-\lambda y}}{P(X>1)} = \frac{\lambda e^{-\lambda y}}{e^{- \lambda}} \quad \text{for } \quad y>1$$

How does the approximate solution do? Here is a comparison when $$\lambda= \frac15$$ of:

• true Monte Carlo empirical pdf of $$Y$$ (the squiggly BLUE curve), versus
• $$\phi_{n=1}(y)$$ - the truncated Exponential pdf (RED DASHED curve)

The fit appears excellent, even here where $$\lambda = \frac15$$ is not particularly small. The only substantive deviation is over $$y \in (1,2)$$, where the correct solution (Blue curve) appears to be approximately Uniform. (The roughly Uniform behaviour occurs for all values of $$\lambda$$, for $$1.)

Smaller values of $$\lambda$$ will yield an even better fit.

The solution thus far is:

• If $$n = 1 \text{:} \quad Y_1 = (X_1 \big| X_1>1)\quad \text{ with pdf } \phi_{n=1}(y)$$

For a better approximation, add more terms:

• If $$n = 2 \text{:} \quad Y_2 = (X_1 \big| X_1 <1) + (X_2 \big| X_2 >1) \quad \text{ which has pdf: }$$

Then, the order 2 approximation is:

$$P(N=1) \phi_{n=1}(y) + P(N>1) \phi_{n=2}(y)$$

where $$P(N=1) = P(X>1) = e^{-\lambda }$$. The small discrepancy over $$y \in (1,2)$$ is now resolved:

... as a zoomed in plot over that region of interest illustrates:

One can continue to add more terms as desired (I did $$n = 3$$ for fun) ... it might get a bit messy algebraically, but certainly possible with a computer algebra system. Fortunately, one does not need to add too many terms, because large $$n$$ is only needed for large $$\lambda$$, and for large $$\lambda$$, a simpler method exists ...

# Case 2: $$\lambda$$ large

If $$\lambda$$ is large, $$P(X_i< 1) \approx 1$$. In the extreme, the problem effectively reduces to finding the sum of $$n$$ iid Exponentials, which is well-known to be Gamma; in particular, that $$Y \sim \text{Gamma}(n, \frac{1}{\lambda})$$ (also known as the Erlang distribution when $$n$$ is an integer, as in our case), with pdf, say $$f(y)$$:

We then need the pmf of $$N=n$$:

\begin{align*}\displaystyle P(N=1) \quad &= \quad P(X_1>1) \\P(N=2)\quad &= \quad P(X_1 \leq 1) P(X_2>1) \\ P(N=3) \quad &= \quad P(X_1 \leq 1) P(X_2 \leq 1) P(X_3>1) \\ \dots \quad &= \quad \dots\\P(N=n) \quad &= \quad P(X \leq 1)^{n-1} P(X>1) = \left(1-e^{-\lambda }\right)^{n-1} e^{-\lambda } = \frac{\left(1-e^{-\lambda }\right)^n}{e^{\lambda }-1} \end{align*}

Let $$g(n)$$ denote the pmf of $$N$$:

Parameter-mix distribution: To solve, we require the expectation $$E_N[ Gamma(N,\frac{1}{\lambda})]$$ where $$N$$ has pmf $$g(n)$$. The pdf of the parameter-mix distribution is simply:

with domain of support on $$y>0$$. Finally, conditioning the latter on $$Y>1$$, we obtain the approximate pdf of $$Y$$, for large $$\lambda$$, as:

$$\phi(y) = \frac{\lambda}{\exp\left(-\lambda e^{-\lambda }\right)} \exp \left(-\lambda \left(1+ e^{-\lambda } y\right)\right) \quad \quad \text{ for} \quad y>1$$

How does the approximate solution do? Here is a comparison when $$\lambda=5$$ of:

• true Monte Carlo empirical pdf of $$Y$$ (the squiggly BLUE curve), versus
• $$\phi_(y)$$ - the Gamma parameter-mix pdf (RED DASHED curve)

... and a 'zoomed-in' version of the same:

Again, the fit appears excellent, even here where $$\lambda = 5$$ is not particularly large. Larger values of $$\lambda$$ will yield an even better fit.

Summary

If $$\lambda$$ is small (e.g. $$\frac15$$ or smaller) or large (e.g. 5 or bigger), then approximate limit solutions that have fairly simple forms appear to yield good solutions. For mid-range values of $$\lambda$$ in-between, a suitable weighted average of the two limit solutions may yield reasonable results, or more terms $$\phi_{n=3}(y)$$ and $$\phi_{n=4}(y)$$ etc can be added.

Notes

• Monte Carlo simulation of the stopped-sum process:

One standard approach would be to write a recursive function such as:

 Func := (rr = MrRandom;  AppendTo[xvals, rr];  If[  rr > 1, xvals, Func])


where

MrRandom := RandomReal[ExponentialDistribution[2]]


and then call 6 runs with, say:

Table[xvals = {}; Func, {6}]


However, for large values of $$\lambda$$, the number of pseudorandom drawings required to attain an $$x_i > 1$$ can be very large indeed, which may cause problems with iteration limits and recursive limits etc. It is also a slow way to proceed. A much nicer way, here using Mathematica, is to generate the samples $$(x_1, \dots, x_n)$$ as a one-liner:

Split[ RandomReal[ExponentialDistribution[.1], 10^7], # < 1 &]


... which generates 10 million Exponential pseudorandom drawings (in one go), and then splits them into separate samples whenever a value greater than 1 is attained.

• The Expect function used above is from the mathStatica package for Mathematica. As disclosure, I should add that I am one of the authors.
• Good!But I don't understand why you felt the mathematical expectation $E_N[ Gamma(N,\frac{1}{\lambda})]$? Can you explain why? Jan 7 '15 at 13:53
• See the section on parameter-mix distributions in Chapter 3 (section 3.4B) of our book mathstatica.com/book/bookcontents.html (free download) Jan 7 '15 at 15:39
• This would be valid for the random sum $\sum\limits_{i=1}^NX_i$ with $N$ independent of $(X_i)$, which is not the case.
– Did
Jan 7 '15 at 19:10
• Thanks - well-spotted: will give it some more thought. Jan 8 '15 at 4:04
• @Did Modified to obtain 2 approximate limit solutions ... Jan 9 '15 at 18:27

There is a subtle dependence betwen $\nu$ and the sums of random variables $\xi_k$, in particular $\nu$ is not independent of $(\xi_k)$. However, note that:

• $g\gt1$ almost surely
• if $\xi_1\gt1$ then $g=\xi_1$
• if $\xi_1\lt1$ then $g=\xi_1+g'$ where $g'$ is independent of $\xi_1$ and distributed like $g$

Thus, for every $x\gt1$, $$P(g\gt x)=P(\xi\gt x)+\int_0^1\mathrm dP_{\xi}(t)P(g\gt x-t),$$ that is, $$P(g\gt x)=\mathrm e^{-\lambda x}+\int_0^1\lambda\mathrm e^{-\lambda t}P(g\gt x-t)\mathrm dt.$$ From this point, one can show that there exists some sequence $(p_n)$ of polynomials such that, for every nonnegative integer $n$ and every $u$ in $[0,1]$, $$P(g\gt n+u)=p_n(u).$$ The sequence $(p_n)$ is uniquely determined by the recursion $$p'_{n+1}(u)=-\lambda\mathrm e^{-\lambda}p_n(u),\quad p_{n+1}(0)=p_n(1),$$ for every $n\geqslant0$, with the initial condition $$p_0(u)=1.$$ Thus, every $p_n(u)$ depends on the parameter $$a=\lambda\mathrm e^{-\lambda},$$ and the generating function of $(p_n)$ is $$\sum_{n\geqslant0}p_n(u)x^n=\frac{\mathrm e^{-axu}}{1-x\mathrm e^{-ax}}.$$ Explicit general formulas for $p_n$ are not so easy to deduce from this but one sees that each polynomial $p_n$ has degree $n$, and even that $$p_n(u)=1-nu+[\text{some monomials from u^2 to u^{n-1}}]+(-1)^n\frac{u^n}{n!}.$$ Nevertheless, note as an example of application that, for every $\lambda\leqslant1$, the smallest positive root of $x=\mathrm e^{ax}$ is $x=\mathrm e^\lambda$ hence $$p_n(u)\approx\mathrm e^{-n\lambda}.$$