# 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? – Chris Culter Jan 8 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$. – Misha Lavrov Jan 8 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. – Brevan Ellefsen Jan 8 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. – Chris Culter Jan 8 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)}.$$

• Thanks! For those of us who aren't familiar with the digamma function, do you think you could explain how it's being used a little more? It would be nice to highlight which steps are nontrivial, where the proof of the identity would require some interesting transformation like integration by parts, versus which steps are just routine bookkeeping of sums and manipulating definitions. Also, it would help to explain where the gamma identity comes from: is that one nontrivial? – Chris Culter Jan 11 at 20:25
• 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. – Maxim Jan 12 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$$.

This is not a complete answer - it doesn't explain why the constant errors cancel out. But it's a start.