4
$\begingroup$

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)$?

$\endgroup$
4
  • $\begingroup$ Is seems like you are implicitly assuming that the $\xi_n$ are ordered $\endgroup$
    – user76844
    Jan 7 '15 at 1:28
  • $\begingroup$ 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/… $\endgroup$
    – Math1000
    Jan 8 '15 at 22:25
  • $\begingroup$ @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). $\endgroup$
    – wolfies
    Jan 14 '15 at 15:34
  • $\begingroup$ 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 :( $\endgroup$
    – Math1000
    Jan 14 '15 at 15:55
1
$\begingroup$

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

$\endgroup$
1
  • $\begingroup$ 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! $\endgroup$
    – wolfies
    Jan 10 '15 at 18:34
1
$\begingroup$

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.
$\endgroup$
5
  • $\begingroup$ Good!But I don't understand why you felt the mathematical expectation $E_N[ Gamma(N,\frac{1}{\lambda})]$? Can you explain why? $\endgroup$
    – user202790
    Jan 7 '15 at 13:53
  • $\begingroup$ See the section on parameter-mix distributions in Chapter 3 (section 3.4B) of our book mathstatica.com/book/bookcontents.html (free download) $\endgroup$
    – wolfies
    Jan 7 '15 at 15:39
  • $\begingroup$ 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. $\endgroup$
    – Did
    Jan 7 '15 at 19:10
  • $\begingroup$ Thanks - well-spotted: will give it some more thought. $\endgroup$
    – wolfies
    Jan 8 '15 at 4:04
  • $\begingroup$ @Did Modified to obtain 2 approximate limit solutions ... $\endgroup$
    – wolfies
    Jan 9 '15 at 18:27
0
$\begingroup$

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

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.