Interactions: Complexity, Stochasticity, Extrema and Rare Events

Contents

• Introduction and overview
• 0. Review: Markov chains, exponential clocks and branching
• 1. Interacting particle systems: Markov processes, semigroups and generators
• 2. The contact process, the voter model and the exclusion process
• 3. Stochastic Ising models
• 4. Complex systems, emergence and self-organized criticality
• 5. Agent-based models: segregation and epidemics
• 6. Agent-based models: markets, socio-technical systems and avalanches
• 7. Markov chain aggregation
• 8. Inference in interacting systems: structure, laws and causes
• 9. Rare events: conspiracy versus catastrophe
• 10. Higher-dimensional extrema: convex hulls and combinatorial lemmas
• 11. Heavy tails in stochastic gradient descent
• 12. Complements for self-study
• 13. Research directions and open problems
• Exercises
• Projects and assessment
• References

Introduction and overview

Audience: a student with a good M1 in probability, statistics or physics who wants to model systems made of many interacting units, simulate them, learn their rules from data, and understand the extreme events they produce.

Empirical first: place two kinds of agents at random on a 40 × 40 grid and let every agent move whenever fewer than 30% of its neighbours are of its own kind. Nobody asks for more than 30%. Yet after ten rounds of moves everyone is satisfied, and the average share of same-kind neighbours has climbed from 0.50 to 0.74 (Schelling’s model, code in ch. 5). Mild individual rules, strong collective pattern: this gap between the micro and the macro level is what the course is about.

The course studies systems whose behaviour comes from the interactions between their parts. It starts with interacting particle systems, the probabilistic theory of such systems, and moves on to agent-based models of societies, epidemics and markets.

It then asks how to infer interactions and causal effects from data, and closes with the extreme and rare events that interactions generate, from avalanches to heavy tails in stochastic gradient descent.

It shares tools with Introduction to Probabilistic Graphical Models and Deep Generative Models, Graphical Models: Discrete Inference and Learning, Graphs in Machine Learning and Statistical Learning with Extreme Values.

What you will be able to do

• Write the generator of an interacting particle system from its flip or jump rates, and use it to find invariant measures and martingales.
• Analyse the contact process, the voter model and the exclusion process with couplings, duality and generator computations.
• Build Metropolis, heat-bath (Glauber) and Kawasaki dynamics for the Ising model and prove that they are reversible.
• Design an agent-based model, simulate it exactly with the Gillespie algorithm, and use it as an in silico experiment.
• Explain emergence and self-organized criticality on the sandpile and on critical branching processes.
• Reduce an agent-based model to a smaller Markov chain by lumping, and compute absorption probabilities and times.
• Recover an interaction network with penalised pseudo-likelihood, and learn an interaction law by least squares from trajectories.
• Separate correlation from causation in data produced by an agent-based model, and adjust for confounding.
• Tell a light-tailed rare event (conspiracy) from a heavy-tailed one (catastrophe), and compute their probabilities.
• Compute the mean perimeter of the convex hull of a random walk with a combinatorial lemma, and the tail index of stochastic gradient descent with the Kesten–Goldie theorem.

Official objectives and syllabus (2026–2027)

The course description on master-mva.com starts from the observation that interactions are part of many natural, social, artificial and virtual systems. Their modelling has long challenged the sciences, from the three-body problem and chaos to feedback loops in climate, particle systems, confounding variables, cities, markets and opinion dynamics.

It states that data and new theoretical and computational methods make stochastic modelling, simulation and machine learning relevant to such systems. The course introduces the general question of modelling systems of interactions: complexity, emergence and self-organized criticality, and agent-based models (ABMs).

ABMs are presented as in silico controlled experiments on complex real systems, and as a means to test statistical and analytical methods. The course examines how to learn the structure and content of interactions and causal relationships with ABMs, and the extreme and rare events that interactions create in stochastic processes, ABMs and other algorithms.

The description announces ten 2.5-hour lectures in four blocks, listed below with their references. Its “Thèmes abordés” section places the course with respect to other MVA courses. The Université Paris-Saclay course detail gives only the assessment type and no programme of its own.

Every official topic is taught on this page. The extensions beyond the syllabus are the review of ch. 0, the self-study complements of ch. 12 (mean-field limits, the hydrodynamics of the exclusion process, rare-event simulation, the universality of persistence) and the research directions of ch. 13.

Where this course fits in the programme

Official prerequisites: probability, linear algebra, basic statistics and programming (for example C, Python or R). Some training in statistical physics is a bonus and not a prerequisite. Ch. 0 reviews the Markov chain and branching-process facts that the course uses.

• Refreshers — Probability and Statistics (Markov chains, Poisson processes, limit theorems), Linear Algebra, Analysis and Optimization (matrix exponentials, least squares) and Scientific Python (vectorised simulation).
• Semester 1 courses that prepare it — Introduction to Probabilistic Graphical Models and Deep Generative Models (Markov chains, MCMC, network models, named by the official description), Computational Statistics (MCMC), Stochastic Calculus in Machine Learning (diffusions and MCMC) and Statistical Learning with Extreme Values (regular variation). Mathematical Methods for Neuroscience covers the McKean–Tanaka–Sznitman theory of interacting particles (Theorem 12.4 here) for neuron models.
• Related courses of semester 2 — Graphical Models: Discrete Inference and Learning and Graphs in Machine Learning, both named by the description. The description also presents the course as complementing a course by S. Mallat; in this catalogue S. Mallat teaches Diffusion and Multiscale Models in Artificial Intelligence.
• Other related courses of semester 2 — Large Random Matrices and Applications to Machine Learning (the Tracy–Widom law, an extreme of a strongly interacting system), Stopping Times and Online Algorithms (stopping times), Learning Physical Phenomena: Metamodelling and Uncertainty Quantification (surrogates of expensive simulators) and Biostatistics (causal inference).
• After the course — the Research Internship: complex systems, computational social science, epidemiology and the theory of stochastic optimisation all use these tools.

Course outline

• 0. Review — discrete and continuous-time Markov chains, competing exponential clocks, the Gillespie algorithm, Galton–Watson branching.
• 1. Interacting particle systems — Markov semigroups and generators, spin systems, the graphical construction, attractiveness and duality.
• 2. Contact, voter and exclusion processes — survival and extinction, consensus and coalescing random walks, invariant product measures.
• 3. Stochastic Ising models — Gibbs measures, Metropolis, heat-bath (Glauber) and Kawasaki dynamics, the Curie–Weiss phase transition.
• 4. Complex systems, emergence, self-organized criticality — generative explanation, scale invariance, critical branching, the sandpile.
• 5. Segregation and epidemics — Schelling’s model, stochastic epidemics, outbreak probability and final size, public health modelling.
• 6. Markets, socio-technical systems, avalanches — herding, heavy-tailed returns, delay propagation and the critical buffer.
• 7. Markov chain aggregation — lumpability, symmetry, exact consensus times.
• 8. Inference — interaction networks from samples, interaction laws from trajectories, causal effects from simulated data.
• 9. Conspiracy versus catastrophe — large deviations for light tails, the single big jump for heavy tails, where heavy tails come from.
• 10. Higher-dimensional extrema — a combinatorial lemma, Kac’s formula, the convex hull of random walks and Lévy processes.
• 11. Heavy tails in SGD — random linear recursions, the Kesten–Goldie theorem, the tail index of stochastic gradient descent.
• 12. Complements for self-study — mean-field limits, hydrodynamics of the exclusion process, stochastic resetting, rare-event simulation.
• 13. Research directions and open problems — a self-study tour of current work.

Session plan

24 h in class = 10 sessions (8 × 2.5 h + 2 × 2 h). The volume of 24 h of lectures comes from the Université Paris-Saclay programme table. master-mva.com announces ten 2.5-hour lectures (25 h) without exercise classes, so the plan keeps ten lectures and shortens the first and the last to 2 h. The split of the content into sessions is a proposal: neither source details the lectures.

The order follows the four official blocks: particle systems (sessions 1–4), agent-based models (5–7), analysis and inference (8–9), extrema and rare events (9–10). Level 1 = foundations, 2 = core models and results, 3 = advanced and research-level material; the level never goes down.

Self-study, outside the 24 h (about 7 h in total): ch. 12 Complements for self-study, with Ex. 16 (≈ 4 h) · ch. 13 Research directions and open problems (≈ 3 h).

How to study this module

• Read ch. 0 before session 1 if continuous-time Markov chains feel distant: every later model is one, and the proofs use its facts without warning.
• After each session, do the exercises of its row, then run the code of the chapter and change one parameter (a rate, a threshold, a buffer, a step size) to watch the collective behaviour change.
• Treat every simulation as an experiment: state what the theory predicts before running it, and compare to two decimals.
• The code uses NumPy, SciPy and scikit-learn only, with synthetic data; each block runs in a few seconds on a laptop.
• Start a project after session 7: the subjects of the last section use the models of ch. 2–6 and the tools of ch. 7–11.

0. Review: Markov chains, exponential clocks and branching

Every model of this course is a Markov chain on a large state space, run in continuous time by exponential clocks. This chapter collects the facts used later, with short proofs.

The Probability and Statistics refresher covers the same ground at a slower pace.

Discrete-time Markov chains

Two standard theorems complete the picture. If P is irreducible and aperiodic, \pi is unique and P^t(x, \cdot) \to \pi from every x. The ergodic theorem adds that time averages converge to averages under \pi.

These facts justify Markov chain Monte Carlo (ch. 3) and long simulations of agent-based models.

Continuous time: exponential clocks

By Lemma 0.3, a state with several possible transitions behaves as a set of independent clocks, one per transition. This gives an exact simulation method.

import numpy as np
from scipy.stats import poisson

rng = np.random.default_rng(0)
a, b = 5.0, 1.0                 # arrival rate, death rate per individual
n, t, t_end = 0, 0.0, 20000.0
occupation = np.zeros(40)       # time spent in each state
while t < t_end:
    rates = np.array([a, b * n])            # birth, death
    total = rates.sum()
    dt = rng.exponential(1 / total)         # waiting time ~ Exp(total rate)
    occupation[min(n, 39)] += dt
    t += dt
    n += 1 if rng.random() < rates[0] / total else -1   # which clock rang
emp = occupation / occupation.sum()
for k in [2, 5, 8]:
    print(f"state {k}: time fraction {emp[k]:.4f}   Poisson(a/b) {poisson.pmf(k, a / b):.4f}")

The output is 0.0863, 0.1752 and 0.0656 against 0.0842, 0.1755 and 0.0653.

Branching processes

Branching processes describe the early phase of an epidemic (ch. 5), avalanches (ch. 4) and herding clusters (ch. 6). In each case the critical value m = 1 separates small events from macroscopic ones.

1. Interacting particle systems: Markov processes, semigroups and generators

An interacting particle system is a Markov process on configurations: every site of a graph carries a state, which changes at a rate depending on the states of nearby sites. This chapter sets up the framework (semigroups, generators, construction) and two tools used throughout: monotone coupling and duality.

The standard references are Liggett (2005) and, for the social-dynamics reading of the same models, Aldous (2013).

Configurations, semigroups and generators

On an infinite graph the sum defining L f is infinite, and existence is a theorem.

The graphical construction

For finite-range systems there is a direct construction (Harris, 1972) that also couples all initial configurations. Let \bar c bound the rates. Give every site x two independent Poisson clocks of rate \bar c, an up clock and a down clock, and attach an independent uniform U to every ring.

At a ring of the up clock of x, a configuration with \eta(x) = 0 flips to 1 if U < c(x, \eta)/\bar c. Down rings act in the same way on sites in state 1. By thinning, each site flips at the right rate. Infinite graphs need one more argument, given in ch. 12.

Duality

Duality turns a question about an infinite configuration into one about a few dual particles. In ch. 2, the voter model is dual to coalescing random walks and the contact process is dual to itself.

The code below writes the generator of the contact process of ch. 2 on a cycle of 4 sites as a 16 \times 16 matrix and computes P_t = e^{tL}. The empty configuration is absorbing, and the probability of having reached it from the full configuration grows to 1.

import itertools
import numpy as np
from scipy.linalg import expm

# Contact process on the cycle with 4 sites: the generator as a 16 x 16 matrix.
S, lam = 4, 1.5
nbrs = {x: [(x - 1) % S, (x + 1) % S] for x in range(S)}
states = list(itertools.product([0, 1], repeat=S))
index = {s: i for i, s in enumerate(states)}
Q = np.zeros((2**S, 2**S))
for s in states:
    for x in range(S):
        rate = 1.0 if s[x] == 1 else lam * sum(s[y] for y in nbrs[x])  # c(x, eta)
        if rate > 0:
            flipped = list(s); flipped[x] = 1 - s[x]
            Q[index[s], index[tuple(flipped)]] += rate
    Q[index[s], index[s]] = -Q[index[s]].sum()

full, empty = index[(1,) * S], index[(0,) * S]
for t in [1.0, 5.0, 20.0]:
    P = expm(t * Q)                                  # P_t = exp(tQ)
    print(f"t={t:5.1f}: rows sum to 1: {np.allclose(P.sum(1), 1)}, "
          f"P_t(full, empty) = {P[full, empty]:.4f}")

The output gives 0.0555, 0.4494 and 0.9301 at times 1, 5 and 20: every finite contact process dies out, and the question of ch. 2 is how long it takes.

2. The contact process, the voter model and the exclusion process

Three models carry most of the theory of interacting particle systems: an infection that spreads and heals, opinions that are copied, and particles that cannot overlap. Each one shows a different tool at work: coupling and a generator bound, a martingale and duality, and an explicit invariant measure.

The contact process

Better bounds are known: 1/(2d-1) \leq \lambda_c \leq 2/d (Liggett, 2005, ch. VI). On \mathbb{Z}, simulations give \lambda_c \approx 1.649. The graphical construction also shows that the contact process is self-dual: reversing time and arrows maps it onto itself.

On the complete graph with N sites and infection rate \lambda / N per pair, the density of infected sites follows, as N \to \infty, the mean-field equation \dot\rho = \lambda \rho (1 - \rho) - \rho . The endemic state \rho^\ast = 1 - 1/\lambda exists and is stable if and only if \lambda > 1. On any finite graph the process still dies out almost surely, but above criticality only after a time that grows very fast with the size: a metastable state.

The voter model

A hub counts as much as all its leaves. On a star with five leaves and a path of three more sites hanging from a leaf, a single opinion-1 voter at the centre wins with probability 5/16 \approx 0.31, not 1/9:

import numpy as np

rng = np.random.default_rng(2)
# A star (centre 0, leaves 1..5) with a path 5-6-7-8 hanging from leaf 5.
edges = [(0, k) for k in range(1, 6)] + [(5, 6), (6, 7), (7, 8)]
n = 9
nbrs = [[] for _ in range(n)]
for u, v in edges:
    nbrs[u].append(v); nbrs[v].append(u)
deg = np.array([len(a) for a in nbrs])

eta0 = np.zeros(n, dtype=int); eta0[0] = 1          # only the centre holds opinion 1
wins, trials = 0, 4000
for _ in range(trials):
    eta = eta0.copy()
    while 0 < eta.sum() < n:
        x = rng.integers(n)                          # every voter rings at rate 1
        eta[x] = eta[nbrs[x][rng.integers(deg[x])]]  # copy a uniform neighbour
    wins += eta[0]
print(f"P(consensus on 1): simulated {wins / trials:.3f}, "
      f"degree-weighted prediction {deg @ eta0 / deg.sum():.3f}, "
      f"unweighted fraction {eta0.mean():.3f}")

The output is 0.313 against the prediction 0.312.

On \mathbb{Z}^d, two lineages meet surely if and only if the random walk is recurrent: opinions form ever larger clusters in d \leq 2 and coexist in d \geq 3 (Theorem 12.2).

The exclusion process

The particles interact, yet one-point means move like independent random walks. For a doubly stochastic p the measures \nu_\rho stay invariant (Liggett, 2005, ch. VIII). The totally asymmetric version and its hydrodynamic limit are in ch. 12.

3. Stochastic Ising models

The Ising model is the simplest system in which interactions alone create a phase transition. Its stochastic dynamics are both models of physical relaxation and Markov chain Monte Carlo samplers, the MCMC of Introduction to Probabilistic Graphical Models.

Gibbs measures and reversible dynamics

Kawasaki dynamics is an exclusion process (ch. 2) with interactions between particles, a lattice gas. The relocation moves of Schelling’s model (ch. 5) are of the same kind.

The Curie–Weiss phase transition

The dynamics follow the landscape f_\beta. Above the critical temperature, heat-bath dynamics on the complete graph mixes in order N \log N steps (Levin, Peres and Wilmer, 2009, ch. 15). Below it, going from +m^\ast to -m^\ast means crossing the barrier N[f_\beta(0) - f_\beta(m^\ast)]: the time is exponential in N, a metastable state again.

On the square lattice with J = 1, Onsager’s critical temperature is T_c = 2 / \ln(1 + \sqrt 2) \approx 2.269. Below it the spontaneous magnetisation is (1 - \sinh(2\beta)^{-4})^{1/8} (Yang, 1952), which equals 0.986 at T = 1.5. The Metropolis simulation below updates the two colours of a checkerboard in turn, because spins of one colour do not interact.

import numpy as np

rng = np.random.default_rng(3)
L = 32
i, j = np.indices((L, L))
colours = [(i + j) % 2 == 0, (i + j) % 2 == 1]       # checkerboard: no two neighbours updated together

def sweep(s, beta):
    for mask in colours:
        nb = np.roll(s, 1, 0) + np.roll(s, -1, 0) + np.roll(s, 1, 1) + np.roll(s, -1, 1)
        dH = 2 * s * nb                               # energy change of a flip (J = 1, h = 0)
        accept = rng.random(s.shape) < np.exp(-beta * dH)   # Metropolis: min(1, e^{-beta dH})
        s[mask & accept] *= -1

for T in [1.5, 2.27, 3.5]:
    s = np.ones((L, L), dtype=int)
    m = []
    for t in range(3000):
        sweep(s, 1 / T)
        if t >= 1000:
            m.append(abs(s.mean()))
    print(f"T = {T:4.2f}: mean |magnetisation| = {np.mean(m):.3f}")

The output is 0.986, 0.636 and 0.063: ordered, critical with large fluctuations, and disordered. Random couplings of both signs give spin glasses (ch. 12).

4. Complex systems, emergence and self-organized criticality

What is a complex system, and what does it mean to explain one of its collective patterns? This chapter gives the working answer of the course: grow the pattern from interacting micro-rules, and recognise the signatures of criticality.

Emergence and generative explanation

Complex systems combine many units, an interaction network, feedback, heterogeneity and no central controller. Sethna (2021) develops their statistical-physics side: order parameters, phase transitions and scale invariance. Krakauer (2023) argues that complexity science and machine learning should be unified: both try to compress high-dimensional interacting systems into predictive models.

Epstein’s generative social science holds that explaining a macro-regularity means growing it in silico from micro-rules: “If you didn’t grow it, you didn’t explain it” (Epstein, 1999).

An ABM is also a laboratory: changing one rule while everything else, random numbers included, stays fixed is a controlled experiment that the real system rarely allows, such as closing schools (ch. 5) or randomising a treatment (ch. 8).

Critical branching: the mean-field avalanche

At a critical point fluctuations have no characteristic size, and event sizes follow power laws (ch. 12 shows why scale invariance forces them). An avalanche in which each unstable unit destabilises a random number of others is a branching process (ch. 0). Its total size is explicit for Poisson offspring.

Self-organized criticality: the sandpile

Nothing is tuned, yet the stationary state is critical. The drive is slow (one grain at a time), relaxation is fast, and grains leave only at the boundary. If avalanches were subcritical, grains would accumulate and make them larger; if supercritical, large avalanches would empty the grid.

In mean field, a toppling destabilises each of 4 neighbours with probability p_3, the fraction of sites at height 3. The branching ratio 4 p_3 is driven to 1, the critical case of Theorem 4.3, with exponent 3/2. In two dimensions the exponent differs and is not known exactly.

import numpy as np

rng = np.random.default_rng(4)
L = 30
z = np.zeros((L, L), dtype=int)          # heights; grains falling off the border are lost

def add_grain(z):
    """Drop one grain at a random site, relax, return the avalanche size (number of topplings)."""
    x, y = rng.integers(L, size=2)
    z[x, y] += 1
    stack, size = [(x, y)] if z[x, y] >= 4 else [], 0
    while stack:
        x, y = stack.pop()
        if z[x, y] < 4:
            continue
        z[x, y] -= 4; size += 1
        for u, v in ((x + 1, y), (x - 1, y), (x, y + 1), (x, y - 1)):
            if 0 <= u < L and 0 <= v < L:
                z[u, v] += 1
                if z[u, v] >= 4:
                    stack.append((u, v))
        if z[x, y] >= 4:
            stack.append((x, y))
    return size

sizes = np.array([add_grain(z) for _ in range(40000)])[10000:]   # discard the transient
print(f"mean height {z.mean():.3f}, fraction of grains causing an avalanche {np.mean(sizes > 0):.3f}")
s = sizes[sizes > 0]
for lo, hi in [(1, 10), (10, 100), (100, 1000)]:
    print(f"P({lo} <= S < {hi} | S > 0) = {np.mean((s >= lo) & (s < hi)):.3f}")
# slope of the log-log histogram on logarithmic bins between 3 and 300 topplings
bins = np.unique(np.logspace(0.5, 2.5, 9).astype(int))
hist, _ = np.histogram(s, bins=bins)
dens = hist / np.diff(bins)
slope = np.polyfit(np.log(np.sqrt(bins[:-1] * bins[1:])), np.log(dens), 1)[0]
print(f"log-log slope of the size density on [3, 300]: {slope:.2f}")

The mean height settles at 2.098, close to the infinite-lattice value 17/8 = 2.125. The avalanche sizes spread evenly over three decades: 41%, 36% and 23% of them fall in [1, 10), [10, 100) and [100, 1000). The fitted slope is -1.06 on this small grid.

5. Agent-based models: segregation and epidemics

Two classic agent-based models, one from social science and one from public health, show how individual rules produce collective outcomes. The second also shows how a model becomes a decision tool.

Schelling’s segregation model

With L = 40, v = 0.1 and the tolerant threshold T = 0.3, the mean share of same-type neighbours rises from 0.50 to 0.74, and all agents are satisfied after ten rounds:

import numpy as np

rng = np.random.default_rng(1)
L, vacancy, T = 40, 0.1, 0.3          # grid side, vacant fraction, tolerance threshold
grid = rng.choice([0, 1, 2], size=(L, L), p=[vacancy, 0.45, 0.45])  # 0 = vacant

def like_fraction(g):
    """For each cell: fraction of occupied neighbours (8-neighbourhood, torus) of the same type."""
    same = np.zeros(g.shape); occ = np.zeros(g.shape)
    for dx in (-1, 0, 1):
        for dy in (-1, 0, 1):
            if dx or dy:
                h = np.roll(np.roll(g, dx, 0), dy, 1)
                same += (h == g) & (h > 0); occ += h > 0
    return np.where(occ > 0, same / np.maximum(occ, 1), 1.0)

def segregation(g):
    f = like_fraction(g)
    return f[g > 0].mean()

print(f"initial mean like-neighbour fraction: {segregation(grid):.3f}")
for sweep in range(60):
    f = like_fraction(grid)
    unhappy = np.argwhere((grid > 0) & (f < T))
    if len(unhappy) == 0:
        break
    for i, j in unhappy[rng.permutation(len(unhappy))]:
        empty = np.argwhere(grid == 0)
        k = empty[rng.integers(len(empty))]
        grid[k[0], k[1]], grid[i, j] = grid[i, j], 0
print(f"rounds of moves: {sweep}, unhappy agents left: {len(unhappy)}")
print(f"final mean like-neighbour fraction  : {segregation(grid):.3f}")

A smoother version of the rules has an explicit stationary law. Let the utility u of an agent be its number of same-type neighbours, and let \Phi count the adjacent pairs of same-type agents.

Each agent maximises its own utility, and the collective state samples a Gibbs measure for the potential \Phi. For a lattice gas with two species, this is a long-range version of Kawasaki dynamics (Proposition 3.4): agents jump to any vacant cell, not only to a neighbouring one. Large \beta means phase separation. Blume (1993) gives the general statement for potential games.

Stochastic epidemics

In the stochastic SIR model, a population of N individuals is split into susceptible, infected and removed. Each infective meets a uniform individual at rate \beta and infects it if it is susceptible; it is removed at rate \gamma. The total rates are \beta S I / N and \gamma I, and the basic reproduction number is R_0 = \beta / \gamma.

For R_0 = 2.5, z = 0.893. Infections stop growing once the immune fraction reaches the herd-immunity threshold 1 - 1/R_0 = 0.6, but the epidemic overshoots it. Reducing transmission by 40% (R_0 = 1.5) brings z down to 0.583.

The simulation of the jump chain below has N = 1000 and one initial case. It finds 41.0% of minor outbreaks, against 1/R_0 = 0.40, and a mean major outbreak of 0.891, against 0.893.

import numpy as np
from scipy.optimize import brentq

rng = np.random.default_rng(5)
N, R0, gamma = 1000, 2.5, 1.0
beta = R0 * gamma

def final_size(S, I):
    """Stochastic SIR (well-mixed): only the jump chain matters for the final size."""
    while I > 0:
        p_inf = beta * S * I / N / (beta * S * I / N + gamma * I)
        if rng.random() < p_inf:
            S, I = S - 1, I + 1
        else:
            I -= 1
    return 1 - S / N

z = np.array([final_size(N - 1, 1) for _ in range(400)])
z_det = brentq(lambda z: z - 1 + np.exp(-R0 * z), 1e-6, 1)   # z = 1 - exp(-R0 z)
print(f"minor outbreaks (< 10%): {np.mean(z < 0.1):.3f}   branching prediction 1/R0 = {1 / R0:.3f}")
print(f"mean size of major outbreaks: {z[z >= 0.1].mean():.3f}   final-size equation: {z_det:.3f}")
print(f"herd-immunity threshold 1 - 1/R0 = {1 - 1 / R0:.3f}")

Public health modelling

Adam (2020) describes the simulations that drove the response to COVID-19: individual-based models of whole countries, with households, schools and workplaces, used to compare interventions. They add what the equations above average out: contact structure, heterogeneity, targeted interventions, and the randomness of the early phase.

Such a model carries three uncertainties: parameters (R_0, the infectious period), stochasticity (Proposition 5.3) and structure (the model itself). Hence ensembles of runs with reported intervals, sensitivity analyses, calibration on data, and scenarios compared with common random numbers, so that a difference between runs measures the intervention and not the noise.

6. Agent-based models: markets, socio-technical systems and avalanches

When agents copy each other or depend on each other, small causes can have large collective effects. This chapter studies herding in markets and the propagation of delays in socio-technical systems, and finds in both a critical point.

Herding: Kirman’s model

Kirman (1993) observed ants that, offered two identical food sources, crowd at one of them and switch from time to time. In his model, N agents choose between two options. An agent switches spontaneously at rate \varepsilon, and each agent of the other option recruits it at rate h/(N - 1).

With k agents on option 1, the count is a birth–death chain with rates k \to k + 1 : \ (N - k)\Big(\varepsilon + \frac{h k}{N - 1}\Big), \qquad k \to k - 1 : \ k \Big(\varepsilon + \frac{h (N - k)}{N - 1}\Big).

Herding wins, a < 1, when spontaneous switching is rarer than recruitment by one given agent. Read as a market, the majority persists for long periods and then flips, in bursts.

import numpy as np
from scipy.stats import betabinom

rng = np.random.default_rng(6)
N, h = 100, 1.0
for eps in [0.004, 0.04]:
    a = eps * (N - 1) / h                      # the parameter of Theorem 6.1
    k, time_in = N // 2, np.zeros(N + 1)
    for _ in range(400000):                     # jump chain + exponential holding times
        up = (N - k) * (eps + h * k / (N - 1))
        down = k * (eps + h * (N - k) / (N - 1))
        time_in[k] += rng.exponential(1 / (up + down))
        k += 1 if rng.random() < up / (up + down) else -1
    emp = time_in / time_in.sum()
    x = np.arange(N + 1)
    exact = betabinom.pmf(x, N, a, a)
    tv = 0.5 * np.abs(emp - exact).sum()
    print(f"a = {a:.2f}: P(k <= 10 or k >= 90) simulated {emp[(x <= 10) | (x >= 90)].sum():.3f}, "
          f"exact {exact[(x <= 10) | (x >= 90)].sum():.3f}; total variation {tv:.3f}")

For a = 0.40 the crowd spends half of the time beyond 90% on one side (0.520 simulated, 0.496 exact); for a = 3.96, about 1% of the time (0.009 and 0.011).

Herding and heavy-tailed returns

The aggregate demand is a sum of cluster sizes with random signs, a sum of heavy-tailed variables. By ch. 9, its large values come from the largest cluster acting alone: one coordinated block moves the price. Empirical returns have lighter tails than this critical model, so the model shows a mechanism rather than the numbers.

Delays in socio-technical systems

In supply chains, transport networks and production schedules, agents need inputs from others and keep buffers of spare time. Delays propagate along the dependencies, and buffers absorb part of them. Smaller buffers are more efficient, until delays accumulate.

Moran and coauthors (2024) call the resulting phase transition timeliness criticality, with avalanches of delays near it. A stylised mean-field model in the same spirit: agent i waits for K inputs drawn at random at each step, absorbs a buffer B, and adds its own random delay \varepsilon_i(t) \sim \mathrm{Exp}(1): \tau_i(t + 1) = \max\Big(0,\ \max_{j \in J_i(t)} \tau_j(t) - B\Big) + \varepsilon_i(t).

import numpy as np
from scipy.optimize import minimize_scalar

rng = np.random.default_rng(7)
K = 2
# speed of the maximum of a branching random walk with K children and Exp(1) steps
speed = minimize_scalar(lambda th: (np.log(K) - np.log(1 - th)) / th,
                        bounds=(1e-3, 1 - 1e-3), method="bounded").fun
print(f"critical buffer B_c = v_K = {speed:.3f}")

N, T = 5000, 600
for B in [2.0, 2.4, 2.6, 2.8, 3.2]:
    tau = np.zeros(N)
    for t in range(T):
        inputs = rng.integers(N, size=(N, K))                 # fresh random inputs (mean field)
        tau = np.maximum(tau[inputs].max(axis=1) - B, 0) + rng.exponential(1.0, N)
        if t == T // 2:
            tau_half = tau.mean()
    print(f"B = {B:.1f}: growth rate of the mean delay {(tau.mean() - tau_half) / (T / 2):.3f}, "
          f"BRW prediction {max(speed - B, 0):.3f}")

With 5000 agents, the delays grow at rates 0.471 and 0.113 for buffers 2.0 and 2.4, and not at all from 2.6: a finite population cannot explore all K^s paths.

The speed at B = 2 is 2.38, 2.47 and 2.58 for N = 500, 5000 and 50000. It approaches v_K slowly, as the (\ln N)^{-2} correction of Brunet and Derrida (1997) predicts.

7. Markov chain aggregation

An ABM is a Markov chain on S^N, but what we observe and want to predict is a macro-variable: a count, a density, a share. This chapter asks when the macro-variable is itself a Markov chain, and uses the answer to compute exact consensus times (Banisch, 2015).

Lumpability

On the complete graph, an ABM with exchangeable agents is invariant under all permutations of the agents, and its orbits are the vectors of counts. On a general network, the automorphism group of the graph gives a finer lumpable partition.

Exact consensus times

The code solves the linear system exactly for N = 100 and checks it against the agent-level model, without any lumping:

import numpy as np

# Voter model on the complete graph K_N, one random agent updated per step.
# Lumped chain on k = number of agents with opinion 1: p_k = k (N - k) / (N (N - 1)) up and down.
N = 100
k = np.arange(1, N)
p = k * (N - k) / (N * (N - 1))
# Expected absorption time tau_k: p_k (tau_{k+1} - 2 tau_k + tau_{k-1}) = -1, tau_0 = tau_N = 0.
A = np.diag(-2 * p) + np.diag(p[1:], -1) + np.diag(p[:-1], 1)
tau = np.linalg.solve(A, -np.ones(N - 1))
x = k / N
H = -(x * np.log(x) + (1 - x) * np.log(1 - x))       # entropy in nats
for kk in [10, 50]:
    print(f"k = {kk}: exact tau = {tau[kk - 1]:8.0f} steps, N(N-1) H(k/N) = {N * (N - 1) * H[kk - 1]:8.0f}")

# Monte Carlo check from k = 50 on the full (unlumped) agent-level model
rng = np.random.default_rng(8)
times = []
for _ in range(300):
    eta, ones, t = [1] * 50 + [0] * 50, 50, 0
    I = rng.integers(N, size=200000)
    J = ((I + 1 + rng.integers(N - 1, size=200000)) % N).tolist()  # agent i copies a uniform other j
    I = I.tolist()
    while 0 < ones < N:
        i, j = I[t], J[t]
        ones += eta[j] - eta[i]; eta[i] = eta[j]; t += 1
    times.append(t)
print(f"agent-level simulation from k = 50: mean {np.mean(times):.0f} steps (s.e. {np.std(times) / np.sqrt(300):.0f})")

The exact times are 3170 and 6813 steps, the approximation gives 3218 and 6862, and the agent-level simulation gives 7068 \pm 308.

When lumping fails

On the path with three sites, the configurations (1, 0, 0) and (0, 1, 0) both have one opinion-1 voter. From the first, k increases only when site 2 copies site 1, at rate 1/2. From the second, it increases when site 1 or site 3 copies site 2, at total rate 2. The count is not lumpable.

Two remedies exist. The automorphisms of the graph (here the reflection) give a finer partition that is lumpable. Approximate aggregation, such as mean-field or pair approximations, keeps a small macro-chain and accepts a quantified error. Either way, the micro-model provides the ground truth against which the reduced model is tested.

8. Inference in interacting systems: structure, laws and causes

The previous chapters went from rules to patterns. Inference goes back: from observed configurations or trajectories to who interacts with whom, how strongly, and with what causal effect. ABMs play a double role here, as the models to fit and as test benches with a known truth.

Inferring interaction structure

Write the Ising measure with site-dependent couplings, the inverse temperature absorbed: \mu(\sigma) \propto \exp\Big( \sum_{x < y} J_{xy}\, \sigma_x \sigma_y + \sum_x h_x \sigma_x \Big). The interaction graph is the set of pairs with J_{xy} \neq 0.

Neighbourhood selection fits, for each site, an \ell_1-penalised logistic regression of \sigma_x on the other spins, and keeps the edge \{x, y\} when both coefficients are non-zero. This pseudo-likelihood method never computes the partition function.

Correlation is not interaction. On a chain with couplings J, spins at distance 2 have correlation \tanh^2 J without any coupling. With 16 spins on a cycle with 4 chords and J = 0.5, from 3000 samples, neighbourhood selection finds the 20 edges and no false ones, while the 20 most correlated pairs contain only 18 edges.

import numpy as np
from sklearn.linear_model import LogisticRegression

rng = np.random.default_rng(9)
p, n, beta = 16, 3000, 0.5
# true interaction graph: a cycle plus 4 random chords, couplings beta on each edge
E = {(i, (i + 1) % p) for i in range(p)}
while len(E) < p + 4:
    i, j = sorted(rng.choice(p, 2, replace=False)); E.add((i, j))
J = np.zeros((p, p))
for i, j in E:
    J[i, j] = J[j, i] = beta

# n independent samples by Gibbs sampling (heat-bath), n chains run in parallel
s = rng.choice([-1, 1], size=(n, p))
for sweep in range(300):
    for i in range(p):
        field = s @ J[i]
        s[:, i] = np.where(rng.random(n) < 1 / (1 + np.exp(-2 * field)), 1, -1)

# neighbourhood selection: l1-penalised logistic regression of each spin on the others
# (scikit-learn >= 1.8 syntax; older versions write penalty="l1")
J_hat = np.zeros((p, p))
for i in range(p):
    others = np.delete(np.arange(p), i)
    clf = LogisticRegression(l1_ratio=1.0, C=0.05, solver="liblinear").fit(s[:, others], s[:, i])
    J_hat[i, others] = clf.coef_[0] / 2           # logit = 2 * sum_j J_ij s_j
found = {(i, j) for i in range(p) for j in range(i + 1, p) if min(abs(J_hat[i, j]), abs(J_hat[j, i])) > 0.1}
true = {tuple(sorted(e)) for e in E}
C = np.corrcoef(s.T)
top = sorted(((abs(C[i, j]), (i, j)) for i in range(p) for j in range(i + 1, p)), reverse=True)[:len(true)]
print(f"true edges {len(true)}, logistic-l1 finds {len(found & true)} of them and {len(found - true)} false ones")
print(f"thresholding correlations (top {len(true)} pairs) finds {len({e for _, e in top} & true)} true edges")
print(f"mean estimated coupling on true edges: {np.mean([J_hat[i, j] for i, j in true]):.2f} (true {beta})")

The estimated couplings average 0.45 instead of 0.5: the \ell_1 penalty shrinks them. From observed trajectories, the likelihood of a spin system is a product of flip rates and waiting times, and its maximisation is again one regression per site.

Learning interaction laws

Lu, Zhong, Tang and Maggioni (2019) learn \phi from M trajectories observed at L times.

The right-hand side f_\phi(x) is linear in \phi. Minimising over a finite-dimensional space \mathcal{H}_n, such as piecewise constant functions on n bins, is therefore a least-squares problem: \hat\phi = \arg\min_{\varphi \in \mathcal{H}_n} \sum_{m, l, i} \big\| \dot x^{(m)}_i(t_l) - f_\varphi\big(x^{(m)}(t_l)\big)_i \big\|^2 . The error is measured on r \mapsto r \phi(r) in L^2(\rho_T), where \rho_T is the law of pairwise distances along the trajectories.

The pairwise structure is what makes the problem one-dimensional: all pairs of all trajectories inform the same function of one variable.

import numpy as np

rng = np.random.default_rng(10)
N, d, dt, L = 10, 2, 0.05, 20                   # agents, dimension, time step, observations per trajectory
phi = lambda r: np.where(r < 1.0, 1.0, 0.1)     # true interaction kernel (opinion-dynamics type)
edges = np.linspace(0, 3, 16)                   # hypothesis space: piecewise constant on 15 bins

def velocity(x, kernel):
    diff = x[None, :, :] - x[:, None, :]        # diff[i, j] = x_j - x_i
    r = np.linalg.norm(diff, axis=2)
    return (kernel(r)[:, :, None] * diff).sum(1) / N

def learn(M, noise=0.05):
    rows, targets, dists = [], [], []
    for _ in range(M):
        x = rng.uniform(0, 2, size=(N, d))
        for _ in range(L):
            v = velocity(x, phi) + noise * rng.standard_normal((N, d))   # noisy observed velocities
            diff = x[None, :, :] - x[:, None, :]
            r = np.linalg.norm(diff, axis=2)
            b = np.clip(np.digitize(r, edges) - 1, 0, len(edges) - 2)
            for i in range(N):
                feat = np.zeros((d, len(edges) - 1))
                for j in range(N):
                    if j != i:
                        feat[:, b[i, j]] += diff[i, j] / N
                        dists.append(r[i, j])
                rows.append(feat); targets.append(v[i])
            x = x + dt * velocity(x, phi)                                 # Euler step of the true system
    A, y = np.concatenate(rows), np.concatenate(targets)
    c = np.linalg.lstsq(A, y, rcond=None)[0]
    dists = np.array(dists)
    b = np.clip(np.digitize(dists, edges) - 1, 0, len(edges) - 2)
    err = np.sqrt(np.mean(((c[b] - phi(dists)) * dists) ** 2) / np.mean((phi(dists) * dists) ** 2))
    return c, err

for M in [5, 50]:
    c, err = learn(M)
    print(f"M = {M:3d} trajectories: relative L2(rho) error {err:.3f}; "
          f"phi_hat on [0.4, 0.6) = {c[2]:.2f}, on [1.4, 1.6) = {c[7]:.2f}")

The relative error is 0.033 with 5 trajectories and 0.011 with 50; the estimated kernel takes the values 1.02 and 0.12, then 1.01 and 0.10, against 1 and 0.1.

Causal inference with agent-based models

The code of an ABM is a structural causal model of its own variables, and running it with a modified rule is an intervention. Manzo (2022) develops this mechanism-based view of causality in the social sciences.

In the model below, cautious agents both wear masks more often and have fewer contacts. The naive comparison of infection rates doubles the effect of masks; adjusting for caution, or randomising masks in silico, gives the true effect.

import numpy as np

rng = np.random.default_rng(11)
n, p = 200_000, 0.02                         # agents, transmission probability per risky contact

def run(mask_rule):
    caution = rng.random(n) < 0.5            # confounder: cautious agents...
    mask = mask_rule(caution)
    contacts = rng.poisson(np.where(caution, 5, 15))   # ...also have fewer contacts
    p_inf = 1 - (1 - p * np.where(mask, 0.5, 1.0)) ** contacts   # a mask halves transmission
    return caution, mask, rng.random(n) < p_inf

# observational world: cautious agents wear masks more often
C, M, Y = run(lambda c: rng.random(n) < np.where(c, 0.8, 0.2))
naive = Y[M].mean() - Y[~M].mean()
adjusted = sum((Y[M & (C == c)].mean() - Y[~M & (C == c)].mean()) * np.mean(C == c) for c in (0, 1))
# in silico experiment: the same model, masks assigned by a coin flip
C2, M2, Y2 = run(lambda c: rng.random(n) < 0.5)
rct = Y2[M2].mean() - Y2[~M2].mean()
print(f"naive difference {naive:+.4f} | back-door adjusted {adjusted:+.4f} | randomised in silico {rct:+.4f}")

The output is -0.159, -0.082 and -0.080.

Interacting systems add a difficulty: interference. My neighbour’s vaccination protects me, so the outcome of an agent depends on the treatments of others. Hudgens and Halloran (2008) separate the direct effect (an agent’s own treatment, at fixed coverage), the indirect effect (the coverage of others, on untreated agents) and their sum, the total effect. An epidemic ABM is a natural test bench for such designs (Exercise 17).

9. Rare events: conspiracy versus catastrophe

A sum of many random contributions, a total loss, a total delay or an aggregate demand, takes an unusually large value. Either all contributions are a little larger than usual, or one of them is huge. Nair, Wierman and Zwart (2022) call these the conspiracy and the catastrophe principles; the tail of the contributions decides which one holds.

Light tails: conspiracy

Conditioned on the rare event, each X_i is distributed approximately as the tilted law of mean a: all summands conspire (the Gibbs conditioning principle). For exponential variables the conditional law is explicit.

Heavy tails: catastrophe

The simulation compares the two worlds for n = 10, at the level that S_n exceeds with probability 2 \times 10^{-4}:

import numpy as np

rng = np.random.default_rng(12)
n, reps = 10, 500_000
laws = {
    "exponential, mean 1": (rng.exponential(1.0, (reps, n)), lambda x: np.exp(-x)),
    "Pareto alpha = 1.5": (rng.pareto(1.5, (reps, n)) + 1, lambda x: x ** -1.5),  # P(X > x) = x^-1.5
}
for name, (X, tail) in laws.items():
    S, Mx = X.sum(1), X.max(1)
    x = np.quantile(S, 1 - 2e-4)                    # a level exceeded with probability 2e-4
    big = S > x
    print(f"{name:20s}: x = {x:6.1f}, P(S>x)/(n P(X>x)) = {np.mean(big) / (n * tail(x)):10.3g}, "
          f"E[max/S | S > x] = {np.mean(Mx[big] / S[big]):.2f}")

For exponentials (x = 24.8), the single-big-jump estimate n \mathbb{P}(X > x) is wrong by a factor of 10^6, and the largest summand carries 28% of the sum (H_{10}/10 = 0.293). For the Pareto law (x = 1397) the ratio is 1.04 and the largest summand carries 99% of the sum.

Where heavy tails come from

The course has met the other mechanisms: critical branching and percolation clusters (ch. 4 and ch. 6), multiplicative recursions (ch. 11) and extremes (ch. 10). Interactions produce heavy tails by coordinating units, as in herding, or by multiplying fluctuations.

10. Higher-dimensional extrema: convex hulls and combinatorial lemmas

The maximum of a random walk is a one-dimensional extreme. In the plane, the natural extreme object is the convex hull of the trajectory, whose extent in every direction is a maximum. A combinatorial lemma computes its mean perimeter exactly, whatever the law of the steps (Majumdar, Comtet and Randon-Furling, 2010).

A combinatorial lemma and Kac’s formula

Let x_1, \dots, x_n be real numbers with partial sums S_k (S_0 = 0) and running maxima M_k = \max(0, S_1, \dots, S_k).

The j-th cyclic shift is the sequence x_{j+1}, \dots, x_n, x_1, \dots, x_j. Write M^{(j)}_{n-1} for its running maximum after n - 1 steps.

Convex hulls of planar random walks

For standard Gaussian steps, \mathbb{E}\, \|S_k\| = \sqrt{\pi k / 2} and \mathbb{E}\, L_n \approx \sqrt{8 \pi n}. For planar Brownian motion on [0, T], each h(\theta) is the maximum of a one-dimensional Brownian motion, of mean \sqrt{2T/\pi}, so the mean perimeter is \sqrt{8 \pi T}.

import numpy as np
from scipy.spatial import ConvexHull

rng = np.random.default_rng(13)
n, reps = 1000, 2000
per, area = [], []
for _ in range(reps):
    S = np.vstack([[0, 0], np.cumsum(rng.standard_normal((n, 2)), axis=0)])   # planar walk from 0
    h = ConvexHull(S)
    per.append(h.area); area.append(h.volume)    # in 2D, scipy's "area" is the perimeter
k = np.arange(1, n + 1)
spitzer_widom = 2 * np.sum(np.sqrt(np.pi * k / 2) / k)    # 2 sum E|S_k| / k, E|S_k| = sqrt(pi k / 2)
print(f"mean perimeter: simulated {np.mean(per):.1f} +- {np.std(per) / np.sqrt(reps):.1f}, "
      f"Spitzer-Widom {spitzer_widom:.1f}, Brownian sqrt(8 pi n) = {np.sqrt(8 * np.pi * n):.1f}")
print(f"mean area     : simulated {np.mean(area):.0f} +- {np.std(area) / np.sqrt(reps):.0f}, "
      f"Brownian pi n / 2 = {np.pi * n / 2:.0f}")
bnb = 0.5 * sum(np.sum(1 / np.sqrt(i * np.arange(1, n - i + 1))) for i in range(1, n))
print(f"discrete mean-area formula (1/2) sum_(i+j<=n) 1/sqrt(ij) = {bnb:.0f}")

For n = 1000, the simulated mean perimeter is 155.2 \pm 0.7 against 154.9 from Theorem 10.3, and the mean area is 1486 \pm 12 against 1480 from Theorem 10.4. The Brownian values, 158.5 and 1571, are reached only slowly.

Lévy processes

The hull of a Lévy flight grows like T^{1/\alpha}, faster than the diffusive \sqrt T: the rare long jumps dominate the extent, the catastrophe principle in two dimensions. The same tools give the mean home range of n independent foragers, the hull of n Brownian paths (Randon-Furling, Majumdar and Comtet, 2009).

11. Heavy tails in stochastic gradient descent

Stochastic gradient descent is itself a stochastic process in which the iterate interacts with random data, multiplicatively. Even with Gaussian data, its iterates can have power-law tails (Gurbuzbalaban, Simsekli and Zhu, 2021). This chapter explains why with random linear recursions.

Random linear recursions

Consider X_{k+1} = A_k X_k + B_k with i.i.d. pairs (A_k, B_k) of real numbers. If \mathbb{E} \ln |A| < 0 and \mathbb{E} \ln^+ |B| < \infty, the products A_1 \cdots A_k decay exponentially and the series X = \sum_{k \geq 0} A_1 \cdots A_k\, B_{k+1} converges almost surely; its law is the unique stationary law, the solution of X \overset{d}{=} A X + B with (A, B) independent of X.

The function p \mapsto \mathbb{E}\, |A|^p is convex, equal to 1 at p = 0 and at p = \alpha, and below 1 in between. Moments of order p < \alpha are finite and, by Theorem 11.1, those of order p \geq \alpha are not.

SGD on least squares

Take one-dimensional data y = a w^\ast + \varepsilon with a, \varepsilon \sim \mathcal{N}(0, 1), and SGD with step \eta and batches of size b. The error X_k = w_k - w^\ast follows a Kesten recursion: X_{k+1} = \Big(1 - \eta\, \frac{1}{b} \sum_{j \leq b} a_j^2\Big) X_k + \eta\, \frac1b \sum_{j \leq b} a_j \varepsilon_j .

For b = 1, \mathbb{E}(1 - \eta a^2)^2 = 1 - 2\eta + 3\eta^2, which equals 1 at \eta = 2/3. The stationary variance of SGD is infinite for 2/3 \leq \eta < 2.42, although the data are Gaussian and the iteration is stable: \mathbb{E} \ln |A| < 0 up to \eta \approx 2.42. Larger steps and smaller batches give heavier tails; Gurbuzbalaban, Simsekli and Zhu (2021) prove this monotonicity for Gaussian linear regression in any dimension.

The tail index is estimated from the iterates with the Hill estimator (Hill, 1975): with |X|_{(1)} \geq |X|_{(2)} \geq \dots the ordered values, \hat\alpha_k = \Big( \frac1k \sum_{i=1}^{k} \ln \frac{|X|_{(i)}}{|X|_{(k+1)}} \Big)^{-1} .

import numpy as np
from scipy.optimize import brentq

rng = np.random.default_rng(14)

def tail_index(eta, b):
    """Solve E|1 - eta * mean of b squared Gaussians|^alpha = 1 (Kesten-Goldie exponent)."""
    A = np.abs(1 - eta * rng.chisquare(b, 400_000) / b)
    assert np.mean(np.log(A)) < 0, "unstable: E log|A| >= 0"
    return brentq(lambda al: np.mean(A ** al) - 1, 1e-3, 50)

def hill(x, k):
    x = np.sort(np.abs(x))[::-1]
    return 1 / np.mean(np.log(x[:k] / x[k]))

for eta, b in [(0.5, 1), (0.9, 1), (1.2, 1), (1.2, 10)]:
    alpha = tail_index(eta, b)
    # 1-D least squares, y = a * 1 + noise: SGD iterate x_{k+1} = (1 - eta a^2) x_k + eta a y
    x, xs = 0.0, []
    a = rng.standard_normal((200_000, b)); noise = rng.standard_normal((200_000, b))
    for t in range(200_000):
        g = np.mean(a[t] * (a[t] * x - (a[t] + noise[t])))   # minibatch gradient
        x -= eta * g
        if t > 1000 and t % 10 == 0:
            xs.append(x - 1)
    print(f"eta = {eta}, b = {b:2d}: Kesten exponent alpha = {alpha:5.2f}, "
          f"Hill estimate on the iterates {hill(np.array(xs), 200):5.2f}")

With b = 1, the exponent drops from 2.88 to 1.30 and 0.78 as \eta grows from 0.5 to 0.9 and 1.2; batches of 10 bring it back to 5.07 at \eta = 1.2. The Hill estimates (3.32, 1.14, 0.82, 4.34) follow the same order; they are noisy, since the iterates are correlated and 200 order statistics are few.

In dimension d, A_k = I - \eta H_k is a random matrix and the exponent solves a matrix version of the same equation (Kesten, 1973). In decentralised SGD, agents on a network average their iterates with their neighbours and take local steps. Gurbuzbalaban and coauthors (2022) study how the network topology, the step size and the batch size then shape the tail index.

12. Complements for self-study

This chapter collects complements to the lectures: the infinite-volume construction, the voter model on lattices, disordered systems, scale invariance, mean-field limits, the hydrodynamics of the exclusion process, stochastic resetting, rare-event simulation and the universality of persistence. Each section depends only on the chapter it complements.

Constructing infinite systems (complement to ch. 1)

The graphical construction of ch. 1 runs every site with two Poisson clocks of rate \bar c. In infinite volume, the state of a site at time t could in principle depend on infinitely many clock rings. The fix is to work on short time windows.

On [0, t_0], a site interacts only with its finite cluster, so its evolution is a finite computation. Repeating on [t_0, 2t_0], [2t_0, 3t_0], and so on, defines the process for all times, with the same clocks coupling all initial conditions.

The voter model on lattices (complement to ch. 2)

Proof sketch. By duality (Proposition 2.5), \eta_t(x) \neq \eta_t(y) requires that the lineages of x and y have not met by time t, and that their ancestors disagree. The difference of the two lineages is a random walk, recurrent if and only if d \leq 2 (Pólya). In d \geq 3 lineages escape each other with positive probability, and the limit law is read off the coalescing walks.

Disordered systems (complement to ch. 3)

If the couplings J_{xy} of the Ising model are themselves random, of both signs, some loops cannot satisfy all their bonds: the system is frustrated. These spin glasses, the Edwards–Anderson model on a lattice and the Sherrington–Kirkpatrick model on the complete graph, have the same Metropolis and heat-bath dynamics as in ch. 3. They have many metastable states and relax very slowly (Mézard and Montanari, 2009).

The statistical physics of disordered systems studies their typical behaviour over the randomness of the couplings. Learning the couplings from samples, the inverse problem, is the subject of ch. 8.

Scale invariance and power laws (complement to ch. 4)

Near a critical point, as in the Ising model at T_c, fluctuations have no characteristic size. The distributions of cluster or event sizes then follow power laws, for a simple reason.

An exponential law e^{-s/s_0} has the scale s_0; a power law has none. A straight line on a log-log plot is weak evidence: Clauset, Shalizi and Newman (2009) give maximum-likelihood estimators and goodness-of-fit tests.

Mean-field limits and propagation of chaos (complement to ch. 2, 3 and 8)

The mean-field equations of ch. 2 (contact process) and ch. 3 (Curie–Weiss) are instances of a general limit. Consider N exchangeable particles in \mathbb{R}^d, driven by independent Brownian motions W_i: dX_i = \frac1N \sum_{j=1}^N b(X_i, X_j)\, dt + \sigma\, dW_i .

Proof sketch. Couple each X_i with the nonlinear process \bar X_i, driven by the same W_i with drift \int b(\bar X_i, y)\, \mu_t(dy). The difference of drifts splits into a Lipschitz term in |X_i - \bar X_i| and a law-of-large-numbers term of order N^{-1/2}.

Gronwall’s lemma then gives \mathbb{E}\, |X_i - \bar X_i| \leq C_T N^{-1/2} on [0, T].

The pairwise kernels learned in ch. 8 are the drifts b(x, y) = \phi(\|y - x\|)(y - x) of this framework. Gradient descent on a wide two-layer neural network is an interacting particle system of the same type, with one particle per neuron (Mei, Montanari and Nguyen, 2018; Chizat and Bach, 2018).

Hydrodynamics of the exclusion process (complement to ch. 2)

In the totally asymmetric exclusion process (TASEP) on \mathbb{Z}, each particle jumps one step to the right at rate 1 if the site is empty. Under the invariant measure \nu_\rho of ch. 2, the current of particles across an edge is j(\rho) = \rho(1 - \rho).

At large scales, local equilibrium holds: around position x = \xi t, the configuration looks like \nu_{\rho} with a density \rho(t, x) that solves the Burgers equation \partial_t \rho + \partial_x \big( \rho (1 - \rho) \big) = 0 .

Why this profile. Density disturbances travel at the characteristic speed j'(\rho) = 1 - 2\rho. In a self-similar solution \rho(\xi), the density present at position \xi t is the one whose speed is \xi, so \xi = 1 - 2\rho. The front of the jam advances at speed 1 into empty space, and the jam melts backwards at speed 1.

Interacting particles, each moving at speed at most 1, produce a deterministic macroscopic flow with shocks and rarefaction fans: the prototype of traffic models.

Stochastic resetting and first passage (complement to ch. 10)

A diffusing searcher looking for a target can do better by restarting from scratch from time to time (Evans and Majumdar, 2011).

Let X be a Brownian motion with diffusion coefficient D, that is with generator D\, \partial_x^2, started at x_0 > 0.

It is reset to x_0 at the times of a Poisson process of rate r, and T is its first passage time at 0.

The simulation is exact: each excursion is a free first-passage time x_0^2 / (2 D Z^2), with Z standard Gaussian, unless the next reset comes first.

import numpy as np
from scipy.optimize import brentq

rng = np.random.default_rng(15)
D, x0 = 1.0, 1.0
z_star = brentq(lambda z: z / 2 - 1 + np.exp(-z), 0.5, 3)        # optimality condition
r_star = D * z_star**2 / x0**2
mfpt = lambda r: (np.exp(x0 * np.sqrt(r / D)) - 1) / r

def first_passage(r):
    """Exact simulation: free first-passage times x0^2 / (2 D Z^2) cut by Exp(r) resets."""
    total = 0.0
    while True:
        T = x0**2 / (2 * D * rng.standard_normal() ** 2)          # free first-passage time
        R = rng.exponential(1 / r)                                 # next reset
        if T < R:
            return total + T
        total += R

for r in [0.5, r_star, 8.0]:
    sim = np.mean([first_passage(r) for _ in range(40000)])
    print(f"r = {r:5.3f}: simulated MFPT {sim:.3f}, formula (e^z - 1)/r = {mfpt(r):.3f}")
print(f"optimal z* = {z_star:.4f}, r* = {r_star:.4f}, minimal MFPT = {mfpt(r_star):.4f} x0^2/D")

The simulated means are 2.052, 1.542 and 2.018 against 2.056, 1.544 and 1.990.

In first-passage resetting, the reset is triggered by the process itself, when it reaches a threshold, rather than by a clock. De Bruyne, Randon-Furling and Redner (2020) study optimisation problems for such processes. The simplest case is already explicit.

The check below uses a lattice walk with steps \pm \delta every \delta^2 / (2D), for which gambler’s ruin and the mean exit time are exactly those of Brownian motion.

import numpy as np

rng = np.random.default_rng(20)
D, x0, L, n = 1.0, 1.0, 2.0, 4000
delta = 0.05                                     # lattice walk: a step of +-delta every delta^2 / (2D)
m, M = round(x0 / delta), round(L / delta)
pos, steps, resets = np.full(n, m), np.zeros(n), np.zeros(n)
alive = np.ones(n, dtype=bool)
while alive.any():
    pos[alive] += rng.choice([-1, 1], size=alive.sum())
    steps[alive] += 1
    hit = alive & (pos == M)
    pos[hit], resets[hit] = m, resets[hit] + 1   # reset to x0 at the threshold L
    alive &= pos > 0                             # absorbed at the target 0
T = steps * delta**2 / (2 * D)
print(f"E T: simulated {T.mean():.3f} (s.e. {T.std() / np.sqrt(n):.3f}), formula x0 L / (2D) = {x0 * L / (2 * D):.3f}")
print(f"E R: simulated {resets.mean():.3f} (s.e. {resets.std() / np.sqrt(n):.3f}), formula x0 / (L - x0) = {x0 / (L - x0):.3f}")

With x_0 = D = 1 and L = 2, the simulation gives 1.002 \pm 0.015 and 1.006 \pm 0.022 against 1 and 1. Evans, Majumdar and Schehr (2020) review the field.

Rare-event simulation (complement to ch. 9)

Crude Monte Carlo needs about 100/p samples to estimate a probability p with 10% relative error. The two principles of ch. 9 lead to two different remedies.

• Light tails: exponential tilting. Sample the X_i from the tilted law of Theorem 9.1, under which the rare event is typical, and reweight by the likelihood ratio e^{-\theta S_n + n \Lambda(\theta)}. With \Lambda'(\theta) = a, the relative error grows only polynomially in n.
• Heavy tails: conditioning on the big jump. Tilting fails, since \Lambda = \infty for \theta > 0. The Asmussen–Kroese estimator uses the catastrophe principle instead: with M_{n-1} and S_{n-1} the maximum and sum of n - 1 summands, \mathbb{P}(S_n > x) = n\, \mathbb{E}\, \bar F\big( \max(M_{n-1},\ x - S_{n-1}) \big).

import numpy as np
from scipy.stats import norm

rng = np.random.default_rng(16)
R = 20000                                         # samples per estimator

# Light tails: P(S_n > n a) for Gaussian steps, by crude Monte Carlo and by exponential tilting
n, a = 50, 0.5
exact = norm.sf(n * a / np.sqrt(n))
crude = (rng.standard_normal((R, n)).sum(1) > n * a)
Xt = rng.standard_normal((R, n)) + a              # tilted law: N(a, 1), theta = a
St = Xt.sum(1)
w = np.exp(-a * St + n * a**2 / 2) * (St > n * a)  # likelihood ratio dP / dP_theta on the event
print(f"Gaussian, exact {exact:.3e}: crude {crude.mean():.3e} (rel. s.e. {crude.std() / crude.mean() / np.sqrt(R):.2f}), "
      f"tilted {w.mean():.3e} (rel. s.e. {w.std() / w.mean() / np.sqrt(R):.3f})")

# Heavy tails: P(S_n > x) for Pareto(1.5) steps, by crude MC and by Asmussen-Kroese conditioning
alpha, n, x = 1.5, 10, 1000.0
tail = lambda y: np.where(y > 1, np.maximum(y, 1) ** -alpha, 1.0)
X = rng.pareto(alpha, (R, n)) + 1
crude = X.sum(1) > x
Y = X[:, :-1]
ak = n * tail(np.maximum(Y.max(1), x - Y.sum(1)))  # n * P(X_n > max(M_{n-1}, x - S_{n-1}))
print(f"Pareto, n P(X > x) = {n * x**-alpha:.3e}: crude {crude.mean():.3e} (rel. s.e. "
      f"{crude.std() / crude.mean() / np.sqrt(R):.2f}), conditional {ak.mean():.3e} "
      f"(rel. s.e. {ak.std() / ak.mean() / np.sqrt(R):.1e})")

With 20000 samples, tilting reduces the relative error from 50% to 1.4% on a probability of 2.0 \times 10^{-4}. The conditional estimator gives 3.295 \times 10^{-4} with a relative error below 10^{-3}, where crude Monte Carlo has 38%. The single-big-jump approximation n \bar F(x) = 3.16 \times 10^{-4} is 4% too low at this level.

Universality of persistence (complement to ch. 10)

The combinatorial lemmas of ch. 10 have a striking consequence for the probability that a walk stays positive.

Proof sketch. Sparre Andersen’s combinatorial lemma, a relative of Lemma 10.1, gives the generating-function identity \sum_{n \geq 0} \mathbb{P}(S_1 > 0, \dots, S_n > 0)\, s^n = \exp\Big( \sum_{k \geq 1} \frac{s^k}{k}\, \mathbb{P}(S_k > 0) \Big). By symmetry and continuity \mathbb{P}(S_k > 0) = 1/2, so the right-hand side is (1 - s)^{-1/2}, whose coefficients are \binom{2n}{n} 4^{-n}.

import numpy as np
from math import comb

rng = np.random.default_rng(17)
n, R = 10, 200_000
steps = {"Gaussian": rng.standard_normal((R, n)),
         "Cauchy": rng.standard_cauchy((R, n)),
         "uniform": rng.uniform(-1, 1, (R, n))}
for name, X in steps.items():
    stays_positive = (np.cumsum(X, axis=1) > 0).all(axis=1)
    print(f"{name:8s}: P(S_1 > 0, ..., S_10 > 0) = {stays_positive.mean():.4f}")
print(f"Sparre Andersen: binom(20, 10) / 4^10 = {comb(2 * n, n) / 4**n:.4f}")

The three estimates, 0.1769, 0.1764 and 0.1768, agree with 0.1762 within Monte Carlo error. Majumdar (2010) surveys these universal first-passage properties and their use in physics.

13. Research directions and open problems

This chapter is a self-study tour of active questions that extend the course. It cites only work whose content it describes; the references give the entry points.

Inference for agent-based models

• Simulation-based inference. The likelihood of an ABM is almost never computable, but the model can be simulated. Approximate Bayesian computation (Beaumont, Zhang and Balding, 2002) accepts parameters whose simulations resemble the data. Neural methods learn posteriors or likelihood ratios from simulations (Cranmer, Brehmer and Louppe, 2020). Open questions: which parameters of an ABM are identifiable from macro-data, and how to detect misspecification.
• Surrogates. When one run is expensive, a statistical emulator of the simulator is fitted on a design of runs and used for calibration and sensitivity analysis; see Learning Physical Phenomena. Building surrogates that keep the rare collective events of ch. 6 and ch. 9 is open.
• Learning the rules. Lu, Maggioni and Tang (2021) extend the learning of interaction kernels to heterogeneous systems of agents observed along multiple trajectories. Sparse regression on a library of candidate terms identifies governing equations from data (Brunton, Proctor and Kutz, 2016). Learning rules of stochastic agents, on unknown networks, from partial observations, is largely open; Krakauer (2023) frames the general programme of unifying complexity science and machine learning.

Interaction structure and causality

• Sample-optimal inverse Ising. Vuffray, Misra, Lokhov and Chertkov (2016) introduce the interaction screening estimator and show that it recovers Ising graphs with a near-optimal number of samples; Nguyen, Zecchina and Berg (2017) review inverse statistical problems from physics to data science. Learning with hidden units, or from dynamics instead of equilibrium samples, raises new identifiability questions.
• Interference. Aronow and Samii (2017) estimate average causal effects under general interference, extending Hudgens and Halloran (2008). How to design experiments on networks with spillovers, and how to validate the designs on ABMs with a known truth (Manzo, 2022), are active topics.

Heavy tails, learning and optimisation

• Heavy tails in deep learning. Simsekli, Sagun and Gurbuzbalaban (2019) report heavy-tailed gradient noise in deep networks. Hodgkinson and Mahoney (2021) show that multiplicative noise produces heavy tails in stochastic optimisation, the mechanism of ch. 11. Simsekli, Sener, Deligiannidis and Erdogdu (2020) relate heavy tails to the Hausdorff dimension of the trajectories and to generalisation bounds. Whether heavy tails help generalisation in practice is debated.
• Decentralised optimisation. Gurbuzbalaban and coauthors (2022) study the tail index of decentralised SGD; how the communication graph should be chosen when tails matter is open.
• Training as an interacting system. The mean-field analyses of Mei, Montanari and Nguyen (2018) and Chizat and Bach (2018) treat neurons as interacting particles (Theorem 12.4). Quantitative propagation of chaos along training, and the corrections at finite width, remain active questions.

Extrema, criticality and resetting

• Convex hulls beyond the mean. Majumdar, Comtet and Randon-Furling (2010) review the mean perimeter and area of hulls of random walks, Brownian paths and Lévy processes. Kabluchko, Vysotsky and Zaporozhets (2017) prove distribution-free formulas for the faces of convex hulls of random walks in any dimension. Full distributions of the perimeter and area, and hulls of interacting walkers, are open.
• Timeliness criticality. Moran and coauthors (2024) find a critical buffer in networks of agents with delays, with avalanches near it. The universality class of delay avalanches on structured networks, finite-size corrections like those of ch. 6, and the control of buffers in real systems are open.
• Resetting and restarts. Evans, Majumdar and Schehr (2020) review stochastic resetting; de Bruyne, Randon-Furling and Redner (2020) study optimisation in first-passage resetting. In computer science, Luby, Sinclair and Zuckerman (1993) give optimal restart strategies for randomised algorithms. Resetting in systems of interacting searchers is a natural next step.

Exercises

The exercises are sorted by level: ★ applies the course, ★★ combines results or proves a short statement, ★★★ is open-ended. The session plan says which ones to do after each session.

1. ★ The two-state chain. A site switches from 0 to 1 at rate a and from 1 to 0 at rate b. Write the rate matrix Q, compute P_t = e^{tQ}, and find the stationary law. Check detailed balance.

   Solution: Q has rows (-a, a) and (b, -b), with eigenvalues 0 and -(a+b). Hence P_t(0, 1) = \frac{a}{a + b} \big(1 - e^{-(a+b)t}\big), \qquad P_t(1, 0) = \frac{b}{a + b} \big(1 - e^{-(a+b)t}\big). The stationary law is \pi = (b, a)/(a + b), and \pi(0)\, a = \pi(1)\, b = ab/(a+b) is detailed balance.

2. ★ Extinction. A Galton–Watson process has offspring law \mathbb{P}(\xi = 0) = \mathbb{P}(\xi = 1) = 1/4, \mathbb{P}(\xi = 2) = 1/2. Compute its mean and its extinction probability from one and from three ancestors.

   Solution: m = 1/4 + 1 = 1.25 > 1. The equation G(s) = \frac14 + \frac{s}{4} + \frac{s^2}{2} = s reads 2s^2 - 3s + 1 = 0, with roots 1/2 and 1. Hence q = 1/2 from one ancestor and q^3 = 1/8 from three (independent subtrees).

3. ★ Independent flips. Each site of a finite set S flips from 0 to 1 at rate a and from 1 to 0 at rate b, independently. Write the generator and show that the product measure \nu = \mathrm{Bernoulli}(a/(a+b))^{\otimes S} is reversible.

   Solution: c(x, \eta) = a (1 - \eta(x)) + b\, \eta(x) and L f(\eta) = \sum_x c(x, \eta) [f(\eta^x) - f(\eta)]. For \eta(x) = 0, \nu(\eta^x)/\nu(\eta) = a/b, which equals c(x, \eta)/c(x, \eta^x) = a/b. This is detailed balance for every flip.

4. ★ The voter model on a star. A star has a centre and n leaves. Compute the probability of consensus on 1 when only the centre holds opinion 1, and when only one leaf does.

   Solution: the degrees are n for the centre and 1 for each leaf, with 2|E| = 2n. By Proposition 2.4 the probabilities are n/(2n) = 1/2 and 1/(2n), far from the unweighted fractions 1/(n+1) and 1/(n+1).

5. ★ Metropolis or heat-bath? For \beta \Delta = 2, compute the Metropolis and heat-bath rates of a flip and of its reverse. Show that the Metropolis rate is always at least the heat-bath rate.

   Solution: Metropolis gives 0.135 and 1; heat-bath gives 1/(1 + e^2) = 0.119 and 0.881. Both ratios equal e^{-2}.

   For u \geq 0, e^{-u} \geq 1/(1 + e^u) because 1 + e^u \geq e^u; for u < 0, the Metropolis rate is 1. By Peskun’s (1973) ordering, Metropolis therefore has an asymptotic variance no larger than heat-bath for every observable.

6. ★ The Borel law. For critical Poisson branching (\mu = 1), compute \mathbb{P}(S = s) for s = 1, 2, 3, and check the first value directly. Compute \mathbb{E}\, S for \mu = 1/2.

   Solution: e^{-1} = 0.368, e^{-2} = 0.135 and \frac32 e^{-3} = 0.0747. Directly, S = 1 means no child, of probability e^{-1}.

   For \mu = 1/2, \mathbb{E}\, S = 1/(1 - \mu) = 2.

7. ★ Final size and vaccination. For R_0 = 2, compute the final size, the herd-immunity threshold and the overshoot. Then vaccinate 30% of the population beforehand and compute the fraction infected.

   Solution: z = 1 - e^{-2z} gives z = 0.797; the threshold is 0.5, so the overshoot is 0.297.

   With 30% immune, the susceptibles are 70% and the effective reproduction number is 1.4. The final fraction y of the susceptibles solves y = 1 - e^{-1.4 y}, so y = 0.511, and 0.7\, y = 0.358 of the population is infected.

8. ★ Kirman’s model at a = 1 (Python). Simulate Kirman’s model with N = 50 for a = 0.5, 1 and 2, and compare the time fractions with the beta-binomial law. Show that the stationary law is uniform on \{0, \dots, N\} when a = 1.

   Solution: with a = 1, \binom{N}{k} B(k + 1, N - k + 1) = 1/(N + 1), so \pi(k) = 1/(N + 1): the transition between unimodal and bimodal laws is a flat law. The simulation below gives total variations of a few percent:

   import numpy as np
   from scipy.stats import betabinom
   
   rng = np.random.default_rng(18)
   N, h = 50, 1.0
   for a in [0.5, 1.0, 2.0]:
       eps = a * h / (N - 1)
       k, time_in = N // 2, np.zeros(N + 1)
       for _ in range(300000):
           up = (N - k) * (eps + h * k / (N - 1))
           down = k * (eps + h * (N - k) / (N - 1))
           time_in[k] += rng.exponential(1 / (up + down))
           k += 1 if rng.random() < up / (up + down) else -1
       emp = time_in / time_in.sum()
       exact = betabinom.pmf(np.arange(N + 1), N, a, a)
       print(f"a = {a}: total variation to beta-binomial {0.5 * np.abs(emp - exact).sum():.3f}, "
             f"max |pi(k) - 1/(N+1)| = {np.abs(exact - 1 / (N + 1)).max():.4f}")

   The output gives total variations 0.033, 0.023 and 0.037, and a deviation from uniform of 0.0000 for a = 1.

9. ★★ Mean field for the contact process. On the complete graph with N sites and infection rate \lambda/N per infected–healthy pair, compute L f for f(\eta) = |\eta|/N. Derive the mean-field equation and study the stability of its fixed points.

   Solution: with \rho = |\eta|/N, L f = \frac{\lambda}{N} \cdot \frac{|\eta| (N - |\eta|)}{N} - \frac{|\eta|}{N} = \lambda \rho (1 - \rho) - \rho . Neglecting fluctuations, which are of order N^{-1/2}, gives \dot\rho = \lambda\rho(1 - \rho) - \rho.

   The derivative of the right-hand side is \lambda - 1 at 0, and 1 - \lambda at \rho^\ast = 1 - 1/\lambda. The endemic state exists and is stable exactly when \lambda > 1, and 0 is then unstable.

10. ★★ Curie–Weiss dynamics lumps. Show that heat-bath dynamics on the complete graph lumps to the number k of + spins. Write its birth and death rates, and check detailed balance for \pi(k) \propto \binom{N}{k} e^{\beta J M^2 / (2N)} with M = 2k - N.

    Solution: the dynamics is invariant under permutations of the sites, so Proposition 7.3 applies. Let u = 2\beta J (M + 1)/N.

    A - spin sees the field (J/N)(M + 1), so k \to k + 1 at rate (N - k)/(1 + e^{-u}).

    In the state k + 1, a + spin sees the same field, so k + 1 \to k at rate (k + 1)/(1 + e^{u}).

    The ratio of the rates is \frac{N - k}{k + 1} e^{u}, and so is \pi(k + 1)/\pi(k), since (M + 2)^2 - M^2 = 4(M + 1).

11. ★★ Critical buffers (Python). Compute v_K of Theorem 6.3 for K = 2, 3, 5 and compare it with the mean of the largest of K independent \mathrm{Exp}(1) delays, H_K = 1 + \frac12 + \dots + \frac1K. Explain the gap.

    Hints and expected result: minimising (\ln K - \ln(1 - \theta))/\theta numerically gives 2.678, 3.289 and 3.994 for K = 2, 3 and 5, against H_K = 1.5, 1.833 and 2.283.

    A buffer equal to the mean largest input delay is far from enough: delays accumulate along the paths of the dependency tree, and the critical buffer is set by the best of K^s paths, a large-deviation quantity.

12. ★★ Lumping the path. For the voter model on the path 1 - 2 - 3, use the reflection 1 \leftrightarrow 3 to build a lumpable partition, and compute the probability of consensus on 1 from (0, 1, 0) and from (1, 0, 0).

    Solution: the orbits are {000}, {111}, {100, 001}, {010}, {110, 011} and {101}, lumpable by Proposition 7.3.

    With degrees (1, 2, 1) and 2|E| = 4, Proposition 2.4 gives 2/4 = 1/2 from 010 and 1/4 from 100.

13. ★★ Chernoff and its prefactor. For i.i.d. \mathrm{Exp}(1) variables, compute I(a) and compare the Chernoff bound with the exact value of \mathbb{P}(S_{20} \geq 40). Explain the remaining factor.

    Solution: \Lambda(\theta) = -\ln(1 - \theta), so I(a) = a - 1 - \ln a and I(2) = 0.307. The bound is e^{-20 \times 0.307} = 2.16 \times 10^{-3}, while S_{20} is Gamma(20, 1) with \mathbb{P}(S_{20} \geq 40) = 1.76 \times 10^{-4}.

    The Bahadur–Rao refinement multiplies the bound by 1/(\theta \sigma \sqrt{2\pi n}), with \theta = 1/2 and tilted standard deviation \sigma = a = 2. This factor is 0.089, giving 1.93 \times 10^{-4}.

14. ★★ Kac’s formula in two steps (Python). For standard Gaussian steps, compute \mathbb{E} \max(0, S_1, S_2) with Theorem 10.2 and check it by simulation.

    Solution: \mathbb{E}\, S_1^+ = 1/\sqrt{2\pi} and \mathbb{E}\, S_2^+ / 2 = \sqrt 2 / (2\sqrt{2\pi}), so \mathbb{E} \max(0, S_1, S_2) = \frac{1 + 1/\sqrt 2}{\sqrt{2\pi}} = 0.6810 .

    import numpy as np
    
    rng = np.random.default_rng(19)
    S = np.cumsum(rng.standard_normal((1_000_000, 2)), axis=1)
    print(f"simulated E max(0, S_1, S_2) = {np.maximum(0, S.max(1)).mean():.4f}, "
          f"Kac {(1 + 1 / np.sqrt(2)) / np.sqrt(2 * np.pi):.4f}")

    The simulation gives 0.6816.

15. ★★ When is the variance of SGD infinite? In the model of ch. 11 with batch size b, find the step size at which \alpha = 2. Evaluate it for b = 1 and b = 10.

    Solution: \chi^2_b / b has mean 1 and variance 2/b, so \mathbb{E}\Big(1 - \eta\, \frac{\chi^2_b}{b}\Big)^2 = 1 - 2\eta + \eta^2 \Big(1 + \frac2b\Big). It equals 1 at \eta = 2b/(b + 2): 2/3 for b = 1 and 5/3 for b = 10.

    From this step size on, the convexity of p \mapsto \mathbb{E}|A|^p gives \alpha \leq 2, and Theorem 11.1 makes the stationary variance infinite.

This holds as long as a stationary law exists, that is while \mathbb{E} \ln|A| < 0: numerically up to \eta \approx 2.42 for b = 1 and 2.43 for b = 10.

16. ★★ Optimal resetting. Prove the optimality condition of Proposition 12.6 and the value of the minimum.

    Solution: \mathbb{E}\, T = (x_0^2/D)\, g(z) with g(z) = (e^z - 1)/z^2. Then g'(z) = 0 gives z e^z = 2(e^z - 1), that is z/2 = 1 - e^{-z}.

    Its positive root is z^\ast = 1.5936, and g(z^\ast) = 3.9215/2.5396 = 1.544.

17. ★★★ Interference in an epidemic ABM. Build an SIR model on a population split into 50 communities, with most contacts inside communities. Randomise each community to a vaccine coverage of 30% or 70%, then individuals within communities. Estimate the direct, indirect and total effects of Hudgens and Halloran (2008), and compare them with a naive individual-level comparison.

    Hints and expected result: the indirect effect, the drop in risk of unvaccinated people between low- and high-coverage communities, is large when R_0 (1 - \text{coverage}) crosses 1 between the two arms. The naive vaccinated-versus-unvaccinated difference mixes the direct effect with the coverage of each community and changes with the design. Report confidence intervals over replications.

18. ★★★ The learning rate of interaction kernels. Repeat the experiment of ch. 8 with the smooth kernel \phi(r) = e^{-r^2}, for M = 5, 10, 20, 50, 100 trajectories, with a number of bins proportional to M^{1/3}. Estimate the empirical rate of the L^2(\rho_T) error and compare it with Theorem 8.4.

    Hints and expected result: for a Lipschitz kernel (s = 1) the theorem predicts an error of order (\ln M / M)^{1/3} when bias and variance are balanced. A log-log fit should give a slope near -1/3. With noise-free velocities the variance vanishes and only the bias, set by the number of bins, remains, which shows why the noise level matters.

19. ★★★ Avalanches of delays. In the model of Theorem 6.3 with N = 5000 and K = 2, choose B slightly above the finite-N threshold. Call an avalanche a maximal run of steps during which the mean delay exceeds a fixed level, and record its size, the summed excess delay. Study the size distribution and how it changes with B and N.

    Hints and expected result: near the threshold, sizes spread over several decades, with a cutoff that grows as B approaches the threshold. Compare the tail with the s^{-3/2} law of critical branching (Theorem 4.3) and discuss why the delay recursion, a maximum over paths, need not share that exponent. Moran and coauthors (2024) study this regime.

Projects and assessment

The official assessment is “Project based” (master-mva.com); the Université Paris-Saclay course detail classifies it as “Evaluation Continue non Intégrale” (partial continuous assessment), and no weighting is displayed. The subjects below are sized for such a project: most are paper studies, in the MVA tradition, with a numerical success criterion. The written exam is only a proposal, for self-assessment or for a written component.

Solo projects

1. [Solo · ≈ 20 h] Learning interaction laws: a paper study. Reproduce the opinion-dynamics experiments of Lu, Zhong, Tang and Maggioni (2019) on synthetic trajectories.
   • Implement: the least-squares estimator of ch. 8 on piecewise-constant and piecewise-linear spaces, with the number of basis functions growing with M; an empirical check of the coercivity constant.
   • Success criterion: a log-log slope of the error within 0.1 of -s/(2s+1) for a smooth kernel, and trajectories predicted by the learned kernel that stay within 5% of the true ones up to the observation horizon.
   • Deliverables: code, a 4–6 page report, a 10-minute talk. Chapters 8, 12.
2. [Solo · ≈ 20 h] Convex hulls of random walks and Lévy flights: a paper study. Study Majumdar, Comtet and Randon-Furling (2010) and Randon-Furling, Majumdar and Comtet (2009).
   • Implement: hulls of planar walks with Gaussian, uniform and \alpha-stable steps; the formulas of Theorems 10.3–10.4; the hull of n independent Brownian paths.
   • Success criterion: mean perimeters and areas within three standard errors of the formulas for n = 10^3; the growth exponent 1/\alpha of Proposition 10.5 measured within 0.05 for \alpha = 1.5.
   • Deliverables: code, a 4–6 page report, a 10-minute talk. Chapters 10, 12.
3. [Solo · ≈ 20 h] Heavy tails in SGD: a paper study. Reproduce the dependence of the tail index on the step size and the batch size found by Gurbuzbalaban, Simsekli and Zhu (2021) for linear regression with Gaussian data in dimension 10.
   • Implement: SGD, the estimation of the Kesten exponent from products of random matrices, the Hill estimator with a stability plot over k.
   • Success criterion: both monotone trends reproduced; the index estimated within 20% where the theory computes it exactly (d = 1, the equation of ch. 11).
   • Deliverables: code, a 4–6 page report, a 10-minute talk. Chapters 9, 11.

Group projects

1. [Group of 2–3 · ≈ 30 h each] An epidemic ABM as a policy laboratory. Build a synthetic population with households, schools and workplaces, and use an SIR ABM to compare interventions, in the spirit of the models described by Adam (2020).
   • Work split: one member writes the exact Gillespie simulator and checks it against Propositions 5.3–5.4 in the well-mixed limit; one calibrates R_0 from a synthetic outbreak by approximate Bayesian computation; one designs the intervention and interference experiments (Exercise 17).
   • Success criterion: outbreak probability and final size within Monte Carlo error of the theory in the well-mixed limit; every scenario reported with intervals over at least 100 runs and common random numbers.
   • Deliverables: code, a 6–8 page report, a 15-minute talk. Chapters 0, 5, 8.
2. [Group of 2–3 · ≈ 30 h each] Schelling, Kawasaki and aggregation. Study the logit Schelling dynamics of Proposition 5.2 as a lattice gas, and reduce small instances exactly (Banisch, 2015).
   • Work split: one member maps the phase diagram in \beta and the vacancy rate; one computes exact stationary laws and hitting times on small grids by lumping with the symmetries of the torus; one compares mean-field and pair approximations with both.
   • Success criterion: on a 3 \times 3 torus, the long-run frequencies within total variation 0.01 of the exact Gibbs law; the segregation index of the large-grid simulations explained by the phase diagram.
   • Deliverables: code, a 6–8 page report, a 15-minute talk. Chapters 3, 5, 7.
3. [Group of 3 · ≈ 30 h each] Timeliness criticality: a paper study. Reproduce the phase transition of delay propagation studied by Moran and coauthors (2024).
   • Work split: one member runs large simulations and the finite-size analysis; one develops the branching-random-walk theory of Theorem 6.3 and its finite-N corrections (Brunet and Derrida, 1997); one studies structured dependency networks (lattices, random regular graphs) and avalanche statistics.
   • Success criterion: v_K of Theorem 6.3 recovered to two decimals by extrapolation in N for K = 2, 3; the avalanche-size distribution near the threshold over at least two decades.
   • Deliverables: code, a 6–8 page report, a 15-minute talk. Chapters 4, 6, 9.

Written exam (proposal)

Format: 2 hours, one handwritten A4 sheet allowed, no electronic device; four independent exercises, 20 points.

1. Course question (4 points). Define a stochastic Ising model. Show that the heat-bath rate satisfies the definition, and that Kawasaki dynamics conserves the magnetisation.

   Expected: Definition 3.2; the ratio of the heat-bath rates of a flip and its reverse is e^{-\beta\Delta}, using \Delta_x(\sigma^x) = -\Delta_x(\sigma); an exchange of two unequal neighbouring spins keeps the number of + spins, and the dynamics is reversible for the measure conditioned on the magnetisation.

2. Computation (6 points). A stochastic SIR epidemic has R_0 = 3. (a) Give the probability of a major outbreak from one and from two initial cases. (b) Give the final-size equation, its root, and the herd-immunity threshold. (c) What fraction must be vaccinated beforehand to prevent major outbreaks?

   Expected: (a) 1 - 1/3 = 2/3 and 1 - 1/9 = 8/9; (b) z = 1 - e^{-3z}, z = 0.940, threshold 2/3; (c) more than 1 - 1/R_0 = 2/3, so that the effective reproduction number falls below 1. Two points per part.

3. Short proof (5 points). Let X_1, X_2 be i.i.d. Pareto with index \alpha. Prove that \mathbb{P}(X_1 + X_2 > x) \sim 2\, \mathbb{P}(X_1 > x). For i.i.d. \mathrm{Exp}(1) variables, compute the same ratio.

   Expected: the lower and upper bounds of Theorem 9.4 with the split at \delta x and (1 - \delta)x (3 points). For exponentials, X_1 + X_2 is Gamma(2, 1), with tail (1 + x)e^{-x}, so the ratio is 1 + x \to \infty: no single big jump (2 points).

4. Problem (5 points). For the voter model on a finite connected graph, prove that the degree-weighted magnetisation is a martingale. On the path with 4 vertices, compute the probability of consensus on 1 when the two middle vertices hold opinion 1.

   Expected: the generator computation of Proposition 2.4 and optional stopping (3 points). Degrees (1, 2, 2, 1), total 6, so the probability is 4/6 = 2/3 (2 points).

Grading rubric (proposal)

References

Official references (master-mva.com)

• D. Adam. The simulations driving the world’s response to COVID-19. Nature, 580(7802), 316–319, 2020.
• D. Aldous. Interacting particle systems as stochastic social dynamics. Bernoulli, 19(4), 1122–1149, 2013.
• S. Banisch. Markov Chain Aggregation for Agent-Based Models. Springer, 2015.
• B. de Bruyne, J. Randon-Furling and S. Redner. Optimization in first-passage resetting. Physical Review Letters, 125(5), 050602, 2020.
• J. M. Epstein. Generative Social Science: Studies in Agent-Based Computational Modeling. Princeton University Press, 2012.
• M. Gurbuzbalaban, Y. Hu, U. Simsekli, K. Yuan and L. Zhu. Heavy-tail phenomenon in decentralized SGD. Preprint arXiv:2205.06689, 2022.
• D. C. Krakauer. Unifying complexity science and machine learning. Frontiers in Complex Systems, 1, 1235202, 2023.
• T. M. Liggett. Interacting Particle Systems. Springer, 2005 (reprint of the 1985 edition).
• F. Lu, M. Zhong, S. Tang and M. Maggioni. Nonparametric inference of interaction laws in systems of agents from trajectory data. Proceedings of the National Academy of Sciences, 116(29), 14424–14433, 2019.
• S. N. Majumdar, A. Comtet and J. Randon-Furling. Random convex hulls and extreme value statistics. Journal of Statistical Physics, 138, 955–1009, 2010.
• G. Manzo. Agent-Based Models and Causal Inference. John Wiley & Sons, 2022.
• J. Moran, M. Romeijnders, P. Le Doussal, F. P. Pijpers, U. Weitzel, D. Panja and J.-P. Bouchaud. Timeliness criticality. Preprint arXiv:2309.15070, 2024 (announced to appear in Nature Physics).
• J. P. Sethna. Statistical Mechanics: Entropy, Order Parameters, and Complexity, 2nd edition. Oxford University Press, 2021.
• Cited in the official session list as “Nair 2022”: J. Nair, A. Wierman and B. Zwart. The Fundamentals of Heavy Tails: Properties, Emergence, and Estimation. Cambridge University Press, 2022.

Textbooks

• A. Dembo and O. Zeitouni. Large Deviations Techniques and Applications, 2nd edition. Springer, 1998.
• J. G. Kemeny and J. L. Snell. Finite Markov Chains. Van Nostrand, 1960.
• D. A. Levin, Y. Peres and E. L. Wilmer. Markov Chains and Mixing Times. American Mathematical Society, 2009.
• M. Mézard and A. Montanari. Information, Physics, and Computation. Oxford University Press, 2009.
• J. R. Norris. Markov Chains. Cambridge University Press, 1997.
• J. Pearl. Causality: Models, Reasoning, and Inference, 2nd edition. Cambridge University Press, 2009.

Particle systems and stochastic Ising models (ch. 0–3, 12)

• M. Dwass. The total progeny in a branching process and a related random walk. Journal of Applied Probability, 6(3), 682–686, 1969.
• D. T. Gillespie. Exact stochastic simulation of coupled chemical reactions. Journal of Physical Chemistry, 81(25), 2340–2361, 1977.
• R. J. Glauber. Time-dependent statistics of the Ising model. Journal of Mathematical Physics, 4(2), 294–307, 1963.
• T. E. Harris. Nearest-neighbor Markov interaction processes on multidimensional lattices. Advances in Mathematics, 9(1), 66–89, 1972.
• R. A. Holley and T. M. Liggett. Ergodic theorems for weakly interacting infinite systems and the voter model. Annals of Probability, 3(4), 643–663, 1975.
• K. Kawasaki. Diffusion constants near the critical point for time-dependent Ising models. I. Physical Review, 145(1), 224–230, 1966.
• N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller and E. Teller. Equation of state calculations by fast computing machines. Journal of Chemical Physics, 21(6), 1087–1092, 1953.
• L. Onsager. Crystal statistics. I. A two-dimensional model with an order-disorder transition. Physical Review, 65(3–4), 117–149, 1944.
• P. H. Peskun. Optimum Monte-Carlo sampling using Markov chains. Biometrika, 60(3), 607–612, 1973.
• H. Rost. Non-equilibrium behaviour of a many particle process: density profile and local equilibria. Zeitschrift für Wahrscheinlichkeitstheorie und verwandte Gebiete, 58(1), 41–53, 1981.
• A.-S. Sznitman. Topics in propagation of chaos. In École d’Été de Probabilités de Saint-Flour XIX — 1989, Lecture Notes in Mathematics 1464, Springer, 1991.
• C. N. Yang. The spontaneous magnetization of a two-dimensional Ising model. Physical Review, 85(5), 808–816, 1952.

Agent-based models, criticality and markets (ch. 4–7)

• P. Bak, C. Tang and K. Wiesenfeld. Self-organized criticality: an explanation of the 1/f noise. Physical Review Letters, 59(4), 381–384, 1987.
• J. D. Biggins. Chernoff’s theorem in the branching random walk. Journal of Applied Probability, 14(3), 630–636, 1977.
• L. E. Blume. The statistical mechanics of strategic interaction. Games and Economic Behavior, 5(3), 387–424, 1993.
• É. Brunet and B. Derrida. Shift in the velocity of a front due to a cutoff. Physical Review E, 56(3), 2597–2604, 1997.
• A. Clauset, C. R. Shalizi and M. E. J. Newman. Power-law distributions in empirical data. SIAM Review, 51(4), 661–703, 2009.
• R. Cont and J.-P. Bouchaud. Herd behavior and aggregate fluctuations in financial markets. Macroeconomic Dynamics, 4(2), 170–196, 2000.
• D. Dhar. Self-organized critical state of sandpile automaton models. Physical Review Letters, 64(14), 1613–1616, 1990.
• J. M. Epstein. Agent-based computational models and generative social science. Complexity, 4(5), 41–60, 1999.
• J. M. Hammersley. Postulates for subadditive processes. Annals of Probability, 2(4), 652–680, 1974.
• J. F. C. Kingman. The first birth problem for an age-dependent branching process. Annals of Probability, 3(5), 790–801, 1975.
• A. Kirman. Ants, rationality, and recruitment. Quarterly Journal of Economics, 108(1), 137–156, 1993.
• M. E. J. Newman. Spread of epidemic disease on networks. Physical Review E, 66(1), 016128, 2002.
• T. C. Schelling. Dynamic models of segregation. Journal of Mathematical Sociology, 1(2), 143–186, 1971.

Inference and causality (ch. 8, 13)

• P. M. Aronow and C. Samii. Estimating average causal effects under general interference, with application to a social network experiment. Annals of Applied Statistics, 11(4), 1912–1947, 2017.
• M. A. Beaumont, W. Zhang and D. J. Balding. Approximate Bayesian computation in population genetics. Genetics, 162(4), 2025–2035, 2002.
• S. L. Brunton, J. L. Proctor and J. N. Kutz. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proceedings of the National Academy of Sciences, 113(15), 3932–3937, 2016.
• K. Cranmer, J. Brehmer and G. Louppe. The frontier of simulation-based inference. Proceedings of the National Academy of Sciences, 117(48), 30055–30062, 2020.
• R. Hegselmann and U. Krause. Opinion dynamics and bounded confidence: models, analysis and simulation. Journal of Artificial Societies and Social Simulation, 5(3), 2002.
• M. G. Hudgens and M. E. Halloran. Toward causal inference with interference. Journal of the American Statistical Association, 103(482), 832–842, 2008.
• F. Lu, M. Maggioni and S. Tang. Learning interaction kernels in heterogeneous systems of agents from multiple trajectories. Journal of Machine Learning Research, 22, 2021.
• H. C. Nguyen, R. Zecchina and J. Berg. Inverse statistical problems: from the inverse Ising problem to data science. Advances in Physics, 66(3), 197–261, 2017.
• P. Ravikumar, M. J. Wainwright and J. D. Lafferty. High-dimensional Ising model selection using ℓ1-regularized logistic regression. Annals of Statistics, 38(3), 1287–1319, 2010.
• M. Vuffray, S. Misra, A. Y. Lokhov and M. Chertkov. Interaction screening: efficient and sample-optimal learning of Ising models. Advances in Neural Information Processing Systems 29 (NeurIPS), 2016.

Extremes, heavy tails and learning (ch. 9–13)

• S. Asmussen and D. P. Kroese. Improved algorithms for rare event simulation with heavy tails. Advances in Applied Probability, 38(2), 545–558, 2006.
• R. R. Bahadur and R. Ranga Rao. On deviations of the sample mean. Annals of Mathematical Statistics, 31(4), 1015–1027, 1960.
• O. Barndorff-Nielsen and G. Baxter. Combinatorial lemmas in higher dimensions. Transactions of the American Mathematical Society, 108(2), 313–325, 1963.
• L. Chizat and F. Bach. On the global convergence of gradient descent for over-parameterized models using optimal transport. Advances in Neural Information Processing Systems 31 (NeurIPS), 2018.
• M. R. Evans and S. N. Majumdar. Diffusion with stochastic resetting. Physical Review Letters, 106(16), 160601, 2011.
• M. R. Evans, S. N. Majumdar and G. Schehr. Stochastic resetting and applications. Journal of Physics A: Mathematical and Theoretical, 53(19), 193001, 2020.
• C. M. Goldie. Implicit renewal theory and tails of solutions of random equations. Annals of Applied Probability, 1(1), 126–166, 1991.
• M. Gurbuzbalaban, U. Simsekli and L. Zhu. The heavy-tail phenomenon in SGD. Proceedings of the 38th International Conference on Machine Learning (ICML), PMLR 139, 2021.
• B. M. Hill. A simple general approach to inference about the tail of a distribution. Annals of Statistics, 3(5), 1163–1174, 1975.
• L. Hodgkinson and M. W. Mahoney. Multiplicative noise and heavy tails in stochastic optimization. Proceedings of the 38th International Conference on Machine Learning (ICML), PMLR 139, 2021.
• Z. Kabluchko, V. Vysotsky and D. Zaporozhets. Convex hulls of random walks, hyperplane arrangements, and Weyl chambers. Geometric and Functional Analysis, 27(4), 880–918, 2017.
• H. Kesten. Random difference equations and renewal theory for products of random matrices. Acta Mathematica, 131, 207–248, 1973.
• M. Luby, A. Sinclair and D. Zuckerman. Optimal speedup of Las Vegas algorithms. Information Processing Letters, 47(4), 173–180, 1993.
• S. N. Majumdar. Universal first-passage properties of discrete-time random walks and Lévy flights on a line: statistics of the global maximum and records. Physica A, 389(20), 4299–4316, 2010.
• S. Mei, A. Montanari and P.-M. Nguyen. A mean field view of the landscape of two-layer neural networks. Proceedings of the National Academy of Sciences, 115(33), E7665–E7671, 2018.
• J. Randon-Furling, S. N. Majumdar and A. Comtet. Convex hull of N planar Brownian motions: exact results and an application to ecology. Physical Review Letters, 103(14), 140602, 2009.
• U. Simsekli, L. Sagun and M. Gurbuzbalaban. A tail-index analysis of stochastic gradient noise in deep neural networks. Proceedings of the 36th International Conference on Machine Learning (ICML), PMLR 97, 2019.
• U. Simsekli, O. Sener, G. Deligiannidis and M. A. Erdogdu. Hausdorff dimension, heavy tails, and generalization in neural networks. Advances in Neural Information Processing Systems 33 (NeurIPS), 2020.
• E. Sparre Andersen. On the fluctuations of sums of random variables II. Mathematica Scandinavica, 2, 195–223, 1954.
• F. Spitzer. A combinatorial lemma and its application to probability theory. Transactions of the American Mathematical Society, 82(2), 323–339, 1956.
• F. Spitzer and H. Widom. The circumference of a convex polygon. Proceedings of the American Mathematical Society, 12(3), 506–509, 1961.