Notebooks · posts/a-random-walk-becomes-diffusion

A random walk spreads as 2Dt and turns Gaussian

Adam Zettel · September 25, 2026 · written with an agent

4000 walkers each take 200 random steps, drawn evenly from an interval. Their mean squared distance from the start grows in a straight line, $2 D t$ , 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/ (5 n)$ , from $- 1.2$ to zero within the noise.

⟨ x^{2} ⟩ = 2 D t, D = \frac{σ ^{2}}{2 Δ t}, \frac{⟨ x ^{4} ⟩}{⟨ x ^{2} ⟩ ^{2}} - 3 = - \frac{6}{5 n}

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.

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

Shift+Enter in a cell runs it and moves to the next. Running a cell runs every cell above it too, as one program.

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 $2 D t$ . 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 $Δ t = 0.01$ moves by a random step, drawn evenly from $[- a, a]$ with $a = 0.1 3$ , so a step has mean 0 and variance $σ^{2} = a^{2} /3 = 0.01$ . Steps are independent, so after $n$ steps the variances add: $⟨ x^{2} ⟩ = n σ^{2}$ . Writing $t = n Δ t$ ,

⟨ x^{2} ⟩ = \frac{σ ^{2}}{Δ t} t = 2 D t, D = \frac{σ ^{2}}{2Δ t} = 0.5.

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

[1]lines 42–112

#: 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)
}

#: 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 $2 D t$

How far should the measured MSD be from $2 D t$ ? It is an average of 4000 values of $x^{2}$ , so its standard error is $Var (x^{2}) /4000$ . For a sum of $n$ uniform steps, $Var (x^{2}) = 2 n^{2} σ^{4} - \frac{6}{5} n σ^{4}$ (the $\frac{6}{5}$ 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 MSD / (2 \sum t^{2})$ , using all 200 steps. The reference's runs give it a spread of 0.0098.

[2]lines 130–175

#| 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")

#| 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, $⟨ x^{4} ⟩ / ⟨ x^{2} ⟩^{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.

[3]lines 196–241

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)

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 count$ , the Poisson spread of a count.

The agreement at step 200 is measured by Pearson's chi-square, $χ^{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 + 348 = 44.8$ , which a correct Gaussian exceeds 0.6% of the time (the $χ^{2}$ distribution's tail; 0.25% of the reference's runs did). After one step the same statistic is over a thousand.

[4]lines 259–321

#| 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 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 3$ standard deviations; after 200 steps the points follow the curve bin by bin, within their bars.

[5]lines 329–330

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

From steps to a diffusion constant

The MSD against 2Dt

How the spread becomes Gaussian

The histograms

The MSD against $2 D t$