Expected number of iterations till a walk on unit line exceeds a distance

Question

Problem

Each timestep $i$, a uniformly distributed next target $X_i$ is sampled in the range $0$ to $1$. We start out at point $X_0$ on the number line. Each timestep, we go from our current point to the next sampled point. What is the expected value of $N$, where $N$ is the first timestep when the total distance covered exceeds $1$ unit? In otherwords, I am looking for $E[N]$ where $$ N = \min\left\{n\in\mathbb{N}: \sum\limits_{i=1}^{n}{|X_{i}-X_{i-1}|}\geq 1\right\} \text{.} $$

I have tried several different approaches already (listed below), to no avail.

Tried Approaches

The obvious approach is observing that the expected absolute difference between $2$ uniform random variables is $\frac{1}{3}$. The reciprocal would give an answer of $3$ iterations. But this doesn't match the simulation results, I guess because it ignores the fact that each consecutive absolute difference is correlated. For example, if the previous absolute difference was high, the next absolute difference is more likely to be higher as well.

For example, we could define $T(x)$ to be the expected number of steps needed to cover $x$ more distance (looking for $T(1)$). But this recurrence relation

$$T(x)=\begin{cases}1+\int_{0}^{1}{T(x-|y-x|)dy} & \text{, }x\gt 0 \\ 0 & \text{, otherwise}\end{cases}$$

doesn't seem solvable. The generating function associated with $T(x)$ doesn't seem to be able to be simplified. Moreover, I don't know how to solve its corresponding differential equation $T'(x)=2(T(x)-T(2x-1))$.

Another approach was to define $Y_i=|X_{i}-X_{i-1}|$, then try to find this conditional CDF $F(y|x)$, differentiate to get the conditional PDF $f(y|x)$, then evaluate $\int_0^1{yf(y|x)dy}$ (representing $E[Y_i|X_i=x]$). But its reciprocal is also definitely not the correct answer.

Simulation

Via simulation, I know the approximate answer is 3.82, but I want a closed form or exact expression.

#!/usr/bin/env python3
import math
import random
NUM_EXPERIMENTS = 100_000
def main():
    sim_answer = simulate_answer_1()
    print('sim: ' + str(sim_answer))
def simulate_answer_1() -> float:
    avg_iter_before_one = 0
    for experiment in range(NUM_EXPERIMENTS):
        dist_covered = 0
        curr = random.random()
        this_n = 0
        while dist_covered < 1:
            next = random.random()
            dist_covered += abs(next - curr)
            this_n += 1
            curr = next
        avg_iter_before_one += this_n
    return avg_iter_before_one / NUM_EXPERIMENTS
if name == 'main':
    main()

Answer 1

Closed form solution

Let $D_n = \sum_{i=1}^n |X_i-X_{i-1}|$. Then we obtain $$P(D_n < 1) = \frac{2^n}{n!} + \frac{2^{n+1}}{(n+1)!} - \frac{2(2n+1)!!}{(n+1)!n!}$$

Since a generating function for this formula is $2e^{2x} -2e^x(I_0(x) + I_1(x)) - 1$ we have the following expression for the expected number of steps: $$\mathbb{E}[N] = 2e^2 - 2e(I_0(1)+I_1(1))-1 = 3.822541014468694\!\ldots$$ where $I_k$ is a Modified Bessel function of the first kind.

Table for the first values of the probability \begin{array} {|c|c|c|c|c|c|c|} \hline n & 1 & 2 & 3 & 4 & 5 & 6 & 7 \\ \hline P & 1 & 5/6 & 13/24 & 133/480 & 331/2880 & 107/2688 & 761/64512\\ \hline \end{array} A proof I have for this is not elegant.

---

1. We start with the idea of Zoe Allen to consider a probability density function $g_n(x, y)$ after $n$ steps, where $x$ is a current position and $y$ is a total distance traveled. Then $$P(D_n<1) = \int_0^1 \int_0^1 g_n(x, y) dxdy$$ We also have the following recurrence: $$g_{n+1}(x, y) = \int_0^1 g_n(t, y - |x-t|)dt$$ We divide the region $[0, 1] \times [0, 1]$ into four quadrants: \begin{align}\color{blue}{\text{I}}&: {y \geq x, y \geq 1-x} &\quad \color{blue}{\text{II}}&: {y \geq x, y \leq 1-x}\\ \color{blue}{\text{III}}&: {y \leq x, y \geq 1-x} &\quad \color{blue}{\text{IV}}&: {y \leq x, y \leq 1-x}\end{align} We are going to use the notation $g^\color{blue}{\text{K}}(x, y)$ which means that $x$ and $y$ are allowed only from the quadrant $\color{blue}{\text{K}}$ (so that we can have different functions for each quadrant). We can then rewrite the recurrence as follows: $$g^\color{blue}{\text{I}}_{n+1}(x, y) = \int_0^{(1-y+x)/2} g^\color{blue}{\text{II}}_n(t, y-x+t) dt + \int_{(1-y+x)/2}^x g^\color{blue}{\text{I}}_n(t, y-x+t) dt + \int_x^{(x+y)/2} g^\color{blue}{\text{I}}_n(t, y+x-t) dt + \int_{(x+y)/2}^1 g^\color{blue}{\text{III}}_n(t, y+x-t) dt$$ $$g^\color{blue}{\text{II}}_{n+1}(x,y) = \int_0^x g^\color{blue}{\text{II}}_n(t, y-x+t) dt + \int_x^{(x+y)/2} g^\color{blue}{\text{II}}_n(t, y+x-t) + \int_{ (x+y)/2}^{x+y} g^\color{blue}{\text{IV}}_n(t, y+x-t)$$ $$g^\color{blue}{\text{III}}_{n+1}(x,y) = \int_{x-y}^{(1-y+x)/2} g^\color{blue}{\text{IV}}_n(t, y-x+t) dt + \int_{(1-y+x)/2}^{x} g^\color{blue}{\text{III}}_n(t, y-x+t) dt + \int_x^1 g^\color{blue}{\text{III}}_n(t, y+x-t)$$ $$g^\color{blue}{\text{IV}}_{n+1}(x,y) = \int_{x-y}^x g^\color{blue}{\text{IV}}_n(t, y-x+t) dt + \int_x^{x+y} g^\color{blue}{\text{IV}}_n(t, y+x-t)$$ Looks dreadful but it is actually much better now. The limits in the integrals are constructed in such a way, that each function remains in its own region during the integration. We can also clearly see that if each of the $g^\color{blue}{\text{K}}_n$ is polynomial for some $n$, then they would remain polynomial after the integration. Now, after the first step we can obtain $g^\color{blue}{\text{I}}_1(x,y) = 0$, $g^\color{blue}{\text{II}}_1(x,y) = 1$, $g^\color{blue}{\text{III}}_1(x,y) = 1$, $g^\color{blue}{\text{IV}}_1(x,y) = 2$. Which means that all $g^\color{blue}{\text{K}}$ in the next steps will be polynomials (which is not true for the whole region). This provides a straightforward tool for computing the next density function by simply integrating polynomials from the previous step, which is already enough to get a non-simulated answer for relatively large $n$.

2. Since the rules of the transformation are relatively simple, we can hope to find the general formula for the polynomials. It is also nice that, due to symmetry of the problem, we have $g^\color{blue}{\text{II}}(x,y) = g^\color{blue}{\text{III}}(1-x,y)$, which means that it is enough to find only one of these polynomials. Another observation is that since $g^\color{blue}{\text{IV}}$ has a recurrence relation only with itself, we could expect a simpler solution for this region. To formulate the problem in terms of the recurrence relationships between the coefficients of polynomials, it is enough to derive the expressions for the transformation of each monomial $x^i y^j$ under the action of these integrals (and the polynomials will follow by linearity). Note, that monomials do not necessary transform to monomials, but can also transform to polynomials, so this is a bit tricky.

Example. Suppose we want to derive how a monomial $xy^2$ from a quadrant $\color{blue}{\text{II}}$ transforms under the action of the integral for the quadrant $\color{blue}{\text{I}}$. Then we can write $$\int_0^{(1-y+x)/2} t(y-x+t)^2 dt = \int_0^{(1-y+x)/2} (t^3 - 2t^2x + 2t^2y + tx^2 -2xyt + ty^2)dt$$ $$=\frac14 t^4 - \frac23 t^3x + \frac23 t^3y + \frac12 t^2 x^2 - xyt^2 + \frac12 t^2 y^2 \biggl|_0^{(1-y+x)/2}$$ $$=\frac{1}{64} (1-y+x)^4 - \frac{1}{12} (1-y+x)^3x + \frac{1}{12} (1-y+x)^3y + \frac18 (1-y+x)^2 x^2 - \frac14 xy(1-y+x)^2 + \frac18 (1-y+x)^2 y^2$$ After expanding the remaining powers of $(1-y+x)$ we obtain a new polynomial in $x, y$.

The key idea which allows to simplify this is to consider polynomials in new variables $f_n(v, w) = g_n(x, y)$ where: $$v = \frac12 (x-y+1),\quad w = \frac12 (x+y-1)$$ This makes the transformation rules much simpler.

Example. Suppose we want to derive how a monomial $v^iw^j$ from a quadrant $\color{blue}{\text{II}}$ transforms under the action of the integral for the quadrant $\color{blue}{\text{I}}$. Then $$\int_0^v \frac{(t-(y-x+t)+1)^i(t+y-x+t-1)^j}{2^{i+j}}=$$ $$=\int_0^v v^i (t-v)^j = v_i \frac{(t-v)^{j+1}}{j+1}\biggl|_0^v = \frac{v^{i+j+1}(-1)^j}{j+1}$$

Similarly, we can find such rules for other integrals: \begin{array} {|c|c|c|c|c|} \hline & \color{blue}{\text{I}} & \color{blue}{\text{II}} & \color{blue}{\text{III}} & \color{blue}{\text{IV}} \\ \hline \color{blue}{\text{I}} & \frac{v^{i}w^{j+1}}{j+1} + \frac{w^j}{2^{i+1}(i+1)} - \frac{v^{i+1}w^j}{i+1} & \frac{v^{i+j+1}(-1)^j}{j+1} & \frac{w^j(1-w)^{i+1}}{i+1} - \frac{w^j}{2^{i+1}(i+1)} & -\\ \hline \color{blue}{\text{II}} &-&\frac{v^iw^{j+1}}{j+1} + \frac{v^{i+j+1} (-1)^j}{j+1} + \frac{w^j}{2^{i+1}(i+1)} - \frac{v^{i+1}w^j}{i+1}&-&\frac{w^j(1+w)^{i+1}}{i+1} - \frac{w^j}{2^{i+1}(i+1)}\\ \hline \color{blue}{\text{III}} & -&- & \frac{v^{i}w^{j+1}}{j+1} + \frac{w^j(1-w)^{i+1}}{i+1} - \frac{v^{i+1}w^j}{i+1}& - \frac{v^i(v-1)^{j+1}}{j+1}\\ \hline \color{blue}{\text{IV}} & -& -&- & \frac{v^iw^{j+1}}{j+1} - \frac{v^i(v-1)^{j+1}}{j+1}+ \frac{w^j(1+w)^{i+1}}{i+1} - \frac{w^j v^{i+1}}{i+1}\\ \hline \end{array} where the rows correspond to the quadrant the monomial comes from and the columns are for which quadrant we are integrating. The terms with $1+w$, $1-w$ and $1-v$ should be additionally expanded using Binomial theorem, but I will not write it here for the sake of brevity.

This is already a substantial improvement over the previous step since now we do not need to integrate anything to find the polynomials (only at the final step to obtain the probability). Fortunately, we also don't need to convert the polynomials back to the original coordinates because we can just write the area integrals for the probability in new coordinates: $$p(D_n \leq 1) = 2 \left( \int_0^{1/2} \int_0^v f_n^\color{blue}{\text{I}} (v, w) + 2\int_0^{1/2} \int_{-v}^0 f_n^\color{blue}{\text{II}} (v, w) + \int_{1/2}^1 \int_{v-1}^0 f_n^\color{blue}{\text{IV}} (v, w)\right)$$ where the factor of $2$ is due to Jacobian determinant and we do not have a quadrant $\color{blue}{\text{III}}$ integral due to symmetry.

3. The table above is essentially recurrence formulas for the coefficients of polynomials. I was able to guess the general solution by examining the dependence on $n$ for monomials with small powers of $i$ and $j$: $$[v^i w^j]f^\color{blue}{\text{I}}_n(v, w) = \frac{(-1)^i}{i! j! (n-1-i-j)!} \left[R_n(i, j, 1) -Q_n(j)\right]$$ $$[v^i w^j]f^\color{blue}{\text{II}}_n(v, w) = \frac{(-1)^i}{i! j! (n-1-i-j)!} R_{n+1}(i, j, 0)$$ $$[v^i w^j]f^\color{blue}{\text{IV}}_n(v, w) = \frac{(-1)^i}{i! j!(n-1-i-j)!} 2^n$$ where $$R_n(i, j, q) = \frac{(-1)^{\min(i,j)}}{2^{n-1-i-j}}\sum_k \begin{pmatrix} \max(i, j+1)-k-1\\k-\min(i, j+1)\end{pmatrix} \begin{pmatrix}2(n-1+k)+q\\n-1+k\end{pmatrix} (-1)^k$$ $$Q_n(j) = \sum_{k=0}^j \begin{pmatrix}n+1\\k\end{pmatrix}$$ From now the problem reduces to calculating grid sums with respect to $(i, j)$ and integrating to find the probability. We can do this explicitly. Here is an example for the quadrant $\color{blue}{\text{IV}}$, which is the easiest: $$p_n^\color{blue}{\text{IV}} = 2^{n+1} \int_{1/2}^1 \int_{v-1}^0 \sum_{0\leq i+j\leq {n-1}} \frac{(-1)^i v^i w^j}{i! j! (n-1-i-j)!} = 2^{n+1} \int_{1/2}^1 \sum_{0\leq i+j\leq {n-1}} \frac{(-1)^{i+1} v^i (v-1)^{j+1}}{i! (j+1)! (n-1-i-j)!} = $$ $$=2^{n+1} \int_{1/2}^1 \sum_{0\leq i+j\leq {n-1}} \sum_{k=0}^{j+1} \frac{(-1)^{i+j-k} v^{i+k}}{i! k! (j+1-k)!(n-1-i-j)!} = 2^{n+1} \sum_{0\leq i+j\leq {n-1}} \sum_{k=0}^{j+1} \frac{(-1)^{i+j-k} v^{i+k+1}}{i! k! (i+k+1)(j+1-k)!(n-1-i-j)!} \biggl|_{1/2}^1 = $$ $$=2^{n+1} \sum_{i=0}^{n-1} \sum_{k=0}^{n-1-i+1} \frac{(-1)^{i-k} v^{i+k+1}}{(i+k+1)i! k!} \sum_{\substack{j=k-1\\j\geq 0}}^{n-1-i} \frac{(-1)^j}{(j+1-k)!(n-1-i-j)!} \biggl|_{1/2}^1 = $$ $$=2^{n+1} \left[\sum_{i=0}^{n-1} \frac{(-1)^{i-1} v^{n+1}}{(n+1)i!(n-1-i+1)!} + \sum_{i=0}^{n-1} \frac{(-1)^i v^{i+1}}{(i+1)!} \frac{1}{(n-1-i+1)(n-1-i)!}\right] \biggl|_{1/2}^1 = $$ $$=\frac{2^{n+1}}{(n+1)!} \left[ 1-(1-v)^{n+1} \right] \biggl|_{1/2}^1 = \frac{1}{(n+1)!}$$ The complete solution is $$p(D_n<1) = \frac{1}{(n+1)!}\left\{1 + \frac{1}{2^{n+1}} \left[\sum_{k=0}^{\lfloor n/2 \rfloor} \sum_{t=0}^{n-1-2k} \begin{pmatrix}2(n-k)\\n-k\end{pmatrix} \begin{pmatrix}t+k\\k\end{pmatrix} (-1)^k 2^{2k+t+1} + \sum_{k=0}^{\lfloor n/2 \rfloor} \begin{pmatrix}2(n-k)\\n-k\end{pmatrix} \begin{pmatrix}n-k\\k+1\end{pmatrix}(-1)^k \right.\right.$$ $$\left.\left.-2\sum_{k=0}^{n-1} \begin{pmatrix}n+1\\k\end{pmatrix} \begin{pmatrix}n+1\\k+2\end{pmatrix} - \sum_{k=0}^{n-1} \sum_{t=0}^{n-1} \begin{pmatrix}n+1\\k\end{pmatrix} (t+4)2^{t} \begin{pmatrix}n-t\\n-k-2-t\end{pmatrix} \right]\right\}$$ which gives an explicit expression for the probabily.
4. We can simplify these sums further with a help of Mathematica and Egorychev method where it fails. I will just state the results: $$\sum_{k=0}^{\lfloor n/2 \rfloor} \sum_{t=0}^{n-2k-1} \begin{pmatrix}2(n-k)\\n-k\end{pmatrix} \begin{pmatrix}t+k\\k\end{pmatrix} (-1)^k 2^{2k+t+1} = 2(n+1)\left(4^n-\begin{pmatrix}2n+1\\ n\end{pmatrix}\right)$$ $$\sum_{k=0}^{\lfloor n/2 \rfloor} \begin{pmatrix}2(n-k)\\n-k\end{pmatrix} \begin{pmatrix}n-k\\k+1\end{pmatrix}(-1)^k = 2 \left( \begin{pmatrix}2n+1\\n\end{pmatrix}-2^n\right)$$ $$\sum_{k=0}^{n-1} \begin{pmatrix}n+1\\k\end{pmatrix} \begin{pmatrix}n+1\\k+2\end{pmatrix}=\frac{2^{n+1}n(n+1)(2n+1)!!}{(n+3)!}$$ $$\sum_{k=0}^{n-1} \sum_{t=0}^{n-1} \begin{pmatrix}n+1\\k\end{pmatrix} \begin{pmatrix}n-t\\n-k-t-2\end{pmatrix} (t+4)2^{t} = \frac{3n+4}{n+2}\begin{pmatrix}2n+2\\n+1\end{pmatrix} - 2\begin{pmatrix}2n+2\\n-1\end{pmatrix} + n\begin{pmatrix}2n+2\\n\end{pmatrix} - 2^{2n+2}$$ After some algebraic manipulations the result simplifies to: $$P(D_n < 1) = \frac{2^n}{n!} + \frac{2^{n+1}}{(n+1)!} - \frac{2(2n+1)!!}{(n+1)!n!}$$

Answer 2

Update:

I am not able to come up with an expression for the desired expected value. One approach that looks useful is splitting the region $D_n<1$ by permutations of $X_i$, and integrating over an $n$-D basis including $\begin{bmatrix}1&1&\cdots&1\end{bmatrix}^T$ and the coefficients for the expression for $D_n$ in that permutation. This isn't too hard to do for a given permutation, but it turns out that permutations with the same number of terms in $D_n$ can have different hypervolumes. It might be possible to find a stricter notion of equivalence for the volumes of permutations, but it would still need to be the case that a pattern emerges for the values of the integrals.

Using simulations, I was able to find that the first few terms of the summation are

$$P[D_n<1]=\frac{a_n}{n!2^{n-3}}$$ $$a_{n,n\geq2}=[1,5,26,133,662,3210,15220,274000\pm8,\dots]$$

Evaluating the sum as described below with these $P[D_n<1]$ up to $n=9$, we get $3.8186$, which is consistent with your simulations.

---

Answer to related problem, assumes $X_0=0$:

I have come up with an analytic expression for the expected value. It might be possible that my expression can be simplified with a greater understanding of group theory, but I highly doubt that the infinite sum can be resolved into a closed-form expression.

First of all, let $D_n=\sum_{i=1}^n\lvert X_i-X_{i-1}\rvert$ be the random variable which is the total distance traveled after $n$ steps. As a step toward finding $E[N]$, we want to find $P[D_n<1]$.

Now, $D_n$ has all those nasty absolute values, but we can get rid of them by splitting $P[D_n<1]$ into every possible ordering-by-magnitude of $X_i$:

$$P[D_n<1] = P[(X_1>X_2>\cdots>X_n) \wedge (D_n<1)] + \cdots +P[(X_{i_1}>X_{i_2}>\cdots>X_{i_n}) \wedge (D_n<1)]+\cdots$$

Let $S_n$ be the set of permutations of $n$ elements, and let $\sigma=[i_1,i_2,\dots,i_n]\in S_n$. Then we can let $I(\sigma)$ be the event that $X_{i_1}>X_{i_2}>\cdots>X_{i_n}$. Also let $\sigma[j]=i_j$.

Let us now define $\epsilon_i(\sigma)$ to have unity magnitude, but be positive whenever $X_i>X_{i-1}$, and negative when $X_{i-1}>X_i$. This can be done by defining $\sigma^{-1}$ to be the inverse permutation of $\sigma$, and defining

$$\epsilon_1(\sigma)=1;\qquad\epsilon_i(\sigma)=\left\{\begin{matrix}1&{\rm when\ }\sigma^{-1}[i]<\sigma^{-1}[i-1]\\-1&{\rm when\ }\sigma^{-1}[i]>\sigma^{-1}[i-1]\end{matrix}\right.,\ 2\leq i\leq n$$

This gives us that for a given $\sigma$, $D_n=\sum_{i=1}^n\epsilon_i(\sigma)(X_i-X_{i-1})$. We can simplify this a little more by defining

$$c_i(\sigma)=\epsilon_i(\sigma)-\epsilon_{i+1}(\sigma),\ 1\leq i<n;\qquad c_n(\sigma)=\epsilon_i(\sigma)$$

Now for a given $\sigma$ we have

$$D_n=\sum_{i=1}^nc_i(\sigma)X_i$$

We also have that

$$P[D_n<1]=\sum_{\sigma\in S_n}P\left[I(\sigma)\wedge\sum_{i=1}^nc_i(\sigma)X_i<1\right]$$

Now to find $P\left[I(\sigma)\wedge\sum_{i=1}^nc_i(\sigma)X_i<1\right]$, we observe that the event is a subset of the $n$-D unit hypercube. It turns out that this region is an $n$-simplex (eg. triangle, tetrahedron, pentachoron) with vertices at $x_0,x_1,\dots,x_n$, where $x_k$ has the coordinates in dimensions $\sigma[j],\ j\leq k$ equal to the same constant, and the rest zero. Additionally, the vertices $x_1,\dots,x_n$ satisfy the equation

$$\sum_{i=1}^nc_i(\sigma)x_{ki}=1,\qquad x_k=\begin{bmatrix}x_{k1}\\x_{k2}\\\vdots\\x_{kn}\end{bmatrix}$$

Thus we have that

$$x_{kj}=\frac{1}{\sum_{i=1}^kc_{\sigma[i]}(\sigma)}\left\{\begin{matrix}1&{\rm when\ }\sigma^{-1}[j]\leq k\\0&{\rm when\ }\sigma^{-1}[j]> k\end{matrix}\right.$$

The hypervolume of the $n$-simplex is known to be $\frac{1}{n!}$ the volume of the parallelotope (eg. parallelogram, parallelepiped) with the same edges $x_k-x_0$, and can thus be calculated as

$$V_\sigma=\frac{1}{n!}\left\lvert\det\begin{bmatrix}x_1&x_2&\cdots&x_n\end{bmatrix}\right\rvert$$

Now this matrix can be transformed with row rearranging into a triangular matrix, with one non-zero element from each $x_k$ on the diagonal. Then the determinant is the product of the elements on the diagonal, so

$$V_\sigma=\frac{1}{n!}\prod_{k=1}^n\frac{1}{\sum_{i=1}^kc_{\sigma[i]}(\sigma)}$$

This volume in the event space is the probability $P\left[I(\sigma)\wedge\sum_{i=1}^nc_i(\sigma)X_i<1\right]$, so we can recombine the permutations to get the probability we have not reached the required distance in $n$ steps:

$$P[D_n<1]=\sum_{\sigma\in S_n}\frac{1}{n!\prod_{k=1}^n\sum_{i=1}^kc_{\sigma[i]}(\sigma)}$$

Now, to determine the expected value of the number of steps needed $N$, we observe that $P[D_n<1]=P[N>n]$. Then

\begin{align*} E[N]&=\sum_{n=1}^\infty nP[N=n] \\ &=1-P[N>1]+\sum_{n=2}^\infty n(P[N>n-1]-P[N>n]) \\ &=1+\sum_{n=1}^\infty P(N>n) \\ &=\boxed{1+\sum_{n=1}^\infty\sum_{\sigma\in S_n}\frac{1}{n!\prod_{k=1}^n\sum_{i=1}^kc_{\sigma[i]}(\sigma)}} \end{align*}

---

Now, I highly doubt that this expression can be reduced to a closed form, but the infinite sum seems to converge fairly quickly. Implementing this expression in MATLAB, I was able to calculate the first 10 terms within a minute, but it seems to run in at least $O(n!)$ time. The results are shown below:

Max n  E[N]
1      2.00000
2      2.75000
3      3.16667
4      3.34896
5      3.41458
6      3.43464
7      3.43996
8      3.44120
9      3.44146
10     3.44151

From this sequence, we can conclude that the true expected value is probably around $E[N]=3.4416$.

Simulating the process for $10^6$ iterations with 10 trials, we get the sorted results:

E[N]
3.4389
3.4399
3.4410
3.4412
3.4416
3.4417
3.4418
3.4420
3.4421
3.4424

which do appear to be centered around the calculated value of 3.4416.

---

Code:

%% Formula Calculation
N = 10;
P = zeros(N,1);
for n=1:N
  p = perms(1:n);
c = zeros(size(p));
  c(:,1)=1;
  for ii = 2:n
    c_plus = zeros(size(p));
    c_plus(:,ii) = 1;
    c_plus(:,ii-1) = -1;
    is_pos = find(p'==ii)<find(p'==ii-1);
    c = c+(2is_pos-1).c_plus;
  end
r = zeros(size(p));
  for ii = 1:size(p,1)
    r(ii,:) = c(ii,p(ii,:));
  end
P(n) = sum(1./prod(cumsum(r,2),2))/factorial(n);
end
E = 1+cumsum(P);
E(end)

%% Simulation

E_est = zeros(10,1);
N_iter=1e6;
for ii = 1:10
  X = [zeros(N_iter,1), rand(N_iter,30)];
  D = abs(diff(X,1,2));
  D_tot = cumsum(D,2);
  Y = sum(D_tot<1,2)+1;
  E_est(ii) = mean(Y);
end
sort(E_est)

Answer 3

There's a fast way to approximate the constant you're looking for that I don't have time to implement.

The size of the next step depends on your current location, so we should consider the probability density of location and how far you've already traveled. Call the joint distribution after $n$ steps $f_n(d, x)$. Then we can find a recursive formula

$$f_{n+1}(d, x) = \int_0^1 f_n(d-|x-y|, y) dy$$

Where $f_n(d, x)$ is considered to be $0$ for $d < 0$. The only region we're interested in is $x, y \in [0,1]$, which we can split up into four quarters along the diagonals of the square. By doing that, the formula for the $d< x, 1-x$ quadrant only depends on the same quadrant the previous step. The formula is the sum of two integrals with bounds that are linear in $x, d$, with no $\max$ or $\min$ required. Similarly the other three quadrants can be written as a sum of two integrals over the previous quadrants without $\max$ or $\min$. And after step $1$, the formula for the density in each quadrants will be a constant hence, in particular, a polynomial in $x, d$, which means after each step all the quadrants will still have polynomial densities. And polynomial multiplication, substitution and integration is fast, so this way you can work out hundreds of terms very quickly.

The way you recover your probability $P(N \le n)$ from this is as the integral of $f_n(d, x)$ over all four quadrants. To implement this you just want a computer algebra package that can integrate and substitute two-variable polynomials. The distribution is also symmetric so you can ignore one of the middle quadrants to simplify things a little.

Answer 4

Not really an answer.

Let $Y_i = |X_{i}-X_{i-1}|$.

The density is given by $f_Y(y) = 2(1-y) [0\le y \le 1]$, with $E[Y]=1/3$ and $\sigma^2_Y=1/18$.

Furthermore, the variables $(Y_i, Y_{i+1})$ are not independent - but the correlation factor is quite low: $\rho = 1/10$. And the other pairs $(Y_i, Y_{i+k})$ are independent for $k>1$.

As a rough approximation, if we assume $Y_i$ are all independent, we get the recursion

$$g(x) = 1 + 2 \int_0^1 g(x-y) 2 (1-y) dy = 1 - 2 \int_{0}^x g(u) (x-1-u) du \tag 1$$

for $0<x<1$, <with $g(x)=0$ for $x<0$, where $g(x)$ is the expected number of tries until the sum reaches $x$:

$$g(x) = \min\left\{n\in\mathbb{N}: \sum\limits_{i=1}^{n}{Y_i}\geq x\right\} $$

Now, $(1)$ is an integral equation that can be solved by Laplace transform. The solution is

$$g(x) = e^x (\sin x + \cos x)$$

with $g(1) = 3.756049227$

This is not far from the true value ($\approx 3.82$). The discrepancy, of course, should be due to the neglected dependency. Furthermore, given that the correlation is positive, and that $P(Y < E[Y]) = 0.555 > \frac12$, it seems reasonable that the correlation tends to decrease the sum, and hence to increase slightly the expected number of iterations.