What is the probability that $\sum_{k=1}^\infty u_k>1$, if $u_1=$ random real number in $(0,1)$, and $u_k=$ random real number in $(0,u_{k-1})$?

Question

Let $u_1=$ i.i.d.$\text{Uniform}(0,1)$-variable.

Let $u_k=$ i.i.d.$\text{Uniform}(0,u_{k-1})$-variable for $k>1$.

What is the probability that $\sum\limits_{k=1}^\infty u_k>1$ ?

The expectation of $\sum\limits_{k=1}^\infty u_k$ is $1$. This can be shown as follows. For $a_i=$ i.i.d.$\text{Uniform}(0,1)$-variable, we have:

$$E\left(\sum\limits_{k=1}^\infty u_k\right)=E\left(\sum\limits_{k=1}^\infty \left(\prod\limits_{i=1}^k a_i\right)\right)=\sum\limits_{k=1}^\infty \left(\prod\limits_{i=1}^k E(a_i)\right)=\sum\limits_{k=1}^\infty \left(\prod\limits_{i=1}^k \left(\frac12\right)\right)=1$$

An Excel simulation suggests that $P\left(\sum\limits_{k=1}^{50} u_k>1\right)$ and $P\left(\sum\limits_{k=1}^{100} u_k>1\right)$ are both roughly $0.44$. Other than that, I do not know how to approach this problem.

Context: This question was inspired by another question involving the sum of random real numbers.

Answer 1

We claim that the desired probability is

$$ \boxed{ 1 - e^{-\gamma} \approx 0.438541 }, $$

where $\gamma$ is the Euler–Mascheroni constant.

---

Step 1. Let $Z_1, Z_2, \ldots$ be i.i.d. $\text{Uniform}(0, 1)$-variables. Then the sum in question can be realized as

$$ S = \sum_{k=1}^{\infty} Z_1 \cdots Z_k. \tag{*} $$

From this, it is easy to check that

$$ S \stackrel{d}= U(1 + S), \tag{1} $$

where $U$ is a $\text{Uniform}(0, 1)$-variable independent of $S$ and $\stackrel{d}=$ stands for equality in distribution. So, if $F(\cdot)$ denotes the CDF of $S$, then for $x > 0$,

\begin{align*} F(x) = \mathbf{P}(S \leq x) = \int_{0}^{1} \mathbf{P}(u(1+S) \leq x) \, \mathrm{d}u \stackrel{(u=x/t)}= x \int_{x}^{\infty} t^{-2} F(t - 1) \, \mathrm{d}t. \tag{2} \end{align*}

Step 2. The equation $\text{(2)}$ shows that $F(x)$ is absolutely continuous. Moreover, for $x \in (0, 1]$,

$$ F(x) = x \int_{1}^{\infty} t^{-2} F(t - 1) \, \mathrm{d}t = F(1)x $$

and hence the density $f(x) = F'(x)$ satisfies $f(x) = F(1)$ for $x$ between $0$ and $1$.

Step 3. Let $\mu_n = \mathbf{E}[S^n]$ be the $n$th moment of $S$, and let $M(\xi) = \mathbf{E}[e^{\xi S}] = \sum_{n=0}^{\infty} \frac{\mu_n}{n!} \xi^n$ denote the moment generating function of $S$. Then

\begin{align*} \mu_n = \mathbf{E}[S^n] = \mathbf{E}[U^n(1 + S)^n] = \frac{1}{n+1} \sum_{k=0}^{n} \binom{n}{k} \mu_k. \end{align*}

From this, we find that $M(\cdot)$ solves the differential equation

$$ (\xi M(\xi))' = e^{\xi} M(\xi). $$

Together with the initial condition $M(0) = 1$, this yields

$$ M(\xi) = \exp \left( \int_{0}^{1} \frac{e^{\xi t} - 1}{t} \, \mathrm{d}t \right). \tag{3} $$

Step 4. Finally, the main claim will follow once we show that $\mathbf{P}(S \leq 1) = e^{-\gamma}$. However,

\begin{align*} \mathbf{P}(S \leq 1) = F(1) &= f(0^+) \\ &= \lim_{r \to \infty} \int_{0}^{\infty} r e^{-rt} f(t) \, \mathrm{d}t \\ &= \lim_{r \to \infty} r \mathbf{E}[e^{-rS}] \\ &= \lim_{r \to \infty} r \exp \biggl( \int_{0}^{1} \frac{e^{-r t} - 1}{t} \, \mathrm{d}t \biggr) \tag{by (3)} \\ &= \lim_{r \to \infty} \exp \biggl( \int_{0}^{r} \frac{e^{-s} - \mathbf{1}_{[0,1]}(s)}{s} \, \mathrm{d}s \biggr) \tag{$s=rt$} \\ &= \exp \biggl( \int_{0}^{\infty} \frac{e^{-s} - \mathbf{1}_{[0,1]}(s)}{s} \, \mathrm{d}s \biggr). \end{align*}

It is not hard to show that the integral in the last line is equal to $-\gamma$, and therefore the desired conclusion follows. $\square$

---

Extra Observations. Just for fun, here is an alternative description of the distribution of $S$. Let $\nu$ be a Poisson point process of unit intensity on $\mathbb{R}^2$, and let $N(A) = \nu(A\times[0,1])$ so that $N$ is a Poisson point process of unit intensity on $\mathbb{R}$. Then the sum $S$ defined in $\text{(*)}$ satisfies

$$ S \stackrel{d}= \int_{(-\infty, 0]} e^{x} \, N(\mathrm{d}x) = \int_{(-\infty, 0]\times[0, 1]} e^{x} \, \nu(\mathrm{d}x\mathrm{d}y). $$

Now there are two ways to extend the distribution of $S$:

1. We may let the vertical direction of the domain of integration vary. So, define $$ L(t) = \int_{(-\infty, 0]\times[0, t]} e^{x} \, \nu(\mathrm{d}x\mathrm{d}y). $$ Then $L$ is a subordinator with drift zero and Lévy measure $x^{-1}\mathbf{1}_{(0,1)}(x) \, \mathrm{d}x$ such that $L(1)$ has the same distribution as $S$. In particular, the distribution of $S$ is infinitely divisible. However, I am not sure if this will help understand the original problem.

2. Instead, we may let the horizontal direction of the domain to vary. So, let $$ X(t) = \int_{(-\infty, t]} e^{x-t} \, N(\mathrm{d}x). $$ Then $X$ is a stationary process (since $X(t) \sim S$ for any $t$) that solves the SDE of the form $$ \mathrm{d}X(t) = -X(t) \, \mathrm{d}t + \mathrm{d}N(t). $$ Slightly abusing notation, this tells that MGF $M(\cdot)$ of $S$ satisfies \begin{align*} M(\xi) &= \mathbf{E}[\exp(\xi X(t+\mathrm{d}t))] \\ &= \mathbf{E}[\exp(\xi (1-\mathrm{d}t)X(t) + \xi \, \mathrm{d}N(t))] \\ &= M(\xi - \xi \, \mathrm{d}t) (1 + (e^{\xi} - 1) \, \mathrm{d}t) \\ &= M(\xi) + [ (e^{\xi} - 1) M(\xi) - \xi M'(\xi) ] \, \mathrm{d}t \end{align*} and hence reproduces the differential equation for $M(\cdot)$, $$ (e^{\xi} - 1) M(\xi) - \xi M'(\xi) = 0. $$ This approach also gives a more transparent view of the phenomenon in which the CDF $F(\cdot)$ of $S$ solves a delay integral equation $\text{(2)}$ and hence is given by different formulas on different intervals of the form $[k-1, k]$, $k = 1, 2, 3, \ldots$. Indeed, this seems related to the occasional jump $\mathrm{d}N(t)$ of unit size.

   For me, this seems hinting that the process $X$ can be used to study the distribution of $S$. But again, I don't have any good idea in this direction now.