ML coding · Probabilistic models
Stationary distribution of a Markov chain
Given a row-stochastic transition matrix P (n × n, each row sums to 1) of an
irreducible, aperiodic Markov chain, return its stationary distribution π: the unique probability vector
with πP = π.
Examples
Input P = [[0.9, 0.1], [0.5, 0.5]]
Output [0.8333, 0.1667]
A calm state that is sticky and a volatile one that is not: 5/6 of the time calm.
Input P = [[0, 1], [0.5, 0.5]]
Output [0.3333, 0.6667]
π = (1/3, 2/3). Check: (1/3)(0, 1) + (2/3)(0.5, 0.5) = (1/3, 2/3).
Input P = [[0.5, 0.25, 0.25], [0.25, 0.5, 0.25], [0.25, 0.25, 0.5]]
Output [0.3333, 0.3333, 0.3333]
A symmetric chain: uniform.
Constraints
- 1 ≤ n ≤ 500
- Rows sum to 1; the chain is irreducible and aperiodic
- Checked to a relative tolerance of 10⁻⁹
Hints
Hint 1
πP = π is the same as (Pᵀ − I)πᵀ = 0. Why can you not just solve that system?
Hint 2
The n equations are linearly dependent, because the rows of P sum to 1. One of them carries no information.
Hint 3
Replace the last equation with Σπ = 1 (a row of ones, right-hand side 1) and call np.linalg.solve.
Where candidates go wrong
- Solving (Pᵀ − I)π = 0 directly: the matrix is singular, and the only solution a solver finds is π = 0.
- Using P instead of Pᵀ: that finds the right eigenvector, which is the all-ones vector.
- Taking an eigenvector of Pᵀ for eigenvalue 1 and forgetting to normalise it, or its sign.
The idea
πP = π transposes to (Pᵀ − I)πᵀ = 0. Because every row of P sums to one, the columns of Pᵀ − I sum to zero, so its rows are dependent and the system has a one-dimensional solution space. Replacing any one equation with the normalisation Σπi = 1 leaves a non-singular system with exactly the stationary distribution as its solution, which one call to a linear solver finds in O(n³).
For very large sparse chains, power iteration π ← πP converges geometrically at a rate set by the second-largest eigenvalue modulus.
Time O(n³) · Space O(n²)
Reference solution
import numpy as np
def stationary(P):
P = np.asarray(P, dtype=float)
n = len(P)
A = P.T - np.eye(n) # pi P = pi, transposed: (P^T - I) pi = 0
A[-1] = 1.0 # the equations are dependent: replace one with sum(pi) = 1
b = np.zeros(n)
b[-1] = 1.0
return np.linalg.solve(A, b)Follow-ups an interviewer asks
- Use power iteration instead and stop when successive vectors differ by less than 10⁻¹².
- Expected return time to state i is 1/π_i. Verify it by simulation.
- PageRank: add a damping factor and teleportation, then find the stationary distribution.
- What goes wrong for a periodic chain such as [[0, 1], [1, 0]]?
Python runs in your browser with NumPy. The first run downloads the runtime, which takes a few seconds.
Why interviewers ask it
Quant interviews ask for this in many disguises: the long-run fraction of days in a high-volatility regime, the share of time a queue is full, the probability a random walk on a graph is at a node. The coding version tests whether you see that the balance equations are dependent and know what to do about it.
The theory behind it is in the lesson Naive Bayes, EM, Viterbi and Markov chains from scratch.
Keep going
Gaussian naive Bayes classifier
Fit class priors, means and variances, then predict by the largest log posterior, vectorised over test points and classes.
Open the problem → Medium20 minOne EM step for a Gaussian mixture
Compute responsibilities in log space, then refit weights, means and variances: one expectation-maximisation step for a 1-D Gaussian mixture, from scratch..
Open the problem → Medium25 minThe Viterbi algorithm for a hidden Markov model
Find the most likely sequence of hidden states in an HMM with dynamic programming in log space, with back-pointers to recover the path.
Open the problem → Easy12 minStationary distribution of a Markov chain
Solve πP = π with Σπ = 1 by replacing one redundant equation, in NumPy.
Open the problem →