Markov Chain Simulation

What the simulation shows

A Markov chain is described by a transition matrix P: entry P_ij is the probability of moving from state i to state j in one step, and every row sums to 1. The lab offers four small chains, each chosen for one lesson:

• Weather — three states, every transition possible. It forgets its starting state within a few steps.
• Periodic 4-cycle — from each state, step to a neighbor with probability ½. The chain alternates between even and odd states forever, so the probability of being in a given state at step t never settles, even though the time spent in each state does.
• Two clusters — A and B talk to each other, C and D talk to each other, and the link between the pairs is weak. The chain stays in one cluster for a hundred steps or more at a time: slow mixing.
• Absorbing — three levels and a finish that is never left. Every path ends there; the question is when.

Simulating a chain needs one uniform random number per step:

1. Start in state X₀ (here, the first state).
2. Draw u ~ U[0, 1). Walk along row X_t of P, adding up probabilities; the next state is the first j whose running sum exceeds u.
3. Count the visit, and repeat.

Alongside the random path, the lab propagates the exact distribution of X_t by repeated matrix multiplication, μ_t+1 = μ_tP. Press Step to follow single moves on the diagram.

The stationary distribution

A distribution π over the states is stationary if one more step leaves it unchanged:

These are linear equations. One of the equations π(P − I) = 0 is redundant, so the lab replaces it with the normalization and solves the system by Gaussian elimination. For the weather chain the exact answer is π = (21/46, 13/46, 12/46) ≈ (0.457, 0.283, 0.261): in the long run 45.7% of days are sunny, whatever the weather today.

Two different convergence statements hide in the pictures. The ergodic theorem says that for an irreducible chain (every state reachable from every other) the fraction of time spent in each state converges to π with probability 1. That holds for the periodic cycle too, and the blue bars and blue line show it; like any Monte Carlo average, the error shrinks roughly like 1/√t, a straight line of slope −½ on the log–log plot. The convergence theorem says the law of X_t itself converges to π, but only if the chain is also aperiodic. The amber curve shows this second kind of convergence, which is geometric and much faster when it happens at all.

Mixing time and the stickiness slider

How far the chain is from equilibrium is measured by the total-variation (TV) distance, the largest possible difference in the probability of any event:

The lab reports the mixing time as the first t with d(t) ≤ ¼, from the chain's starting state.

The stickiness s replaces P with a “lazy” chain that stays put with extra probability s:

Because πP′ = (1 − s)πP + sπ = π, the stationary distribution does not move: the targets in the bar chart stay fixed while everything else slows down. Every eigenvalue λ of P becomes (1 − s)λ + s, closer to 1, so the chain forgets its start more slowly and runs of the same state get longer, as the trace shows. On the two-cluster chain the mixing time grows from 69 steps at s = 0 to 687 at s = 0.9, and the time spent after 2,000 steps is visibly wrong.

Stickiness also cures periodicity. The 4-cycle has an eigenvalue −1, which is what makes it oscillate; laziness moves it to 2s − 1, inside the unit circle. Any s > 0 makes the law converge, which is why theoretical work often assumes a lazy chain with s = ½. The effect is not monotone here: from state 1, s = ½ gets within ¼ of π in a single step, while s = 0.1 needs 4 and s = 0.9 needs 9.

Absorbing chains and the fundamental matrix

A state that cannot be left is absorbing. If it can be reached from everywhere, π puts all its mass there, and the interesting question is how long the chain takes to arrive. Write Q for the transitions among the transient states. The fundamental matrix N = (I − Q)⁻¹ counts expected visits, and its row sums are the expected absorption times:

For the absorbing preset the answer from L1 is 10.56 steps. Stickiness multiplies I − Q by exactly (1 − s), so the expected time becomes 10.56/(1 − s). A single path can be absorbed much earlier or later than expected, so the lab also runs 1,000 independent copies and reports their mean.

Example result

With the weather chain, no stickiness and 2,000 steps (seed 42), the stationary distribution over sunny, cloudy and rainy days is .457 .283 .261, and the path spent fractions .466 .285 .249 of its time in the three states, a total-variation distance of 0.0119. The exact law of the chain was within ¼ of π after 1 step, and by step 2,000 the difference is 3.1e-16, rounding error.

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.

// Stationary distribution: πP = π with Σπ = 1. Transpose to (Pᵀ − I)π = 0 and
// replace the last (redundant) equation by the normalization Σπᵢ = 1.
export function stationary(P) {
  const n = P.length;
  const A = [], b = [];
  for (let i = 0; i < n; i++) {
    A.push(P.map((row, j) => (i === n - 1 ? 1 : row[i] - (i === j ? 1 : 0))));
    b.push(i === n - 1 ? 1 : 0);
  }
  return solve(A, b);
}

// One step of the chain: pick the next state j with probability P[i][j]
// by walking the cumulative sums of row i.
export function nextState(row, u) {
  let acc = 0;
  for (let j = 0; j < row.length; j++) {
    acc += row[j];
    if (u < acc) return j;
  }
  return row.length - 1; // guards against rounding when the row sums to 1 − ε
}

// Exact law of the chain one step later: μ′ⱼ = Σᵢ μᵢ Pᵢⱼ.
export function evolve(mu, P) {
  return P[0].map((_, j) => { let s = 0; for (let i = 0; i < mu.length; i++) s += mu[i] * P[i][j]; return s; });
}

// Total-variation distance: half the L1 distance between two distributions.
export function tv(a, b) {
  let s = 0;
  for (let i = 0; i < a.length; i++) s += Math.abs(a[i] - b[i]);
  return s / 2;
}

def stationary(P):
    """Solve pi P = pi with sum(pi) = 1: (P^T - I) pi = 0, last row replaced by ones."""
    n = len(P)
    A = [[1 if i == n - 1 else P[j][i] - (1 if i == j else 0) for j in range(n)] for i in range(n)]
    b = [1 if i == n - 1 else 0 for i in range(n)]
    return solve(A, b)


def next_state(row, u):
    """Pick the next state by walking the cumulative sums of the current row."""
    acc = 0
    for j, pj in enumerate(row):
        acc += pj
        if u < acc:
            return j
    return len(row) - 1


def evolve(mu, P):
    """Exact law one step later: mu'_j = sum_i mu_i P_ij."""
    out = []
    for j in range(len(P[0])):
        s = 0
        for i in range(len(mu)):
            s += mu[i] * P[i][j]
        out.append(s)
    return out


def tv(a, b):
    """Total-variation distance."""
    s = 0
    for x, y in zip(a, b):
        s += abs(x - y)
    return s / 2

Run the full model from a terminal:

node simulation.js --chain weather --stickiness 0 --seed 42
python simulation.py --chain weather --stickiness 0 --seed 42

Download the model

Assumptions & limitations

• The Markov property. The next state depends only on the current one. Weather, customers and molecules usually remember more; a common fix is to enlarge the state (for example, today's and yesterday's weather).
• Time-homogeneous. The same matrix applies at every step: no seasons, no learning.
• Finite, small state spaces. Exact solutions by elimination and matrix products are cheap here; chains with millions of states need iterative or sampling methods, as in PageRank and Metropolis–Hastings.
• Mixing time from one start. The reported mixing time starts from the first state, not the worst case used in the textbook definition.
• Floating-point sampling. Rows of P′ can sum to 1 − 10⁻¹⁶; the sampler then returns the last state, a bias far below anything measurable.

Related simulations

• OS-012Random walkCoin-flip steps on a line or grid: distance grows like √n, and in 1-D and 2-D every walker comes back.
• OS-005Metropolis–HastingsA Markov chain that samples any density known up to a constant; tune the step size and watch it mix.
• OS-033PageRankRank pages by where a random surfer spends its time: power iteration vs a simulated walk.
• OS-016Gambler's ruinBet 1 unit until broke or at N: exact odds, expected game length, and why a small edge matters.

Background: How Accurate Is a Monte Carlo Estimate?

References

1. Norris, J. R. (1997). Markov Chains. Cambridge University Press.
2. Levin, D. A., Peres, Y. & Wilmer, E. L. (2017). Markov Chains and Mixing Times, 2nd ed. American Mathematical Society.
3. Grinstead, C. M. & Snell, J. L. Introduction to Probability, chapter 11 (Markov chains). American Mathematical Society. math.dartmouth.edu/~prob/prob/prob.pdf