Random Walk Simulation

What the simulation shows

Every walker starts at the origin and takes n independent steps. On the line a step is +1 with probability p and −1 otherwise; on the square lattice it goes north, south, east or west with probability ¼ each. The four panels answer four different questions about the same walkers:

• Paths. Walker no. 1 (orange) wanders, crosses zero and stays away for long stretches, looking nothing like the average. The dashed envelopes are the theoretical mean ±1 and ±2 standard deviations; roughly 68% and 95% of walkers sit inside them at any moment.
• Distribution. In 1-D the final position is 2B − n with B binomial, so the histogram should match the exact binomial expectation drawn over it; for large n that is a normal bell. In 2-D the distance from the origin follows, approximately, a Rayleigh distribution.
• Spread. The root-mean-square (RMS) distance tracks √t: after 100 steps a typical walker is only about 10 steps from home, after 10,000 steps about 100.
• Returns. The fraction of walkers that have revisited the origin, against the exact probability.

Things to try: switch to the 2-D lattice and watch the returns curve rise far more slowly; set p = 0.6 and see the cloud drift, the envelopes tilt and the returns level off at 80%; drop to a single walker to see how little one path says about the ensemble.

Why the spread grows like √n

Write the 1-D position as a sum of independent steps, X_n = ξ₁ + … + ξ_n, with each ξ = +1 with probability p and −1 with probability q = 1 − p. Each step has mean p − q and variance 4pq, and variances of independent steps add:

For the unbiased walk (p = ½) the mean is 0 and E[X_n²] = n, so the RMS displacement is exactly √n. The 2-D lattice walk gives the same E[R_n²] = n, because each unit step has mean zero and is independent of the current position: the cross term 2R·ξ averages to zero and every step adds exactly 1 to E[R²]. A bias adds a drift that grows like n and soon dominates the √n spread, which is why a small edge in a game of chance beats luck in the long run.

From walk to diffusion. Shrink the step to Δx and the time per step to Δt, keeping D = Δx²/(2Δt) fixed. The central limit theorem turns the binomial into a Gaussian with variance 2Dt, which is the solution of the diffusion equation ∂u/∂t = D ∂²u/∂x². Donsker's theorem says the whole rescaled path converges to Brownian motion. The diffusion simulation runs walkers and the partial differential equation side by side, and the Brownian motion simulation starts from the continuum limit.

Coming back: Pólya’s recurrence theorem

In 1921 George Pólya proved that a simple random walk returns to its starting point with probability 1 on the line and on the plane lattice, but not in three or more dimensions. In 3-D the probability of ever returning is only about 0.34.

The lab computes the exact probability of a return by step t. Let u_k be the probability of being at the origin at step 2k and f_k the probability that the first return happens then. A visit at step 2k is a first return at some step 2j followed by a return in the remaining 2(k − j) steps:

Solving for f and summing gives the dashed curve in the returns panel. For large k, u_k falls like 1/√(πk) in 1-D and 1/(πk) in 2-D. Both sums diverge, which is the criterion for recurrence; in 3-D u_k falls like k^−3/2, the sum converges and the walk is transient. The divergence in 2-D is only logarithmic, so returns are certain but slow: after 1,000 steps about 68% of 2-D walkers have come back, after 10,000 steps still only about 74%, while in 1-D the figures are 97% and 99%. Even in 1-D the expected time to return is infinite. With a bias, the 1-D return probability drops to 1 − |p − q|.

Algorithm and implementation

1. Set all W positions to the origin and mark every walker as “not yet returned”.
2. For each step, for each walker: draw u ~ U[0, 1). In 1-D move +1 if u < p, else −1; in 2-D pick one of four directions from ⌊4u⌋.
3. If the walker is at the origin and has not been there since the start, count a return. Track the largest distance reached.
4. Record the mean position, the mean squared distance and the returned fraction at each step.
5. Compare with theory: n(2p − 1), √E[R²] and the renewal-equation return probability.

Positions are integers, so the sums are exact and the JavaScript and Python versions agree to the last digit for a given seed. The renewal sum costs about n²/8 multiplications.

Example result

With 200 unbiased walkers taking 1,000 steps each (seed 42), the mean final position is 0.77 against a theoretical 0.00, and the RMS displacement is 33.36 against √1000 = 31.62. 96.5% of the walkers returned to the origin at least once, compared with the exact probability 97.5%. The walker that strayed farthest reached 96.0 steps from home.

Implementation

The page runs simulation.js from the package below; the Python file is a line-by-line port with the same random-number generator, so a given seed produces the same result in both languages. The excerpts show the core of each.

// Advance every walker by one step. 1-D: +1 with probability p, else −1.
// 2-D: move to one of the four lattice neighbours, each with probability 1/4.
export function step(s, p, rng) {
  let sumX = 0, sumR2 = 0;
  for (let i = 0; i < s.n; i++) {
    if (s.dims === 1) {
      s.x[i] += rng.uniform() < p ? 1 : -1;
    } else {
      const k = rng.int(4);
      if (k === 0) s.x[i]++; else if (k === 1) s.x[i]--; else if (k === 2) s.y[i]++; else s.y[i]--;
    }
    const r2 = s.x[i] * s.x[i] + s.y[i] * s.y[i];
    if (r2 === 0 && !s.back[i]) { s.back[i] = 1; s.returned++; } // first return to the origin
    if (r2 > s.maxR2) s.maxR2 = r2;
    sumX += s.x[i];
    sumR2 += r2;
  }
  s.t++;
  s.sumX = sumX;
  s.sumR2 = sumR2;
}

// Exact P(first return to the origin by step 2k), k = 0 … ⌊steps/2⌋.
// u[k] = P(at origin at step 2k): C(2k,k)(pq)^k in 1-D, (C(2k,k)/4^k)² in 2-D.
// First returns f[k] follow from the renewal equation u[k] = Σ_j f[j]·u[k−j].
export function returnCurve(steps, dims, p) {
  const m = Math.floor(steps / 2);
  const u = new Float64Array(m + 1), f = new Float64Array(m + 1), F = new Float64Array(m + 1);
  const q4 = 4 * p * (1 - p);
  let a = 1, c = 1; // a = C(2k,k)/4^k, c = (4pq)^k
  u[0] = 1;
  for (let k = 1; k <= m; k++) {
    a = a * (2 * k - 1) / (2 * k);
    c = c * q4;
    u[k] = dims === 1 ? a * c : a * a;
    let conv = 0;
    for (let j = 1; j < k; j++) conv += f[j] * u[k - j];
    f[k] = u[k] - conv;
    F[k] = F[k - 1] + f[k];
  }
  return F;
}

def step(s, p, rng):
    """Advance every walker one step: 1-D +1 w.p. p else -1; 2-D one of four neighbours."""
    x, y, back = s["x"], s["y"], s["back"]
    sum_x = sum_r2 = 0
    for i in range(s["n"]):
        if s["dims"] == 1:
            x[i] += 1 if rng.uniform() < p else -1
        else:
            k = rng.int(4)
            if k == 0:
                x[i] += 1
            elif k == 1:
                x[i] -= 1
            elif k == 2:
                y[i] += 1
            else:
                y[i] -= 1
        r2 = x[i] * x[i] + y[i] * y[i]
        if r2 == 0 and not back[i]:  # first return to the origin
            back[i] = True
            s["returned"] += 1
        if r2 > s["max_r2"]:
            s["max_r2"] = r2
        sum_x += x[i]
        sum_r2 += r2
    s["t"] += 1
    s["sum_x"], s["sum_r2"] = sum_x, sum_r2


def return_curve(steps, dims, p):
    """Exact P(first return by step 2k) from the renewal equation u[k] = sum_j f[j] u[k-j]."""
    m = steps // 2
    u, f, F = [0.0] * (m + 1), [0.0] * (m + 1), [0.0] * (m + 1)
    q4 = 4 * p * (1 - p)
    a = c = 1.0  # a = C(2k,k)/4^k, c = (4pq)^k
    u[0] = 1.0
    for k in range(1, m + 1):
        a = a * (2 * k - 1) / (2 * k)
        c = c * q4
        u[k] = a * c if dims == 1 else a * a
        conv = 0.0
        for j in range(1, k):
            conv += f[j] * u[k - j]
        f[k] = u[k] - conv
        F[k] = F[k - 1] + f[k]
    return F

Run the full model from a terminal:

node simulation.js --mode 1d --walkers 200 --seed 42
python simulation.py --mode 1d --walkers 200 --seed 42

Download the model

Assumptions & limitations

• Independent, identical steps. Each step ignores the past (a Markov chain with no memory of direction). Walks that avoid their own trail, persist in one direction or interact with each other behave differently; self-avoiding polymers, for example, spread faster than √n.
• Unit steps on a lattice, discrete time, no interaction. Real particles move continuously and collide; the lattice walk matches diffusion only after rescaling, and many walkers are just independent repetitions of one walk.
• Finite runs. Recurrence is a statement about infinite time. A simulation can show the return probability rising toward 1, never reaching it.
• Theory curves. The binomial and the renewal-equation return probability are exact; the 2-D Rayleigh curve is a large-n approximation.

Related simulations

• OS-013Brownian motionParticles jiggle with Gaussian steps; their mean squared displacement grows as 4Dt, as Einstein predicted.
• OS-037DiffusionThe heat equation solved two ways, by finite differences and by random walkers, against the exact Gaussian.
• OS-016Gambler's ruinBet 1 unit until broke or at N: exact odds, expected game length, and why a small edge matters.
• OS-010Galton boardBalls bouncing left or right through rows of pegs pile up into a binomial bell curve.

Background: Pseudo-Random Numbers and Seeds in Simulation

References

1. Pólya, G. (1921). Über eine Aufgabe der Wahrscheinlichkeitsrechnung betreffend die Irrfahrt im Straßennetz. Mathematische Annalen, 84, 149–160. doi.org/10.1007/BF01458701
2. Feller, W. (1968). An Introduction to Probability Theory and Its Applications, Vol. 1, 3rd ed., chapters III and XIV. Wiley.
3. Grinstead, C. M. & Snell, J. L. Introduction to Probability, chapter 12 (random walks). American Mathematical Society. math.dartmouth.edu/~prob/prob/prob.pdf