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)$?
 A: 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).   
A: 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<y<\approx 2$.)
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.

A: 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}.$$
