Chains
Discrete-time Markov chains answered by linear algebra — stationary distributions, reversibility, absorption, branching, urns and fitting a chain to a path.
This module covers Markov chains in discrete time with a finite set of states, together with Galton–Watson branching, the Ehrenfest model and Pólya’s urn. Most of the questions a course answers by simulation are answered here in closed form. The stationary distribution is a linear solve, the expected time to absorption and the absorption probabilities come from one matrix inverse, and the extinction probability of a branching process is a fixed-point iteration. A simulated chain is still available, one step at a time, and its long-run behaviour can be held against the solved answer as an independent check.
Conventions. A transition matrix p is n×n and row-stochastic: \(p_{ij}\) is the probability of moving from state \(i\) to state \(j\), and every row sums to one. A distribution over states is a vector n that sums to one; a state is written as a one-hot vector. Randomness is never drawn inside the module: a simulated step takes a uniform draw u in \([0, 1)\) from the caller. In a process declaration (see the Chain, Absorbing, Urn and Transition sample models) the draws come from the declared stream.
Limitations. stationary assumes the chain has a unique stationary distribution (a single closed communicating class). The mixing time uses the eigenvalues from Linear algebra, which are computed by an unshifted QR iteration designed for symmetric matrices: it is reliable for reversible chains and for chains with real, well-separated eigenvalues, and not for chains with complex eigenvalues. There is no n-step transition matrix \(P^n\) in this module; use afterOne repeatedly, or stationary for the limit.
open Prelude.Chains
Where a chain is after one step
A distribution \(\mu\) over the states moves to \(\mu P\) after one step.
afterOne
afterOne(start, p) · start distribution vector n · p transition matrix n×n → distribution vector n
Let \(\mu \in \mathbb{R}^n\) (the argument start) be a distribution over \(n\) states and let \(P \in \mathbb{R}^{n \times n}\) (the argument p) be a row-stochastic transition matrix. Then
The distribution after one step from start. Passing a one-hot start returns the row of p for that state.
Intuition Picture probability as sand spread over the states. One step pours each state’s sand out along its row of \(P\), and the new amount at state \(t\) is everything that flowed in. It moves the whole distribution at once, rather than following one random path.
Use it when you want a short-horizon forecast of how a population splits across states: next month’s share of active, dormant and churned customers, or tomorrow’s weather odds. Call it repeatedly for \(n\) steps, since there is no \(P^n\) here. For where the chain ends up in the long run use stationary, and to simulate a single trajectory use landed or walk.
state-transition-modellingforecasting
def afterOne(start, p) ≔ einsum("s,st->t", start, p)
Building blocks
composed(a, b)— Let \(A \in \mathbb{R}^{n \times m}\) and \(B \in \mathbb{R}^{m \times k}\). Then \(\operatorname{composed}(A, B)_{ik} = \sum_{j} A_{ij} B_{jk} = (AB)_{ik}\). The matrix product of two matrices, written as one contraction.
The stationary distribution, by a solve
A stationary distribution satisfies \(\pi P = \pi\) with \(\sum_i \pi_i = 1\). The first condition alone is singular, because every row of \(P\) sums to one, so the module adds a matrix of ones \(J\), which enforces the normalisation and makes the system non-singular:
\[(P^\top - I + J)\,\pi = \mathbf{1}\]The solver in Linear algebra is conjugate gradient, which is valid only for symmetric positive-definite matrices, and \(P^\top - I + J\) is not symmetric unless the chain is reversible. The system is therefore solved through its normal equations \(A^\top A\,\pi = A^\top \mathbf{1}\), whose matrix is symmetric positive-definite and whose answer is the same.
stationary
stationary(p) · p transition matrix n×n → stationary distribution vector n
Let \(P \in \mathbb{R}^{n \times n}\) (the argument p) be a row-stochastic matrix with a unique stationary distribution, and let
Then
\[\operatorname{stationary}(P) = \operatorname{solve}\big(A^{\top} A,\; A^{\top} \mathbf{1}_n\big),\]where solve is the preconditioned conjugate gradient of Linear algebra, run for \(n\) turns. In exact arithmetic this is \(\pi = (A^{\top} A)^{-1} A^{\top} \mathbf{1}_n = A^{-1} \mathbf{1}_n\), the vector with \(\pi^{\top} P = \pi^{\top}\) and \(\sum_i \pi_i = 1\).
The stationary distribution \(\pi\) of the chain, found by one linear solve rather than by running the chain. The chain must have a unique stationary distribution.
Intuition \(\pi\) is the spread of sand that one step leaves unchanged: every state receives exactly as much as it sends out. It is also the fraction of time a long run spends in each state, whatever the starting point, so a linear solve replaces a long simulation.
Use it when you need long-run shares: steady-state market share, PageRank-style importance of pages, the availability of a machine that fails and gets repaired. The chain needs a single closed class; for a chain with absorbing states use absorption instead. Conjugate gradient runs a fixed \(n\) turns on the normal equations, so on a badly conditioned chain check that afterOne(stationary(p), p) returns the same vector.
long-run-behaviourrankingstate-transition-modelling
def stationary(p) ≔ solve(normalOf(balance(p)), rhsOf(balance(p), everyOne(p)))
Building blocks
balance(p)— Let \(P \in \mathbb{R}^{n \times n}\). Then \(\operatorname{balance}(P) = P^{\top} - I_n + \mathbf{1}_n \mathbf{1}_n^{\top}\), that is \(\operatorname{balance}(P)_{ij} = P_{ji} - [\,i = j\,] + 1\). The non-singular balance matrix Pᵀ − I + J.everyOne(p)— Let \(P \in \mathbb{R}^{n \times n}\). Then \(\operatorname{everyOne}(P)_i = \frac{1}{n} \sum_{j=1}^{n} 1 = 1\), so \(\operatorname{everyOne}(P) = \mathbf{1}_n\). A vector of ones with one entry per state, the right-hand side of the balance system.normalOf(a)— Let \(A \in \mathbb{R}^{m \times n}\). Then \(\operatorname{normalOf}(A)_{ij} = \sum_{k} A_{ki} A_{kj} = (A^{\top} A)_{ij}\). The normal-equations matrix AᵀA.rhsOf(a, b)— Let \(A \in \mathbb{R}^{m \times n}\) and \(b \in \mathbb{R}^m\). Then \(\operatorname{rhsOf}(A, b)_i = \sum_{k} A_{ki} b_k = (A^{\top} b)_i\). The normal-equations right-hand side Aᵀb.
Reversibility
A chain is reversible when it satisfies detailed balance, \(\pi_i p_{ij} = \pi_j p_{ji}\) for every pair of states. Detailed balance implies stationarity but not the other way round, so it is worth measuring. The time-reversed chain has transitions
\[\tilde{p}_{ij} = \frac{\pi_j\, p_{ji}}{\pi_i}\]and equals \(P\) exactly when the chain is reversible. How quickly a chain forgets its starting state is governed by the second-largest eigenvalue modulus \(\lvert \lambda_2 \rvert\), through the relaxation time \(1 / (1 - \lvert \lambda_2 \rvert)\).
flow
flow(p, dist) · p transition matrix n×n · dist distribution vector n → probability flow n×n
Let \(P \in \mathbb{R}^{n \times n}\) (the argument p) be a transition matrix and \(\pi \in \mathbb{R}^n\) (the argument dist) a distribution over the states. Then
The matrix of probability flows: how much probability moves from \(i\) to \(j\) in one step when the chain is distributed as dist.
Intuition \(\pi_i P_{ij}\) is the traffic on the one-way street from \(i\) to \(j\): the share of all steps that are an \(i \to j\) move when the chain is spread as dist. The row sums give what leaves each state and the column sums what arrives.
Use it when you want to know which transitions carry the volume, such as the busiest hand-offs in a customer journey or the links that bottleneck a workflow, or to compare the traffic in both directions of a pair. Pass stationary(p) for long-run traffic, or the current distribution for the next step’s. It is the input that detailedBalance and reverseOf are built on.
flow-analysisstate-transition-modellinglong-run-behaviour
def flow(p, dist) ≔ einsum("i,ij->ij", dist, p)
detailedBalance
detailedBalance(p, dist) · p transition matrix n×n · dist distribution vector n → largest imbalance scalar
Let \(P \in \mathbb{R}^{n \times n}\) (the argument p) be a transition matrix and \(\pi \in \mathbb{R}^n\) (the argument dist) a distribution. Then
The largest violation of detailed balance. Pass stationary(p) as dist; an answer of zero (to rounding) means the chain is reversible.
Intuition Stationarity only asks that each state’s total inflow equal its total outflow. Detailed balance asks more: every pair of states trades equally in both directions, so there is no net circulation. Such a chain looks the same played backwards, like molecules jiggling at equilibrium rather than a conveyor belt going round.
Use it when you verify an MCMC sampler (Metropolis–Hastings is built to satisfy detailed balance), decide whether to trust mixing (its eigenvalues are reliable for reversible chains), or look for a one-way cycle, such as a lifecycle that always circulates in the same direction. The answer is in probability units, and stationary is an iterative solve, so compare against a small tolerance rather than exact zero.
model-validationsampler-validation
def detailedBalance(p, dist) ≔ maximum(abs(flow(p, dist) − einsum("ij->ji", flow(p, dist))))
reverseOf
reverseOf(p, dist) · p transition matrix n×n · dist distribution vector n → reversed transition matrix n×n
Let \(P \in \mathbb{R}^{n \times n}\) (the argument p) be a transition matrix and \(\pi \in \mathbb{R}^n\) (the argument dist) a distribution, normally the stationary one. Then
The transition matrix of the chain run backwards in time. States with zero probability are guarded against division by zero.
Intuition Run the film backwards: given the chain is at \(i\) now, where was it one step ago? It is Bayes’ rule with \(\pi\) as the prior, the flow into \(i\) from each \(j\) divided by the chance of being at \(i\). It equals \(P\) exactly when the chain is reversible.
Use it when you reason backwards from an observation: the page a visitor most likely came from, or the stage a ticket was in before it reached this one. It also shows which pairs break reversibility when you compare it with p entry by entry, where detailedBalance gives only the worst gap. A row for a state with zero probability divides by \(10^{-300}\) and is meaningless; ignore it.
backward-inferencemodel-validation
def reverseOf(p, dist) ≔ einsum("ji,i->ij", flow(p, dist), 1 / nonzero(dist))
nonzero
nonzero(v) · v any value tensor → guarded value tensor
Let \(v\) be a tensor of any shape. Then, entrywise,
\[\operatorname{nonzero}(v) = v + [\,\lvert v \rvert < 10^{-300}\,] \cdot 10^{-300} .\]Every entry of magnitude below \(10^{-300}\) (in practice, every exact zero) is increased by \(10^{-300}\), so the result can safely be used as a divisor or inside ln. Other entries are unchanged. Jumps and Points use it as their guard too.
Intuition A nudge that only exact zeros (and anything smaller than \(10^{-300}\)) feel. Dividing by the result gives a huge finite number instead of infinity or NaN, and \(0 \cdot \ln(\operatorname{nonzero}(0))\) comes out as a clean 0. That matters here because a mask computes both of its sides, so one NaN would spoil the whole expression.
Use it when you divide by a total, count or probability that can be zero (normalising a row that was never visited, a rate over an empty window) or take the log of a probability that can be zero. It is not smoothing: to regularise counts, add a pseudo-count before normalising. Paths has an identical guard named guarded.
numerical-stabilitynormalisation
def nonzero(v) ≔ v + (abs(v) < 1e-300) · 1e-300
secondOf
secondOf(p) · p transition matrix n×n → second eigenvalue modulus scalar
Let \(P \in \mathbb{R}^{n \times n}\) (the argument p) and let \(\lambda_1, \dots, \lambda_n\) be the values returned by \(\operatorname{eigenvalues}(P)\) from Linear algebra. Then
The largest eigenvalue modulus below one, \(\lvert \lambda_2 \rvert\). Every eigenvalue of modulus \(0.9999999\) or more is excluded, so for a periodic chain this is the largest of the others; if every eigenvalue is excluded the result is 0. See the limitation on eigenvalues above.
Intuition Any starting distribution is the stationary one plus components that shrink by a factor \(\lvert \lambda \rvert\) every step. The largest modulus below one belongs to the slowest-dying memory of the start. A value near 1 means a bottleneck: groups of states the chain rarely moves between.
Use it when you diagnose slow convergence or near-separate clusters, such as communities in a network or regimes that rarely switch, or compare two designs of a chain. Most of the time you want mixing, which turns this into a number of steps. Trust it for reversible chains; with complex eigenvalues the QR iteration is unreliable. A periodic chain never settles, even though this still returns a value below one.
convergence-diagnosticslong-run-behaviour
def secondOf(p) ≔ maximum(abs(spectrum(p)) · (abs(spectrum(p)) < 0.9999999))
mixing
mixing(p) · p transition matrix n×n → relaxation time scalar
Let \(P \in \mathbb{R}^{n \times n}\) (the argument p) and \(\lambda_\ast = \operatorname{secondOf}(P)\). Then
The relaxation time \(1 / (1 - \lvert \lambda_2 \rvert)\): roughly how many steps the chain needs before its starting state stops mattering.
Intuition The start’s influence fades roughly like \(\lvert \lambda_2 \rvert^{t}\), which shrinks by a factor \(e\) about every \(1/(1 - \lvert \lambda_2 \rvert)\) steps. After a few of these relaxation times the chain has forgotten where it began.
Use it when you choose the burn-in for an MCMC run or a simulation (discard several relaxation times), estimate how many cycles a system takes to return to steady state after a shock, or compare how fast two card shuffles randomise. It is an order of magnitude, not a guaranteed bound. It inherits secondOf’s limits: meaningless for a periodic chain, unreliable with complex eigenvalues, and 1 when every eigenvalue was excluded.
convergence-diagnosticssimulationlong-run-behaviour
def mixing(p) ≔ 1 / nonzero(1 − secondOf(p))
Building blocks
spectrum(p)— Let \(P \in \mathbb{R}^{n \times n}\). Then \(\operatorname{spectrum}(P) = \operatorname{eigenvalues}(P) \in \mathbb{R}^n\), the eigenvalues from Linear algebra’s QR iteration. The eigenvalues of the transition matrix; used bysecondOf.
Absorption
Order the states so that the transient ones come first. The transition matrix then has a transient block \(Q\) (transient to transient) and a block \(R\) (transient to absorbing). The fundamental matrix
\[N = (I - Q)^{-1}\]holds in \(N_{ij}\) the expected number of visits to transient state \(j\) starting from transient state \(i\). Its row sums are the expected number of steps before absorption, and \(NR\) gives the probability of ending in each absorbing state.
fundamental
fundamental(q) · q transient block t×t → fundamental matrix t×t
Let \(Q \in \mathbb{R}^{t \times t}\) (the argument q) be the transient-to-transient block of an absorbing chain, so that \(I_t - Q\) is non-singular. Then
where inverse is forty Newton–Schulz turns; to that accuracy \(N = (I_t - Q)^{-1}\).
The fundamental matrix, computed with inverse from Linear algebra.
Intuition \((I - Q)^{-1} = I + Q + Q^2 + \cdots\): entry \((i, j)\) adds up, over every step, the chance of standing on \(j\), which is the expected number of visits to \(j\) before the chain falls into an absorbing state. The series converges because some probability leaks out of the transient states at every step.
Use it when you need the expected time spent in each intermediate stage: days a patient spends in each ward before discharge, months a loan spends in each arrears bucket before it is repaid or written off, visits to each page before a user leaves. For the two usual summaries use expectedSteps and absorption. When absorption is very slow, \(I - Q\) is nearly singular and forty Newton–Schulz turns may fall short, so check that \((I - Q)N \approx I\).
absorption-analysistime-in-state
def fundamental(q) ≔ inverse(identity(q) − q)
expectedSteps
expectedSteps(q) · q transient block t×t → expected steps vector t
Let \(Q \in \mathbb{R}^{t \times t}\) (the argument q) be the transient block and \(N = \operatorname{fundamental}(Q)\). Then
The expected number of steps before absorption from each transient state.
Intuition The time until absorption is the time spent in each transient state, added up. Summing a row of the visit counts \(N\) gives the expected number of steps a chain started there takes before it is caught.
Use it when you want a mean time-to-completion or time-to-failure: how many turns a game lasts, how many attempts a process with retries needs to finish, how many stages of wear a component passes through before it fails. It gives only the mean, not the spread. For a \(\pm 1\) walk between two barriers, the closed form durationOf in Paths needs no matrix at all.
hitting-timesabsorption-analysisreliability-analysis
def expectedSteps(q) ≔ einsum("ij->i", fundamental(q))
absorption
absorption(q, r) · q transient block t×t · r transient-to-absorbing block t×a → absorption probabilities t×a
Let \(Q \in \mathbb{R}^{t \times t}\) (the argument q) be the transient block, \(R \in \mathbb{R}^{t \times a}\) (the argument r) the transient-to-absorbing block, and \(N = \operatorname{fundamental}(Q)\). Then
The probability that the chain started in transient state \(i\) is finally absorbed in absorbing state \(k\). Each row sums to one.
Intuition Every transient state eventually drains into one of the absorbing sinks. \((NR)_{ik}\) counts the expected visits to each transient \(j\) and multiplies by the chance of jumping from \(j\) straight into \(k\), so it is like water from tap \(i\) finding which drain it goes down.
Use it when you need the probability of each final outcome: a lead converts or churns, a loan is repaid or defaults, a game is won or lost, a patient recovers or does not. Every row should sum to one, which is a free check that q and r were split correctly. For the plain \(\pm 1\) walk between two barriers, ruinOf in Paths gives the same answer in closed form.
absorption-analysisoutcome-predictionrisk-measurement
def absorption(q, r) ≔ einsum("ij,jk->ik", fundamental(q), r)
One step of a simulated chain
A categorical draw needs no search. The running totals of a probability vector split \([0, 1)\) into intervals, and the chosen category is the one whose interval contains the uniform draw \(u\):
\[\text{chosen}_k = [\, C_k > u \,]\,[\, C_k - v_k \le u \,], \qquad C_k = \sum_{j \le k} v_j\]runningTotal
runningTotal(v) · v values vector n → cumulative sums vector n
Let \(v \in \mathbb{R}^n\) with components \(v_0, \dots, v_{n-1}\) (positions counted from zero). Then
\[\operatorname{runningTotal}(v)_k = \sum_{j=0}^{n-1} [\,k \ge j\,]\, v_j = \sum_{j=0}^{k} v_j , \qquad k = 0, \dots, n-1 .\]The cumulative sum \(C_k\). It builds an n×n mask, so it suits vectors of up to a few thousand entries. Points uses it to turn gaps into arrival times.
Intuition Lay the values end to end on a ruler: entry \(k\) is where the \(k\)-th piece ends. For probabilities the pieces partition \([0, 1)\) into consecutive intervals, and for waiting gaps the ends are timestamps.
Use it when you turn a probability vector into a cumulative distribution for sampling (as landedOn does), turn gaps into arrival times (as timesOf in Points does), or need cumulative totals such as sales to date. It builds an n×n mask, so memory grows as \(n^2\) and a few thousand entries is the practical limit.
running-totalssimulationevent-arrivals
def runningTotal(v) ≔ einsum("ij,j->i", held(positions(v)) ≥ laid(positions(v)), v)
landedOn
landedOn(v, u) · v probabilities vector n · u uniform draw scalar in [0, 1) → chosen category one-hot vector n
Let \(v \in [0, 1]^n\) be a probability vector with \(\sum_k v_k = 1\), let \(u \in [0, 1)\) be a uniform draw, and let \(C = \operatorname{runningTotal}(v)\). Then
\[\operatorname{landedOn}(v, u)_k = [\, C_k > u \,]\,[\, C_k - v_k \le u \,] , \qquad k = 0, \dots, n-1 .\]Draws one category from the probabilities v using the single uniform u, and returns it as a one-hot vector: the entry is 1 for the one \(k\) with \(u \in [C_k - v_k, C_k)\).
Intuition A roulette wheel with slices sized by v: the uniform draw is where the ball stops, and the chosen slice is the one that contains it. Two comparisons per slice replace any search or branch.
Use it when you draw one option from given probabilities: the next state of a chain, the segment a synthetic customer belongs to, the action a stochastic policy takes, any weighted random choice. The same u always picks the same category, so pass a fresh draw each time. v must sum to one: if it sums to less and u lands past the total, the result is all zeros, meaning no category was chosen. To move a chain state use landed.
categorical-samplingsimulation
def landedOn(v, u) ≔ (runningTotal(v) > u) · (runningTotal(v) − v ≤ u)
landed
landed(at, p, u) · at current state one-hot vector n · p transition matrix n×n · u uniform draw scalar → next state one-hot vector n
Let \(x \in \{0, 1\}^n\) (the argument at) be a one-hot state, \(P \in \mathbb{R}^{n \times n}\) (the argument p) a transition matrix, and \(u \in [0, 1)\) a uniform draw. Then
One step of the chain from the state at: the row of p for that state is used as the probabilities of a categorical draw.
Intuition Moving a token on a board game: read the row of \(P\) for the square the token is on, then spin that row’s roulette wheel to find the next square.
Use it when you step one simulated trajectory whose state is already a valid one-hot vector: a customer’s journey through a funnel, a weather sequence, a board game. Inside a process whose state starts as zeros use walk, which supplies the opening state. To move the whole population in expectation, with no randomness, use afterOne.
simulationstate-transition-modelling
def landed(at, p, u) ≔ landedOn(nextRow(at, p), u)
walk
walk(at, start, p, u) · at current state vector n · start opening state one-hot vector n · p transition matrix n×n · u uniform draw scalar → next state one-hot vector n
Let \(x \in \mathbb{R}^n\) (the argument at) be the current state, either one-hot or all zeros, let \(x^{(0)} \in \{0, 1\}^n\) (the argument start) be the opening one-hot state, \(P \in \mathbb{R}^{n \times n}\) a transition matrix, and \(u \in [0, 1)\) a uniform draw. Then
One step of a simulated chain, suitable as the advance body of a process. A process weight starts as zeros, which is not a state, so while at sums to zero the step starts from start instead.
Intuition landed plus one rule: “if I have not started yet, stand on start first”. A process weight begins as zeros, which is no state at all, and this is the step that notices it.
Use it when you write the advance body of a process that generates chain trajectories, for example to hold the empirical time in each state against stationary, or to feed a simulated path to fittedChain and check that it recovers p. Any all-zero state restarts from start, including one produced when a row of p sums to less than one and the draw lands past its total, so keep rows stochastic. Outside a process, landed does the same job.
simulationstate-transition-modellingmodel-validation
def walk(at, start, p, u) ≔ landed(begun(at, start), p, u)
Building blocks
nextRow(at, p)— Let \(x \in \mathbb{R}^n\) and \(P \in \mathbb{R}^{n \times n}\). Then \(\operatorname{nextRow}(x, P)_t = \sum_s x_s P_{st} = (x^{\top} P)_t\). The row of transition probabilities out of the current state; the same contraction asafterOne.begun(at, start)— Let \(x, x^{(0)} \in \mathbb{R}^n\). Then \(\operatorname{begun}(x, x^{(0)}) = x + [\, \sum_s x_s = 0 \,]\, x^{(0)}\).atitself, orstartwhileatstill sums to zero.
Galton–Watson branching
Each individual has \(k\) children with probability \(f_k\), independently. With generating function \(G(s) = \sum_k f_k s^k\), the extinction probability is the smallest fixed point of \(G\) in \([0, 1]\), and it equals one exactly when the mean offspring \(m = \sum_k k f_k\) is at most one. Iterating \(s \leftarrow G(s)\) from \(s = 0\) converges to that fixed point from below.
offspringMean
offspringMean(f) · f offspring distribution vector K → mean offspring scalar
Let \(f = (f_0, \dots, f_{K-1}) \in \Delta^{K-1}\) be an offspring distribution, where \(f_k\) (positions counted from zero) is the probability of \(k\) children. Then
\[\operatorname{offspringMean}(f) = m = \sum_{k=0}^{K-1} k\, f_k .\]The mean number of children.
Intuition \(m\) is the growth factor per generation: on average each individual is replaced by \(m\) others. Below one the line shrinks, above one it can explode. It plays the role of \(R_0\) in an epidemic.
Use it when you take a first look at any cascade: the reproduction number of an infection, how many sign-ups each new user’s invitations bring, a chain reaction, whether a family name survives. Entry \(k\) of f (counting from zero) is the probability of exactly \(k\) children. It decides subcritical and drives expectedSize; for the probability that the line actually dies out use extinctionOf.
population-dynamicsepidemic-modelling
def offspringMean(f) ≔ einsum("k,k->", f, positions(f))
pgf
pgf(f, s) · f offspring distribution vector K · s argument scalar → generating function value scalar
Let \(f = (f_0, \dots, f_{K-1}) \in \Delta^{K-1}\) be an offspring distribution and \(s \in \mathbb{R}\). Then
\[\operatorname{pgf}(f, s) = G(s) = \sum_{k=0}^{K-1} f_k\, s^k .\]The probability generating function.
Intuition One polynomial that packs the whole offspring distribution. Read \(G(s)\) as “the chance that every child has some property, when each has it independently with probability \(s\)”. With “its line dies out” as the property, the extinction probability \(q\) must satisfy \(q = G(q)\).
Use it when you set up that fixed-point equation, check a distribution (\(G(1)\) should be 1 and \(G(0)\) is the chance of no children), or compose generations: \(G(G(s))\) is the generating function two generations out. Most users call extinctionOf, which iterates this function for them.
extinction-riskpopulation-dynamics
def pgf(f, s) ≔ einsum("k,k->", f, s ^ positions(f))
extinctionOf
extinctionOf(f) · f offspring distribution vector K → extinction probability scalar
Let \(f \in \Delta^{K-1}\) be an offspring distribution with generating function \(G = \operatorname{pgf}(f, \cdot)\). Let \(s_0 = 0\) and \(s_{j+1} = G(s_j)\) for \(j = 0, 1, \dots\). Then
\[\operatorname{extinctionOf}(f) = s_{200} .\]The probability that a line started by one individual eventually dies out, from 200 iterations of \(G\) starting at zero. Convergence is fast away from the critical case; when the mean is very close to one the iteration converges slowly and 200 turns can stop short of the true value.
Intuition A line dies out exactly when every child’s line dies out. Starting from \(s_0 = 0\), each iterate \(s_j\) is the probability of being extinct by generation \(j\), and it climbs toward the true answer. So the result is the chance of dying out within 200 generations.
Use it when you want the chance that a new strain, a rumour or a referral chain fizzles out by itself even though its mean says it could grow. With \(z\) independent founders the answer is \(q^z\). Near the critical mean of one the climb is slow and 200 turns undershoot the true value. For only a yes/no on certain extinction use subcritical.
extinction-riskpopulation-dynamicsrisk-measurement
def extinctionOf(f) ≔ fixedLoop(200, f, 0)
subcritical
subcritical(f) · f offspring distribution vector K → flag 1 or 0
Let \(f \in \Delta^{K-1}\) be an offspring distribution with mean \(m = \operatorname{offspringMean}(f)\). Then
\[\operatorname{subcritical}(f) = [\, m \le 1 \,] .\]1 when the mean offspring is at most one (extinction is certain, unless every individual has exactly one child), 0 otherwise. The critical case \(m = 1\) is included.
Intuition The one threshold that matters: if each individual leaves at most one replacement on average, every line eventually dies, although at exactly \(m = 1\) it can take a very long time.
Use it when you need a quick go/no-go: is an outbreak self-limiting (\(R \le 1\)), can a viral loop sustain itself, is a reactor below criticality. It says nothing about how likely extinction is when \(m > 1\) (use extinctionOf) or how fast it comes. One degenerate case escapes the theorem: if every individual has exactly one child, \(m = 1\) and this returns 1, yet the line never dies.
extinction-riskepidemic-modelling
def subcritical(f) ≔ offspringMean(f) ≤ 1
expectedSize
expectedSize(f, start, n) · f offspring distribution vector K · start initial population scalar · n generations scalar → expected population scalar
Let \(f \in \Delta^{K-1}\) be an offspring distribution with mean \(m = \operatorname{offspringMean}(f)\), let \(Z_0 \in \mathbb{R}\) (the argument start) be the initial population and \(n \in \mathbb{N}\) the number of generations. Then
The expected population after n generations.
Intuition Each generation multiplies the expected population by \(m\), so the mean compounds like interest: exponential growth above one, exponential decay below.
Use it when you project expected cases after \(n\) infection generations, expected users after \(n\) referral rounds, or sanity-check a simulated branching process. Treat the mean with care when \(m > 1\): many runs still die out (with probability extinctionOf(f)), so the mean averages a lot of zeros with a few very large populations and describes no typical run.
population-dynamicsforecasting
def expectedSize(f, start, n) ≔ start · offspringMean(f) ^ n
Building blocks
fixedLoop(k, f, s)— Let \(f \in \Delta^{K-1}\) with generating function \(G\), \(k \in \mathbb{N}\) and \(s \in \mathbb{R}\). Then \(\operatorname{fixedLoop}(0, f, s) = s\) and \(\operatorname{fixedLoop}(k, f, s) = \operatorname{fixedLoop}(k - 1, f, G(s))\), so \(\operatorname{fixedLoop}(k, f, s) = G^{k}(s)\). Applies the generating functionktimes starting froms; used byextinctionOf.
Two named chains
The Ehrenfest model has \(N\) molecules in two urns and moves one molecule, chosen uniformly, at each step. Its stationary distribution is Binomial\((N, \tfrac12)\), and the expected time to return to the state with every molecule on one side is \(2^N\).
Pólya’s urn draws a ball, returns it, and adds reinforce more balls of the same colour. The share of red is a martingale: it does not settle at the starting share but at a random, Beta-distributed limit that differs from run to run. With \(r\) red, \(b\) blue and reinforcement \(c\), the limit is Beta\((r/c,\, b/c)\).
ehrenfestUp
ehrenfestUp(left, molecules) · left molecules in the left urn scalar · molecules total N scalar → probability scalar
Let \(N > 0\) (the argument molecules) be the number of molecules and \(\ell \in \{0, 1, \dots, N\}\) (the argument left) the number in the left urn. Then
The probability that the next move increases the left urn, \((N - \ell)/N\).
Intuition Pick one of the \(N\) molecules at random; it jumps to the left urn only if it was on the right, which happens with chance \((N - \ell)/N\). The fuller one side gets, the likelier the next move empties it, so the chain is pulled back toward the middle.
Use it when you build the Ehrenfest transition matrix (from \(\ell\), up with this probability and down with \(\ell/N\)) to study diffusion and equilibration, model a count that reverts to its mean, such as load shared between two servers, or need a test chain with a known answer for stationary (Binomial\((N, \tfrac12)\)) and mixing.
diffusion-modellingstate-transition-modelling
def ehrenfestUp(left, molecules) ≔ (molecules − left) / nonzero(molecules)
recurrence
recurrence(molecules) · molecules total N scalar → expected return time scalar
Let \(N \in \mathbb{N}\) (the argument molecules) be the number of molecules. Then
The expected number of steps to return to the state with all molecules in one urn.
Intuition Kac’s lemma: the mean return time to a state is one over its stationary probability, and “all on one side” has probability \(2^{-N}\). The gas is certain to unmix again, but only after \(2^N\) steps, which is why nobody ever sees it happen.
Use it when you estimate how rarely a symmetric system revisits an extreme state, or explain why a reversible process still looks irreversible. The lemma applies to any chain: the expected return time to state \(i\) is 1 / stationary(p) at \(i\). The count grows so fast that it overflows a double to infinity beyond about \(N = 1023\).
hitting-timeslong-run-behaviourrare-events
def recurrence(molecules) ≔ 2 ^ molecules
urnShare
urnShare(red, blue) · red red balls scalar · blue blue balls scalar → share of red scalar
Let \(r \ge 0\) (the argument red) and \(b \ge 0\) (the argument blue) be the numbers of red and blue balls. Then
The current share of red balls, which is also the probability that the next draw is red.
Intuition The share of red balls is also the chance the next draw is red, and every red draw raises it: the rich get richer. With one ball added per draw, it is exactly the predictive probability of a Beta–Bernoulli model whose prior counts are the initial balls.
Use it when you read the current probability in a self-reinforcing process: market share under network effects, the popularity of an option that gets recommended more the more it is chosen, a Bayesian success rate updated with pseudo-counts. With no balls at all the guard makes the share 0 rather than ½. To make a draw from it use urnDrawn.
reinforcement-dynamicsbayesian-updating
def urnShare(red, blue) ≔ red / nonzero(red + blue)
urnDrawn
urnDrawn(red, blue, u) · red red balls scalar · blue blue balls scalar · u uniform draw scalar → drew red 1 or 0
Let \(r, b \ge 0\) (the arguments red and blue) be the ball counts and \(u \in [0, 1)\) a uniform draw. Then
1 when the draw u picks a red ball, 0 when it picks blue.
Intuition A coin weighted by the current share of red: the uniform draw falls below the share with exactly that probability.
Use it when you want one draw from explicit ball counts, for instance to build your own urn variant (different reinforcement per colour, balls removed rather than added, a third colour handled in two draws). For the standard Pólya urn with its bookkeeping already done, use urnAt.
simulationreinforcement-dynamics
def urnDrawn(red, blue, u) ≔ u < urnShare(red, blue)
urnAt
urnAt(count, red, blue, reinforce, u) · count balls added so far vector 2 (red, blue) · red initial red scalar · blue initial blue scalar · reinforce balls added per draw scalar · u uniform draw scalar → updated added counts vector 2
Let \(a = (a_0, a_1) \in \mathbb{R}^2\) (the argument count) hold the red and blue balls added so far, let \(r, b \ge 0\) (the arguments red and blue) be the initial balls, \(c \ge 0\) (the argument reinforce) the balls added per draw, and \(u \in [0, 1)\) a uniform draw. Let
Then
\[\operatorname{urnAt}(a, r, b, c, u) = a + c\,\big(d\, e_0 + (1 - d)\, e_1\big) = \big(a_0 + c\,d,\;\; a_1 + c\,(1 - d)\big) ,\]where \(e_0, e_1\) are the basis vectors at positions 0 and 1.
One draw of Pólya’s urn. count holds the red and blue balls added so far (starting at zeros); the draw is made from the initial balls plus those, and the drawn colour gains reinforce.
Intuition The urn is the initial balls plus what count has added, and each draw adds reinforce balls of the drawn colour. Early luck therefore gets locked in: a run that starts with a few reds keeps drawing reds, and the share settles at a value that differs from run to run.
Use it when you simulate path dependence as a process step, starting count at zeros: winner-take-most technology adoption, preferential attachment in networks, popularity that feeds on itself. Run many seeds and look at the spread of final shares, which should follow Beta\((r/c,\, b/c)\); concentration tells you in advance how wide it will be.
simulationreinforcement-dynamicspath-dependence
def urnAt(count, red, blue, reinforce, u) ≔ count + reinforce · ((positions(count) = 0) · urnDrawn(red + elementAt(count, 0), blue + elementAt(count, 1), u) + (positions(count) = 1) · (1 − urnDrawn(red + elementAt(count, 0), blue + elementAt(count, 1), u)))
concentration
concentration(red, blue, reinforce) · red initial red scalar · blue initial blue scalar · reinforce balls added per draw scalar → concentration scalar
Let \(r, b \ge 0\) (the arguments red and blue) be the initial balls and \(c > 0\) (the argument reinforce) the balls added per draw. Then
The concentration \((r + b)/c\) of the Beta limit of the red share: the larger it is, the closer the limit stays to the starting share.
Intuition The initial balls counted in units of one draw’s reinforcement: how many “votes” the starting mix has. With many, each draw nudges the share a little and the limit stays near the start; with few, the first draws decide the outcome.
Use it when you want to know how unpredictable an urn’s outcome is before simulating it. With starting share \(s = r/(r+b)\), the limit has mean \(s\) and variance \(s(1-s)/(\text{concentration} + 1)\). Use it to judge how path-dependent a market is, or to choose how strongly a prior should resist early evidence.
uncertainty-quantificationreinforcement-dynamicspath-dependence
def concentration(red, blue, reinforce) ≔ (red + blue) / nonzero(reinforce)
Fitting a chain to a path
An observed path is a T×S matrix of one-hot rows, one row per time step. The transition counts \(n_{ij} = \sum_t [x_{t-1} = i][x_t = j]\) come from one contraction of the path against itself lagged by one step, and the maximum-likelihood transition matrix divides each row of counts by its total:
\[\hat{p}_{ij} = \frac{n_{ij}}{\sum_k n_{ik}}, \qquad \ell = \sum_{i,j} n_{ij} \ln \hat{p}_{ij}\]moveCounts
moveCounts(path) · path observed states T×S one-hot → transition counts S×S
Let \(X \in \{0, 1\}^{T \times S}\) (the argument path) be an observed path whose row \(t\), for \(t = 0, \dots, T-1\) (counted from zero), is the one-hot state at time \(t\). Then
The number of observed moves from each state to each state.
Intuition A tally sheet: for every pair of consecutive observations, put a mark in the cell (from, to). Matching the path against itself shifted by one step makes all the marks at once.
Use it when you turn logs into transition data: page-to-page clicks, status changes of support tickets, credit-rating migrations. Keep the counts, not just the probabilities, to see how much evidence each row has. To combine several observed paths, add their count matrices; do not stack the paths into one matrix, which would count a fake move where one path ends and the next begins.
parameter-estimationevent-countingstate-transition-modelling
def moveCounts(path) ≔ einsum("ti,tj->ij", laggedRows(path), path)
fittedChain
fittedChain(path) · path observed states T×S one-hot → estimated transition matrix S×S
Let \(X \in \{0, 1\}^{T \times S}\) (the argument path) be an observed one-hot path and \(n = \operatorname{moveCounts}(X)\). Then
The maximum-likelihood transition matrix. A state that is never left in the path gets a row of zeros.
Intuition The obvious estimate: of all the times the path left state \(i\), what fraction went to \(j\)? That frequency is also the maximum-likelihood transition matrix.
Use it when you learn a transition model from observed sequences and then ask it questions with stationary, absorption or mixing. Rows backed by few departures are noisy. A state never left, including one seen only as the last observation, gets a row of zeros, so the matrix is not stochastic and stationary does not apply until you fix that row. For smoothing, add pseudo-counts to moveCounts and normalise the rows yourself.
parameter-estimationstate-transition-modelling
def fittedChain(path) ≔ einsum("ij,i->ij", moveCounts(path), 1 / nonzero(einsum("ij->i", moveCounts(path))))
chainLikelihood
chainLikelihood(path) · path observed states T×S one-hot → log-likelihood scalar
Let \(X \in \{0, 1\}^{T \times S}\) (the argument path) be an observed one-hot path, \(n = \operatorname{moveCounts}(X)\) and \(\hat{p} = \operatorname{fittedChain}(X)\). Then
The maximised log-likelihood of the path under its fitted chain, conditional on the first state. The guard makes a pair never observed contribute \(0 \cdot \ln 10^{-300} = 0\).
Intuition How well the fitted chain explains its own path: the log-probability of every observed move, added up. It is never positive, and closer to zero means the moves were more predictable.
Use it when you compare models of the same sequence, such as a Markov chain against independent draws from the state frequencies in a likelihood-ratio test. The value is maximised on the very path it scores, so a richer model always looks at least as good; charge the chain for its \(S(S-1)\) free parameters (AIC, BIC). It cannot score a matrix you supply, since it always uses fittedChain.
likelihood-evaluationmodel-selection
def chainLikelihood(path) ≔ einsum("ij,ij->", moveCounts(path), ln(nonzero(fittedChain(path))))
Building blocks
laggedRows(path)— Let \(X \in \mathbb{R}^{T \times S}\) with rows at times \(t = 0, \dots, T-1\). Then \(\operatorname{laggedRows}(X)_{tj} = \sum_{s=0}^{T-1} [\, t - s = 1 \,]\, X_{sj}\), so row \(t\) is row \(t - 1\) of \(X\) for \(t \ge 1\) and row \(0\) is zeros. The path shifted back one step: row t holds the state at t − 1, and the first row is zeros.
Examples
Where a chain settles
open Prelude.Chains
fn longRun(p) ↦ stationary(p)
fn afterStep(start, p) ↦ afterOne(start, p)
fn relaxTime(p) ↦ mixing(p)
fn balanceGap(p) ↦ detailedBalance(p, stationary(p))
using Examples;
var p = new double[,] { { 0.9, 0.1 }, { 0.5, 0.5 } };
var pi = Markov.LongRun(p);
var next = Markov.AfterStep(new[] { 1d, 0d }, p);
var relax = Markov.RelaxTime(p);
var gap = Markov.BalanceGap(p);
Output (verified, rounded to 6 decimal places):
stationary = [0.833333, 0.166667]
afterOne = [0.9, 0.1]
mixing = 1.666667
balanceGap = 0
Balance between the two states requires \(0.1\,\pi_0 = 0.5\,\pi_1\), so \(\pi = (5/6, 1/6)\). The eigenvalues of \(P\) are 1 and 0.4, giving a relaxation time of \(1/0.6\). Every two-state chain is reversible, so the detailed-balance gap is zero.
Absorption in a short gambler’s ruin
States 0 and 3 are absorbing; from states 1 and 2 the chain moves one step left or right with probability ½.
open Prelude.Chains
fn stepsToAbsorb(q) ↦ expectedSteps(q)
fn whereAbsorbed(q, r) ↦ absorption(q, r)
using Examples;
var q = new double[,] { { 0, 0.5 }, { 0.5, 0 } }; // 1→1, 1→2 / 2→1, 2→2
var r = new double[,] { { 0.5, 0 }, { 0, 0.5 } }; // 1→0, 1→3 / 2→0, 2→3
var steps = Markov.StepsToAbsorb(q);
var ending = Markov.WhereAbsorbed(q, r);
Output (verified, rounded to 6 decimal places):
steps = [2, 2]
absorption = shape [2, 2] values [0.666667, 0.333333, 0.333333, 0.666667]
\(I - Q\) has inverse \(\tfrac{4}{3}\begin{pmatrix}1 & \tfrac12 \\ \tfrac12 & 1\end{pmatrix}\), whose row sums are 2. Starting from state 1, the chain ends at 0 with probability 2/3.
Extinction of a branching process
Each individual has no children with probability ¼, one with probability ¼ and two with probability ½.
open Prelude.Chains
fn dieOut(f) ↦ extinctionOf(f)
using Examples;
var extinction = Markov.DieOut(new[] { 0.25, 0.25, 0.5 });
Output (verified):
extinction = 0.5
The fixed points of \(G(s) = 0.25 + 0.25 s + 0.5 s^2\) are \(s = \tfrac12\) and \(s = 1\); the mean is 1.25, so extinction is not certain and the answer is the smaller root.
One categorical draw
open Prelude.Chains
fn pickState(v, u) ↦ landedOn(v, u)
using Examples;
var chosen = Markov.PickState(new[] { 0.2, 0.5, 0.3 }, 0.6);
Output (verified):
landed = [0, 1, 0]
The running totals are 0.2, 0.7 and 1.0, and 0.6 falls in the second interval \([0.2, 0.7)\).
Fitting a chain to an observed path
The path visits states 0, 0, 0, 0, 1, 1, 0.
open Prelude.Chains
fn fitMoves(path) ↦ fittedChain(path)
using Examples;
var path = new double[,] { { 1, 0 }, { 1, 0 }, { 1, 0 }, { 1, 0 }, { 0, 1 }, { 0, 1 }, { 1, 0 } };
var fitted = Markov.FitMoves(path);
Output (verified):
fitted = shape [2, 2] values [0.75, 0.25, 0.5, 0.5]
From state 0 the path moves to 0 three times and to 1 once; from state 1 it moves to 1 once and to 0 once.