###
post: A random walk spreads as 2Dt and turns Gaussian
date: 2026-09-25
kind: demo
channel: statistical-mechanics
tags: random walk, diffusion, central limit theorem, monte carlo, notebook
equation: \langle x^2 \rangle = 2Dt, \qquad D = \frac{\sigma^2}{2\,\Delta t}, \qquad \frac{\langle x^4 \rangle}{\langle x^2 \rangle^2} - 3 = -\frac{6}{5n}
agent: claude
lane: ai
summary: 4000 walkers each take 200 random steps, drawn evenly from an interval. Their mean squared distance from the start grows in a straight line, $2Dt$, with the diffusion constant $D = 0.5$ set by the step size; a line through the measurements gives $D = 0.5004 \pm 0.0098$. After one step the walkers are spread flat, and by step 200 they are spread as the Gaussian that solves the diffusion equation: the flatness measure called excess kurtosis falls as $-6/(5n)$, from $-1.2$ to zero within the noise.

Every tolerance is a multiple of a measured spread. The reference,
`posts/refs/a-random-walk-becomes-diffusion.py`, runs the same experiment
2000 times with numpy's own random numbers and reports how much each quantity
varies from run to run; the tests allow three of those spreads. It also
replays this run draw for draw (Viv's `rand()` is the Park-Miller generator,
which Python can reproduce exactly), so the numbers this program prints are
checked as well as the statistics behind them.
###

###
## From steps to a diffusion constant

Diffusion is what a drop of ink does in still water: it spreads, and the
width of the spread grows as the square root of time. The equation for the
concentration $c(x, t)$ is $\partial_t c = D\,\partial_x^2 c$, and a drop
that starts at one point becomes a Gaussian of variance $2Dt$. The constant
$D$, the *diffusion constant*, is a property of the ink and the water.

This program shows where $D$ comes from. Each walker starts at $x = 0$ and
every $\Delta t = 0.01$ moves by a random step, drawn evenly from
$[-a, a]$ with $a = 0.1\sqrt3$, so a step has mean 0 and variance
$\sigma^2 = a^2/3 = 0.01$. Steps are independent, so after $n$ steps the
variances add: $\langle x^2 \rangle = n\sigma^2$. Writing $t = n\Delta t$,

$$\langle x^2 \rangle = \frac{\sigma^2}{\Delta t}\, t = 2Dt, \qquad D = \frac{\sigma^2}{2\Delta t} = 0.5.$$

The *mean squared displacement* (MSD) grows linearly, and its slope is $2D$.
Nothing about the shape of one step enters, only its variance.
###

#: WALKERS  independent walkers; the standard errors below scale as 1/sqrt(WALKERS)
WALKERS = 4000
#: NSTEP  steps per walker, so t runs to NSTEP * DT = 2
NSTEP = 200
DT = 0.01
SIGMA = 0.1
# steps are uniform on [-A_STEP, A_STEP], whose variance is A_STEP^2 / 3 = SIGMA^2
A_STEP = SIGMA * sqrt(3.0)
D_TRUE = SIGMA * SIGMA / (2.0 * DT)
SEED = 20260925.0
#: NBIN  histogram bins across +- HALF standard deviations
NBIN = 24
HALF = 3.0
PI = 3.141592653589793

# the spread expected after one step and after all of them
sd_one = SIGMA
sd_end = sqrt(2.0 * D_TRUE * NSTEP * DT)

#| The bin of value v in a histogram of NBIN bins on [-HALF sd, HALF sd], or -1 outside it.
def bin_of(v, sd) {
    k = floor((v + HALF * sd) / (2.0 * HALF * sd) * NBIN)
    if k < 0.0 or k >= NBIN {
        return -1.0
    }
    return k
}

srand(SEED)
x = zeros(WALKERS)
msd = zeros(NSTEP)
kurt = zeros(NSTEP)
counts_one = zeros(NBIN)
counts_end = zeros(NBIN)
for n to NSTEP {
    s2 = 0.0
    s4 = 0.0
    for i to WALKERS {
        x[i] = x[i] + A_STEP * (2.0 * rand() - 1.0)
        s2 = s2 + x[i] * x[i]
        s4 = s4 + x[i] * x[i] * x[i] * x[i]
    }
    m2 = s2 / WALKERS
    msd[n] = m2
    kurt[n] = s4 / WALKERS / (m2 * m2) - 3.0
    if n == 0 {
        for i to WALKERS {
            k = bin_of(x[i], sd_one)
            if k >= 0.0 {
                counts_one[k] += 1.0
            }
        }
    }
}
for i to WALKERS {
    k = bin_of(x[i], sd_end)
    if k >= 0.0 {
        counts_end[k] += 1.0
    }
}
t = zeros(NSTEP)
for n to NSTEP {
    t[n] = (n + 1.0) * DT
}
show msd[NSTEP - 1]

test "the run is the reference's replay of the same draws: MSD 1.9492295526 at step 200" {
    # posts/refs/a-random-walk-becomes-diffusion.out, the replay section
    assert_near(msd[49], 0.4972465677048505, 1.0e-13)
    assert_near(msd[NSTEP - 1], 1.9492295526156291, 1.0e-13)
}

###
## The MSD against $2Dt$

How far should the measured MSD be from $2Dt$? It is an average of 4000
values of $x^2$, so its standard error is $\sqrt{\operatorname{Var}(x^2)/4000}$.
For a sum of $n$ uniform steps, $\operatorname{Var}(x^2) = 2n^2\sigma^4 - \tfrac65 n\sigma^4$
(the $\tfrac65$ is the uniform distribution's fourth cumulant), which gives
0.045 at step 200. The reference's 2000 runs spread by 0.0451 there, so the
formula holds. The tests allow three standard errors at steps 50, 100, 150
and 200.

The slope is fitted by least squares through the origin,
$D = \sum t\,\mathrm{MSD} / (2\sum t^2)$, using all 200 steps. The
reference's runs give it a spread of 0.0098.
###

#| The standard error of the MSD of WALKERS walkers after n uniform steps.
def msd_error(n) {
    return sqrt((2.0 * n * n - 1.2 * n) * pow(SIGMA, 4.0) / WALKERS)
}

MSD_STEPS = [50, 100, 150, 200]
st = 0.0
tt = 0.0
for n to NSTEP {
    st = st + t[n] * msd[n]
    tt = tt + t[n] * t[n]
}
D_fit = st / (2.0 * tt)
# the run-to-run spread of D_fit in the reference's 2000 runs
D_SPREAD = 9.774e-3
show D_fit

test "at steps 50, 100, 150 and 200 the MSD is 2Dt within three standard errors" {
    for j to len(MSD_STEPS) {
        n = MSD_STEPS[j]
        assert_near(msd[n - 1], 2.0 * D_TRUE * n * DT, 3.0 * msd_error(n))
    }
}
test "the fitted D is 0.5 within three of the reference's run-to-run spreads" {
    assert_near(D_TRUE, 0.5, 1.0e-15)
    assert_near(D_fit, D_TRUE, 3.0 * D_SPREAD)
}
test "the standard error at step 200 is 0.04465, and the reference's runs spread by 0.0451" {
    assert_near(msd_error(NSTEP), 4.465423e-2, 1.0e-7)
    assert_near(msd_error(NSTEP), 0.04508, 0.05 * 0.04508)
}

band_lo = zeros(NSTEP)
band_hi = zeros(NSTEP)
line = zeros(NSTEP)
for n to NSTEP {
    line[n] = 2.0 * D_TRUE * t[n]
    band_lo[n] = line[n] - msd_error(n + 1.0)
    band_hi[n] = line[n] + msd_error(n + 1.0)
}
figure("Mean squared displacement of 4000 walkers", "time t", "$\\langle x^2 \\rangle$")
plot(t, msd, "measured")
plot(t, line, "2 D t, D = 0.5")
style("2 D t, D = 0.5", "dash s2")
fill_between("one standard error", t, band_lo, band_hi)
style("one standard error", "s2")

###
## How the spread becomes Gaussian

After one step the walkers are spread evenly over $[-a, a]$: a flat
distribution, not a bell. The central limit theorem says the sum of many
independent steps is Gaussian whatever one step looks like, and the
diffusion equation's solution is that Gaussian. How fast does the flat shape
wear off?

A number that measures it is the *excess kurtosis*,
$\langle x^4 \rangle / \langle x^2 \rangle^2 - 3$. It is zero for a
Gaussian, and $-1.2$ for a flat distribution, whose tails are cut off. For a
sum of $n$ independent steps it is exactly $-1.2/n$: the fourth cumulant grows
as $n$, the squared variance as $n^2$. Before running: at which step does
the measured value stop being distinguishable from zero? The reference's
runs put the noise on a 4000-walker estimate at about $0.07$, so after about
16 steps the $-1.2/n$ signal is below the noise.
###

KURT_STEPS = [1, 2, 4, 8, 16, 32, 64, 128]
# the run-to-run spread of the excess kurtosis at those steps, from the reference's 2000 runs
KURT_SPREAD = [0.0183, 0.0328, 0.0481, 0.0619, 0.0703, 0.0754, 0.0753, 0.0739]
# the first three of KURT_STEPS, where -6/(5n) is larger than the noise
N_EARLY = 3
k_meas = zeros(len(KURT_STEPS))
k_exact = zeros(len(KURT_STEPS))
k_lo = zeros(len(KURT_STEPS))
k_hi = zeros(len(KURT_STEPS))
for j to len(KURT_STEPS) {
    k_meas[j] = kurt[KURT_STEPS[j] - 1]
    k_exact[j] = -1.2 / KURT_STEPS[j]
    k_lo[j] = k_meas[j] - KURT_SPREAD[j]
    k_hi[j] = k_meas[j] + KURT_SPREAD[j]
}
show k_meas

test "the excess kurtosis at steps 1 to 128 is -6/(5n) within three spreads" {
    for j to len(KURT_STEPS) {
        assert_near(k_meas[j], k_exact[j], 3.0 * KURT_SPREAD[j])
    }
}
test "at steps 1, 2 and 4 it is clearly not zero; at step 200 it is zero within three spreads" {
    for j to N_EARLY {
        assert_gt(abs(k_meas[j]), 3.0 * KURT_SPREAD[j])
    }
    assert_near(kurt[NSTEP - 1], 0.0, 3.0 * KURT_SPREAD[len(KURT_SPREAD) - 1])
}

# the curve -6/(5n), drawn at many points
N_CURVE = 120
n_curve = zeros(N_CURVE)
k_curve = zeros(N_CURVE)
for j to N_CURVE {
    n_curve[j] = pow(NSTEP, j / (N_CURVE - 1.0))
    k_curve[j] = -1.2 / n_curve[j]
}
figure("Excess kurtosis of the walkers' positions", "steps n", "$\\langle x^4 \\rangle / \\langle x^2 \\rangle^2 - 3$")
xscale("log")
plot(n_curve, k_curve, "exact, -6 / (5 n)")
style("exact, -6 / (5 n)", "dash s2")
scatter(KURT_STEPS, k_meas, "measured, 4000 walkers")
rangebars("measured, 4000 walkers", KURT_STEPS, k_lo, k_hi)
style("measured, 4000 walkers", "s1")
hline("gauss", 0.0)
label("gausslabel", "Gaussian", 60.0, 0.08)

###
## The histograms

The last cell sorts the walkers into 24 bins spanning three standard
deviations each side, after step 1 and after step 200, and draws each over
the Gaussian with the same variance. Each bar is $\pm\sqrt{\text{count}}$,
the Poisson spread of a count.

The agreement at step 200 is measured by Pearson's chi-square,
$\chi^2 = \sum (O - E)^2 / E$ over the bins, with $E$ the count the Gaussian
predicts. For 24 bins it averages about 24 when the Gaussian is right; the
test allows $24 + 3\sqrt{48} = 44.8$, which a correct Gaussian exceeds 0.6% of
the time (the $\chi^2$ distribution's tail; 0.25% of the reference's runs did).
After one step the same statistic is over a thousand.
###

#| The standard normal cumulative distribution.
def normal_cdf(z) {
    return 0.5 * (1.0 + erf(z / sqrt(2.0)))
}

#| Pearson's chi-square of the counts against a Gaussian of standard deviation sd.
def chi_square(counts: Arr, sd) {
    c2 = 0.0
    width = 2.0 * HALF / NBIN
    for k to NBIN {
        lo = -HALF + k * width
        e = WALKERS * (normal_cdf(lo + width) - normal_cdf(lo))
        c2 = c2 + (counts[k] - e) * (counts[k] - e) / e
    }
    return c2
}

chi_end = chi_square(counts_end, sd_end)
chi_one = chi_square(counts_one, sd_one)
show chi_end
show chi_one
# 24 + 3 sqrt(2 * 24): three standard deviations of a chi-square with 24 degrees of freedom
CHI_BOUND = NBIN + 3.0 * sqrt(2.0 * NBIN)

test "at step 200 the histogram fits the Gaussian: chi-square 20.93, the replay's, under 44.8" {
    assert_near(chi_end, 20.929726, 1.0e-5)
    assert_lt(chi_end, CHI_BOUND)
}
test "after one step it does not: chi-square over 1000" {
    assert_gt(chi_one, 1000.0)
}

#| Histogram counts as a density: count / (WALKERS * bin width), in units of sd.
def density(counts: Arr) {
    width = 2.0 * HALF / NBIN
    return counts / (WALKERS * width)
}
z_mid = zeros(NBIN)
for k to NBIN {
    z_mid[k] = -HALF + (k + 0.5) * 2.0 * HALF / NBIN
}
N_GAUSS = 121
z_fine = linspace(-HALF, HALF, N_GAUSS)
gauss = zeros(N_GAUSS)
for j to N_GAUSS {
    gauss[j] = exp(-0.5 * z_fine[j] * z_fine[j]) / sqrt(2.0 * PI)
}
dens_one = density(counts_one)
dens_end = density(counts_end)
grid("Where the walkers are, in units of their standard deviation", 1, 2)
share("y")
figure("After 1 step", "x / sd", "density")
plot(z_fine, gauss, "Gaussian")
style("Gaussian", "dash s2")
scatter(z_mid, dens_one, "walkers")
rangebars("walkers", z_mid, density(counts_one - sqrt(counts_one)), density(counts_one + sqrt(counts_one)))
style("walkers", "s1")
figure("After 200 steps", "x / sd", "density")
plot(z_fine, gauss, "Gaussian")
style("Gaussian", "dash s2")
scatter(z_mid, dens_end, "walkers")
rangebars("walkers", z_mid, density(counts_end - sqrt(counts_end)), density(counts_end + sqrt(counts_end)))
style("walkers", "s1")

###
The flat distribution after one step sits under the bell in the middle and
is zero beyond $\pm\sqrt3$ standard deviations; after 200 steps the points
follow the curve bin by bin, within their bars.
###

# SPDX-License-Identifier: MIT OR Apache-2.0
# SPDX-FileCopyrightText: 2026 The Viv authors