12 - Markov Chains and Stochastic Processes

It is raining today. What will the weather be in ten days, and how many rainy days do you expect in a year? If you assume that tomorrow's weather depends only on today's, you have already built a Markov chain, and both questions have an answer that comes from a matrix multiplication and a linear system. This article, the twelfth in the Probability and Statistics series, starts from a three-state weather model and arrives at three different uses of the same idea: long-run behaviour, the gambling game that ends in ruin, and PageRank. It is also the foundation of the Metropolis sampler seen in the article on Bayesian estimation. Every result is computed exactly and confirmed with a fixed-seed simulation.

What a Markov chain is and how to write it

A stochastic process is a sequence of random variables X₀, X₁, X₂, … describing a phenomenon over time. It is called a Markov chain (in discrete time and with a finite number of states) if the distribution of the next step, given the whole past, depends only on the current state: P(X_t+1 = j | X_t = i, past) = P(X_t+1 = j | X_t = i). If this probability does not change with t, the chain is homogeneous and is summarised by a transition matrix P, where the entry in row i and column j is P_ij, the probability of moving from i to j. Every row is a probability distribution: it contains non-negative numbers that sum to 1.

In the weather model the three states are sunny, cloudy, rainy. After a sunny day the sun returns with probability 0.7, cloud with 0.2, rain with 0.1; the other two rows are chosen plausibly but arbitrarily. The listing shows three facts. The first is the Chapman-Kolmogorov equation: n-step probabilities are powers of the matrix, because to go from i to j in m + n steps you pass through an intermediate state, and the sum over intermediate states is exactly the matrix product. The second is that after about ten steps the rows of Pⁿ are almost equal to one another: the starting point has been forgotten. The third is the value they converge to, the stationary distribution π, which satisfies π P = π and sums to 1. The program obtains it by solving the system in exact fractions, so that the result does not depend on rounding, and then compares it with the visit frequency of a long simulated trajectory.

from fractions import Fraction as F
import numpy as np

states = ["sunny", "cloudy", "rainy"]
P = np.array([[0.7, 0.2, 0.1],
              [0.3, 0.4, 0.3],
              [0.2, 0.3, 0.5]])
assert np.allclose(P.sum(axis=1), 1)          # every row is a distribution

# n-step probabilities: Chapman-Kolmogorov says P^(m+n) = P^m P^n
P2, P10 = np.linalg.matrix_power(P, 2), np.linalg.matrix_power(P, 10)
print("P^2 row 'sunny'     :", np.round(P2[0], 4))
print("P^10 row 'sunny'    :", np.round(P10[0], 4))
print("P^10 row 'rainy'    :", np.round(P10[2], 4))

# Stationary distribution pi (pi P = pi, sum 1), in exact fractions:
# I drop the last equation (redundant) and replace it with the normalisation
M = [[F(str(P[j][i])) - (1 if i == j else 0) for j in range(3)] for i in range(3)]  # (P^T - I)
M[2] = [F(1), F(1), F(1)]
b = [F(0), F(0), F(1)]
for c in range(3):                              # Gauss-Jordan elimination
    r = next(r for r in range(c, 3) if M[r][c] != 0)
    M[c], M[r], b[c], b[r] = M[r], M[c], b[r], b[c]
    piv = M[c][c]
    M[c] = [x / piv for x in M[c]]; b[c] /= piv
    for r in range(3):
        if r != c:
            f = M[r][c]
            M[r] = [x - f * y for x, y in zip(M[r], M[c])]; b[r] -= f * b[c]
print("exact pi            :", [str(x) for x in b], "=", [round(float(x), 5) for x in b])
pi = np.array([float(x) for x in b])
print("check pi P = pi     :", np.allclose(pi @ P, pi))

# Eigenvalues: 1 and two values of modulus below 1
print("eigenvalues         :", np.round(np.sort(np.linalg.eigvals(P).real)[::-1], 4))

# Fixed-seed simulation: fraction of days in each state along a single trajectory
rng = np.random.default_rng(12)
N = 300_000
x = 0
visits = np.zeros(3)
cum = P.cumsum(axis=1)
u = rng.random(N)
for t in range(N):
    x = int(np.searchsorted(cum[x], u[t]))
    visits[x] += 1
print("simulated frequencies:", np.round(visits / N, 5))
print("maximum gap          :", f"{np.abs(visits / N - pi).max():.5f}")

# Output:
# P^2 row 'sunny'     : [0.57 0.25 0.18]
# P^10 row 'sunny'    : [0.4568 0.2825 0.2607]
# P^10 row 'rainy'    : [0.4562 0.2827 0.2611]
# exact pi            : ['21/46', '13/46', '6/23'] = [0.45652, 0.28261, 0.26087]
# check pi P = pi     : True
# eigenvalues         : [1.     0.4732 0.1268]
# simulated frequencies: [0.45557 0.28303 0.2614 ]
# maximum gap          : 0.00095

From a sunny day, after two days the probability of sun falls from 0.7 to 0.57, and after ten days to 0.4568: the rainy row gives 0.4562, almost the same. The exact solution is π = (21/46, 13/46, 6/23) = (0.4565, 0.2826, 0.2609): in the long run 45.7% of days are sunny and 26.1% rainy, however the series began. On 300,000 simulated days the maximum gap from the exact frequencies is 0.00095, consistent with the error expected for 300,000 correlated steps, which I estimate at around 0.0015 (successive days are correlated, so the error is larger than that of an independent sample). The eigenvalues of P are 1, 0.4732 and 0.1268: the eigenvalue 1 corresponds to π, and the second in order of size governs how fast the starting point is forgotten, as you will see shortly.

When the chain forgets where it started

Not all chains behave this way, and there are two conditions. A chain is irreducible if from every state you can reach every other in a finite number of steps: otherwise the set of states splits into parts that do not communicate. It is aperiodic if returning to a state is not tied to multiples of an integer greater than 1. The convergence theorem for finite chains (Levin, Peres and Wilmer, Theorem 4.9) states that an irreducible aperiodic chain has a unique stationary distribution π and that Pⁿ(i, ·) tends to π for every starting state i. If the chain is only irreducible, π exists and is unique, but the convergence of Pⁿ may fail.

The last part of the next listing is the smallest counterexample: two states A and B that always swap (A goes to B with probability 1 and vice versa). The chain is irreducible but has period 2, and the powers of P do not converge: A^10 is the identity, A^11 the swap. What does converge is the time average: the mean of the first 1000 powers is 0.5 in every entry, that is, the stationary (1/2, 1/2). This is why in the weather model, and in general when you use a chain to simulate, it is the fraction of time spent in each state that always converges, while the distribution at step n converges only in the absence of periodicity.

How long it takes to forget, and how long it takes to return

The program below measures two quantities. The first is the mean return time: starting from a state, how many steps on average it takes to pass through it again. For an irreducible chain on a finite space it equals 1/π_i (this is Kac's lemma, in the text by Levin, Peres and Wilmer): if 45.7% of days are sunny, between one sunny day and the next an average of 1/0.4565 = 2.19 days go by. The second is the total variation distance from the stationary distribution, TV(n), that is, half the sum of the absolute values of the differences between Pⁿ(i, ·) and π, taken at the worst starting state: it measures how far a distribution at step n can still be told apart from π.

import numpy as np

P = np.array([[0.7, 0.2, 0.1],
              [0.3, 0.4, 0.3],
              [0.2, 0.3, 0.5]])
pi = np.array([21/46, 13/46, 6/23])             # stationary, computed in the previous listing

# 1) Mean return time: for an irreducible chain it is 1/pi_i (Kac's formula)
rng = np.random.default_rng(5)
cum = P.cumsum(axis=1)
def return_time(state):
    x, t = state, 0
    while True:
        x = int(np.searchsorted(cum[x], rng.random())); t += 1
        if x == state:
            return t
for i, name in enumerate(["sunny", "cloudy", "rainy"]):
    r = np.array([return_time(i) for _ in range(40_000)])
    print(f"return to {name:6s}: simulated {r.mean():.3f} +- {r.std() / np.sqrt(len(r)):.3f}   1/pi = {1/pi[i]:.3f}")

# 2) How fast the starting point is forgotten: total variation distance
#    TV(n) = max over rows of (1/2) * sum |P^n[i, j] - pi[j]|
lam2 = np.sort(np.abs(np.linalg.eigvals(P)))[-2]
print(f"second eigenvalue in modulus: {lam2:.4f}")
Pn = np.eye(3)
for n in range(1, 9):
    Pn = Pn @ P
    tv = 0.5 * np.abs(Pn - pi).sum(axis=1).max()
    print(f"n={n}  TV={tv:.6f}  TV/lam2^n={tv / lam2**n:.4f}")

# 3) Detailed balance: pi_i P_ij == pi_j P_ji ?
flow = pi[:, None] * P
print("weather chain reversible:", np.allclose(flow, flow.T))
Q = np.array([[0.5, 0.5, 0.0, 0.0],             # walk on 4 nodes in a row, with sticky ends
              [0.25, 0.5, 0.25, 0.0],
              [0.0, 0.25, 0.5, 0.25],
              [0.0, 0.0, 0.5, 0.5]])
w, v = np.linalg.eig(Q.T)
q = np.real(v[:, np.argmax(np.real(w))]); q /= q.sum()
print("stationary of Q:", np.round(q, 4), " reversible:", np.allclose(q[:, None] * Q, (q[:, None] * Q).T))

# 4) Periodicity: two states that alternate
A = np.array([[0.0, 1.0], [1.0, 0.0]])
print("A^10 =", np.linalg.matrix_power(A, 10).tolist(), " A^11 =", np.linalg.matrix_power(A, 11).tolist())
print("mean of the first 1000 powers:", np.round(sum(np.linalg.matrix_power(A, k) for k in range(1, 1001)) / 1000, 3).tolist())

# Output:
# return to sunny : simulated 2.194 +- 0.013   1/pi = 2.190
# return to cloudy: simulated 3.537 +- 0.018   1/pi = 3.538
# return to rainy : simulated 3.869 +- 0.024   1/pi = 3.833
# second eigenvalue in modulus: 0.4732
# n=1  TV=0.256522  TV/lam2^n=0.5421
# n=2  TV=0.126522  TV/lam2^n=0.5650
# n=3  TV=0.060522  TV/lam2^n=0.5712
# n=4  TV=0.028722  TV/lam2^n=0.5728
# n=5  TV=0.013602  TV/lam2^n=0.5733
# n=6  TV=0.006438  TV/lam2^n=0.5734
# n=7  TV=0.003047  TV/lam2^n=0.5734
# n=8  TV=0.001442  TV/lam2^n=0.5734
# weather chain reversible: False
# stationary of Q: [0.1667 0.3333 0.3333 0.1667]  reversible: True
# A^10 = [[1.0, 0.0], [0.0, 1.0]]  A^11 = [[0.0, 1.0], [1.0, 0.0]]
# mean of the first 1000 powers: [[0.5, 0.5], [0.5, 0.5]]

The simulated return times (2.194 ± 0.013; 3.537 ± 0.018; 3.869 ± 0.024, with 40,000 returns each) agree with 1/π within two standard errors: 2.190, 3.538, 3.833. The TV distance falls from 0.2565 at the first step to 0.0014 at the eighth, and the ratio TV(n)/|λ₂|ⁿ, with λ₂ = 0.4732 the second eigenvalue in modulus, settles at 0.5734: the decay is geometric and its ratio is the second eigenvalue. How many steps does it take for TV to fall below 0.01 (a usual definition of mixing time)? Six are enough, and the number grows logarithmically if you ask for greater precision. For chains with many states, a second eigenvalue close to 1 signals slow convergence: this is why MCMC samplers with too small a step stay for a long time near their starting point.

Detailed balance: the chain you can reverse in time

A chain with stationary distribution π is reversible if detailed balance holds: π_i P_ij = π_j P_ji for every pair of states, that is, the average flow from i to j equals that from j to i. The condition implies π P = π and is the property that is built on purpose when designing a sampler: in the Metropolis algorithm the acceptance probability is chosen precisely so that detailed balance holds with respect to the target distribution (Levin, Peres and Wilmer, chapter 3). The third part of the previous listing checks it on two chains. The weather model is not reversible: π_sunny P_{sunny→cloudy} = 0.0913 against π_cloudy P_{cloudy→sunny} = 0.0848. The walk on four nodes in a row, by contrast, is, with stationary (1/6, 1/3, 1/3, 1/6): every birth-death chain, in which you can only move to neighbours, is.

This explains an asymmetry that surprises beginners: a non-reversible chain can converge perfectly well, like the weather model. Reversibility is not necessary for a stationary distribution to exist, it is only a sufficient and convenient condition to check, because checking it requires looking at pairs of states one at a time, instead of solving a system of equations.

Chains with absorbing states: the gambler's ruin

A state is absorbing if once reached it is never left. In a chain with absorbing states it is convenient to reorder the matrix in blocks: Q contains the transitions among transient states and R those towards absorbing states. The fundamental matrix N = (I − Q)⁻¹ has as its (i, j) entry the expected number of visits to transient state j starting from i. The sum of a row of N is the expected number of steps before absorption, and N R gives the probability of ending in each absorbing state. The classic case is the gambler's ruin: you start with 3 euros, in each round you win one euro with probability p or lose one with probability 1 − p, and you stop at 0 or at 10.

import numpy as np

# Gambler's ruin: capital 0..N, each round +1 with prob. p, -1 with prob. 1-p; 0 and N absorb
def ruin(N, p):
    n_trans = N - 1                              # transient states: 1..N-1
    Q = np.zeros((n_trans, n_trans)); R = np.zeros((n_trans, 2))   # R: [0, N]
    for i in range(1, N):
        if i + 1 < N: Q[i-1, i] = p
        else:         R[i-1, 1] = p
        if i - 1 > 0: Q[i-1, i-2] = 1 - p
        else:         R[i-1, 0] = 1 - p
    Nf = np.linalg.inv(np.eye(n_trans) - Q)      # fundamental matrix: expected visits to each state
    return Nf @ R, Nf.sum(axis=1)                # absorption probabilities, expected duration

N = 10
for p in (0.5, 0.45):
    absorb, duration = ruin(N, p)
    print(f"p={p}: from 3 euros P(reaching {N}) = {absorb[2, 1]:.5f}, expected duration = {duration[2]:.4f} rounds")
print("p=0.5, closed formulas: i/N =", 3/N, " and i(N-i) =", 3 * (N - 3))

rng = np.random.default_rng(31)
def game(i, N, p):
    t = 0
    while 0 < i < N:
        i += 1 if rng.random() < p else -1; t += 1
    return i == N, t
R_ = np.array([game(3, N, 0.45) for _ in range(100_000)])
v = R_[:, 0].mean(); d = R_[:, 1].astype(float)
print(f"p=0.45 simulated: P(reaching) = {v:.5f} +- {np.sqrt(v*(1-v)/len(R_)):.5f}, duration = {d.mean():.3f} +- {d.std()/np.sqrt(len(d)):.3f}")

# Output:
# p=0.5: from 3 euros P(reaching 10) = 0.30000, expected duration = 21.0000 rounds
# p=0.45: from 3 euros P(reaching 10) = 0.12825, expected duration = 17.1748 rounds
# p=0.5, closed formulas: i/N = 0.3  and i(N-i) = 21
# p=0.45 simulated: P(reaching) = 0.12860 +- 0.00106, duration = 17.261 +- 0.055

With a fair coin (p = 0.5) the matrix calculation returns a probability of 0.30000 of reaching 10 and an expected duration of 21 rounds, exactly the closed formulas i/N and i(N − i) (Levin, Peres and Wilmer, §2.1). With p = 0.45, a house edge of only five percentage points below 50%, the probability of reaching 10 collapses to 0.12825 and the expected duration drops to 17.17 rounds. The simulation with 100,000 games confirms it: 0.12860 ± 0.00106 and 17.261 ± 0.055. The practical message is that a small asymmetry per round becomes a large asymmetry in the final result, and it is the reason why insisting on playing against a house with a slight edge almost always ends in ruin.

PageRank: a chain to rank pages

PageRank (Brin and Page, 1998) assigns each web page the importance that a random surfer would give it: from a page it follows one of its links with probability d (0.85 in the original work), or jumps to any page with probability 1 − d. The importance is the stationary distribution of this chain. The random jump serves a precise mathematical reason: it guarantees that the chain is irreducible and aperiodic even if the web graph has isolated pages or groups of pages that do not link to each other, and therefore that π exists, is unique, and can be computed by iterating the product x ← x M (the power method) from any starting point. The listing does this on five pages and checks the result against the left eigenvector.

import numpy as np

# PageRank on 5 pages: link[i] = pages that i points to
link = {0: [1, 2], 1: [2], 2: [0], 3: [2, 4], 4: [3]}
n, d_ = 5, 0.85                                   # d_: probability of following a link
G = np.zeros((n, n))
for i, out in link.items():
    for j in out: G[i, j] = 1 / len(out)
M = d_ * G + (1 - d_) / n                         # with prob. 1-d_ we jump at random: the chain is irreducible and aperiodic
x = np.full(n, 1 / n)
for it in range(1, 200):
    new = x @ M
    if np.abs(new - x).sum() < 1e-12: break
    x = new
w, v = np.linalg.eig(M.T)
e = np.real(v[:, np.argmax(np.real(w))]); e /= e.sum()
print(f"power iteration: {it} steps", np.round(x, 5))
print("left eigenvector:", np.round(e, 5), " equal:", np.allclose(x, e))

# Output:
# power iteration: 55 steps [0.32731 0.16911 0.34977 0.08689 0.06693]
# left eigenvector: [0.32731 0.16911 0.34977 0.08689 0.06693]  equal: True

The iteration converges in 55 steps, with tolerance 10^-12 on the sum of the differences, and coincides with the eigenvector: page 2, linked by 0, 1 and 3, has the greatest importance (0.3498), followed by 0 (0.3273). Pages 3 and 4, which link only to each other and receive no link from the other three, stay at the bottom (0.0869 and 0.0669). The number of steps depends on d: the closer d is to 1 the slower the convergence, because the weight of the random jump that reshuffles decreases. This is a model of PageRank, not the ranking system of a search engine, which uses many other signals and is not described in this article.

A worked exercise and a test to do on your own

A machine has two states, working (F) and broken (G). Every day a working machine breaks with probability 0.1 and a broken one is repaired with probability 0.4. The stationary distribution satisfies π_F · 0.1 = π_G · 0.4 (detailed balance, valid for every two-state chain): π_G = π_F / 4, hence π_F = 0.8 and π_G = 0.2. The machine is broken one day in five. The mean return time to the broken state is 1/0.2 = 5 days, while each breakdown lasts on average 1/0.4 = 2.5 days.

What you can do now

You can write the transition matrix of a phenomenon, compute n-step probabilities with powers, obtain the stationary distribution in exact form, and recognise when irreducibility and aperiodicity guarantee convergence. You can use the return time 1/π and the second eigenvalue to estimate the mixing speed, check detailed balance, and solve a chain with absorbing states with the fundamental matrix. These tools appear in the MCMC samplers of the previous article, and in the Markov decision process of reinforcement learning. For the machine learning side see also Markov chains: the property that forgets the past. The next article in the series moves on to non-parametric tests.

The same series also includes Exploratory Data Analysis (EDA) with Python.