average length of "harmonic walk" mod $n$ is $\ln(n)$ -- why?

Let a harmonic walk be a stochastic process where at turn $$t \ge 2$$ the probability of ending the walk is $$\frac{1}{t}$$ and the first turn never ends the walk. A harmonic walk can be viewed as a distribution over $$\mathbb{N}_{\ge 1}$$ without finite mean.

Empirically, it seems to be the case that the average length of a harmonic walk mod $$n$$ is $$\ln\left(n\right)$$ . I'm wondering why this should be the case. The definition of the harmonic series (1) looks similar to the integral definition of the natural log, $$\ln(x) \stackrel{\text{def}}{=} \int_{1}^{x} 1/s \;\text{ds}$$, so there might be reason to suspect that a relationship might exist, but that resemblance isn't remotely convincing.

The harmonic series $$H$$ is given below (1) with the symbol "$$H$$" referring to the formal sum, not the value.

$$H \stackrel{\text{def}}{=} \sum_{k=1}^{\infty} \frac{1}{k} \;\;\;\;\;\;\;\;\text{and H diverges} \tag{1}$$

We can write an equivalent formal sum of partial products $$H'$$ (2a). I don't know the exact term for the relationship between $$H$$ and $$H'$$ as formal expressions, but $$\frac{1}{k}$$ is finite and definitely equal to $$\prod_{l=2}^{k} \frac{l-1}{l}$$.

$$H' \stackrel{\text{def}}{=} \sum_{k=1}^{\infty} \prod_{l=2}^{k} \frac{l-1}{l} \tag{2a}$$ $$H' = 1 + \left(\frac{1}{2}\right) + \left(\frac{1}{2} \times \frac{2}{3} \right) \dots \tag{2b}$$

Squinting at the definition of $$H'$$, we can come up with a stochastic process $$\Psi$$ .

• At every turn in $$\Psi$$, we either stop immediately or add 1 to $$t$$ .

• If $$t = 1$$, then advance to the next turn with probability $$1$$.

• If $$t \ne 1$$, then stop with probability $$\frac{1}{t}$$ and continue to the next turn with probability $$\frac{t-1}{t}$$ .

If we take the output of $$\Psi \;\;\text{mod}\;\; n$$ for various values of $$n$$, it seems to match $$\ln\left(n\right)$$ .

n        Ψ mod n     ln(n)
2          0.693     0.693
3          1.099     1.099
4          1.386     1.386
5          1.609     1.609
6          1.793     1.792
7          1.947     1.946


Why should this be the case?

Here is the python code used to produce samples:

import random
import math
import os

STOP = "STOP"
GO   = "GO"

def step(n):
assert (n > 1)
if random.randint(1, n) == 1:
return STOP
else:
return GO

def walk_length():
n = 2
len_ = 1
while GO == step(n):
len_ += 1
n += 1
return len_

for x in range(1000000):
print(walk_length() % int(os.getenv("MODULUS")))


and the summary script from here:

#! /bin/sh

Rscript -e 'summary (as.numeric (readLines ("stdin")));'


A sample run.

> env MODULUS=7 python process_steps.py | summary
Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
0.000   1.000   1.000   1.947   3.000   6.000

• Could you please use a different letter for the turn number and the modulus? Jan 8, 2019 at 21:43
• @BrevanEllefsen Aside from the fact that it takes a bit of work to get to the harmonic numbers from here, they alone don't explain why the expectation here is much closer to $\ln n$ than to $\ln n + \gamma$. Jan 8, 2019 at 22:03
• @MishaLavrov I know, and your answer fleshes this out much better than my comment. I'm not quite sure either why the constant cancels, but I thought that fact might help the OP a bit. Now that your answer supersedes my comment entirely, I'll remove it. Jan 8, 2019 at 22:06
• It feels like you'll want to introduce a new function of $s=t+x$ with $x\in[0,1)$, simplify its derivative, and then integrate the derivative over $x$, converting the integral over $x$ of a sum over $t$ into a single integral over $s$, and ending up with $\int_1^n\frac1s\,ds$. But I'm not sure of the details, to put it mildly. Jan 8, 2019 at 23:02

We have $$\operatorname P[L = n] = \frac 1 2 \frac 2 3 \cdots \frac {n - 1} n \frac 1 {n + 1} = \frac 1 {n (n + 1)} [n \geq 1], \\ \operatorname P[(L \bmod m) = n] = \sum_{k \geq 0} \operatorname P[L = m k + n] \,[0 \leq n < m], \\ \operatorname E[L \bmod m] = \sum_{0 < n < m} \sum_{k \geq 0} \frac n {(m k + n) (m k + n + 1)} = \\ \sum_{0 < n < m} \frac n m \left( \psi {\left( \frac {n + 1} m \right)} - \psi {\left( \frac n m \right)} \right) = \\ \psi(1) - \frac 1 m \sum_{0 < n \leq m} \psi {\left( \frac n m \right)} = \ln m.$$ The last identity follows from $$\Gamma(m z) = m^{m z - 1/2} (2 \pi)^{(1 - m)/2} \prod_{0 < n \leq m} \Gamma {\left( z + \frac {n - 1} m \right)}.$$
• The starting point is the Weierstrass definition of the gamma function as an infinite product. Taking the logarithm and differentiating gives $$\psi(z) = -\frac 1 z - \gamma - \sum_{k \geq 1} \left( \frac 1 {k + z} - \frac 1 k \right).$$ Therefore we get $\sum_{k \geq 0} (1/(k + a) - 1/(k + b)) = \psi(b) - \psi(a)$. The $\Gamma(m z)$ formula is called the multiplication-theorem of Gauss and Legendre in Whittaker and Watson. Taking the logarithm, differentiating and letting $z = 1/m$ gives the desired identity. Jan 12, 2019 at 0:57
If the $$k^{\text{th}}$$ turn has a probability of $$\frac1k$$ of ending the walk, then the overall probability that the walk ends on turn $$k$$ and not before is $$\frac12 \cdot \frac23 \cdots \frac{k-2}{k-1} \cdot \frac1k = \frac{1}{k(k-1)}$$. So the probability that the walk ends on a turn that's $$k$$ modulo $$n$$ is $$\frac{1}{k(k-1)} + \frac{1}{(n+k)(n+k-1)} + \frac{1}{(2n+k)(2n+k-1)} + \dots$$ for $$k=2,\dots,n-1$$, dropping the first term for $$k=0$$ and $$k=1$$.
We have $$\frac{1}{(jn+k)(jn+k-1)} \le \frac{1}{j^2(n-1)^2}$$ so the sum of all the terms except the first can be upper-bounded by $$\frac{1}{(n-1)^2}(\frac11 + \frac14 + \frac19 + \dots) = \frac{\pi^2}{6(n-1)^2}$$. So if $$X$$ is the value of the walk mod $$n$$ then its expected value satisfies $$\sum_{k=2}^{n-1} k \cdot \frac{1}{k(k-1)} \le \mathbb E[X] \le \sum_{k=2}^{n-1} k \cdot \frac{1}{k(k-1)} + \sum_{k=0}^{n-1} k \cdot \frac{\pi^2}{6(n-1)^2}.$$ The lower bound simplifies to the $$(n-2)^{\text{th}}$$ harmonic number $$H_{n-2}$$, which is well-approximated by $$\ln n$$ up to a constant error. The error term in the upper bound simplifies to $$\binom n2 \frac{\pi^2}{6(n-1)^2}$$, which approaches $$\frac{\pi^2}{12}$$ as $$n \to \infty$$.