Computational Statistics

Contents

• Introduction and overview
• 0. Review: Monte Carlo, convexity and martingales
• 1. The stochastic gradient algorithm
• 2. Bayesian inference and its computational challenges
• 3. Majorize–Minimize and the EM algorithm
• 4. Stochastic approximation
• 5. Simulation of random variables
• 6. Markov chain Monte Carlo: Metropolis–Hastings and Gibbs
• 7. Stochastic versions of EM
• 8. Adaptive MCMC and parallel tempering
• 9. Approximate Bayesian computation
• 10. Research directions and open problems
• Exercises
• Projects and assessment
• References

Introduction and overview

Audience: an M2 student with a solid M1 in probability and optimisation who must fit models with hidden variables, compute Bayesian estimators, or sample from a distribution known only up to a constant.

Empirical first: EM fitted to a two-component Gaussian mixture raises its likelihood at every iteration, yet still gains 13 log-likelihood units between iterations 20 and 200 (ch. 3). A random-walk Metropolis chain started in one mode of a bimodal target never visits the other mode in 20 000 steps (ch. 8). The course explains both behaviours, and then repairs them.

The course is about the algorithms that turn a statistical model into numbers: stochastic gradients, EM and its stochastic versions, stochastic approximation, simulation and Markov chain Monte Carlo, and likelihood-free inference. It shares tools with Convex Optimization and Applications in Machine Learning, Introduction to Probabilistic Graphical Models and Deep Generative Models and Stochastic Calculus in Machine Learning, and it prepares Bayesian Machine Learning and Biostatistics.

What you will be able to do

• Prove the convergence rates of the stochastic gradient algorithm, and diagnose its step-size and conditioning bottlenecks.
• Write a posterior distribution, derive Bayes estimators for standard losses, and identify which integral makes a Bayesian computation hard.
• Derive EM for a latent-variable model, prove that it increases the likelihood, and predict its speed from the fraction of missing information.
• Recognise a Majorize–Minimize algorithm and build a surrogate for a new objective.
• Prove the almost-sure convergence of a Robbins–Monro algorithm, and choose step sizes and averaging to reach the optimal variance.
• Simulate random variables by inversion, transformation and rejection, and design an importance-sampling estimator for a rare event.
• Implement Metropolis–Hastings and Gibbs samplers, prove their invariance, and assess them with autocorrelations, effective sample sizes and several chains.
• Implement SAEM and MCMC-SAEM, and explain why one simulation per iteration suffices.
• Tune a sampler by adaptive MCMC without destroying its ergodicity, and sample multimodal targets with parallel tempering.
• Perform approximate Bayesian computation, and choose its tolerance and summary statistics knowingly.

Official objectives and syllabus (2026–2027)

The course description on master-mva.com states that the course details statistical computational methods and Bayesian inference. It starts with the stochastic gradient algorithm, its theoretical properties and its numerical bottlenecks. After a brief introduction to Bayesian inference, it discusses the numerical challenges of computing Bayesian estimators, which lead to the EM algorithm, its stochastic versions and, more generally, Majorize–Minimize algorithms.

Stochastic approximation is then detailed, to understand the convergence conditions of these algorithms. The second part is devoted to the simulation of random variables: from the rejection method to the classical MCMC samplers (Metropolis–Hastings and Gibbs), then recent advances in MCMC: adaptive MCMC, adaptive parallel tempering and approximate Bayesian computation. The description announces 20 h of lectures and 20 h of training and lab sessions, and lists the course for the Master 2 MVA and Data Science.

The Université Paris-Saclay course detail gives the volume (20 h of lectures, 20 h of labs, 5 ECTS) and the nature of the assessment, but no programme items. The table below has one row per topic of the master-mva.com description and per announced element of the organisation.

Beyond the syllabus, the page adds importance sampling (ch. 5.4), MCMC diagnostics (ch. 6.4) and, in self-study, Langevin and Hamiltonian samplers, stochastic-gradient MCMC and neural simulation-based inference (ch. 10).

Where this course fits in the programme

Official prerequisites: not displayed. The material assumes M1 probability (conditional expectation, laws of large numbers, central limit theorem, finite Markov chains), M1 optimisation (gradients, convexity) and fluent NumPy.

• Preparation — the refreshers Probability and Statistics for MVA, Linear Algebra, Analysis and Optimization for MVA and Scientific Python, NumPy and PyTorch for MVA; ch. 0 restates what this course uses.
• Related courses in semester 1 — Convex Optimization and Applications in Machine Learning (deterministic counterpart of ch. 1), Introduction to Probabilistic Graphical Models and Deep Generative Models (EM and variational inference), Stochastic Calculus in Machine Learning (Langevin sampling, ch. 10), Introduction to Statistical Learning and, in the Santé track, Medical Image Analysis Based on Generative, Geometric and Biophysical Models.
• Courses that build on it in semester 2 — Bayesian Machine Learning, Inverse Problems and Imaging: Statistical and Stochastic Approaches, Biostatistics (mixed-effects models, ch. 7), Reinforcement Learning and Sequential Learning (stochastic approximation, ch. 4), Graphical Models: Discrete Inference and Learning (Gibbs sampling), Diffusion and Multiscale Models in Artificial Intelligence (Langevin dynamics, ch. 10).

Course outline

• 0. Review: Monte Carlo, convexity and martingales — the Monte Carlo error bar, Kullback–Leibler divergence, exponential families, smooth and strongly convex functions, martingale convergence.
• 1. The stochastic gradient algorithm — unbiased gradients, rates in the strongly convex and convex cases, and the numerical bottlenecks.
• 2. Bayesian inference and its computational challenges — posterior, Bayes estimators, large-sample behaviour, and the integrals that make everything hard.
• 3. Majorize–Minimize and the EM algorithm — surrogate optimisation, EM as its probabilistic instance, Gaussian mixtures, convergence and speed.
• 4. Stochastic approximation — Robbins–Monro, the Robbins–Siegmund lemma, almost-sure convergence, the ODE method, rates and averaging.
• 5. Simulation of random variables — inversion, transformations, rejection and importance sampling.
• 6. Markov chain Monte Carlo: Metropolis–Hastings and Gibbs — invariance, ergodicity, the two classical samplers and their diagnostics.
• 7. Stochastic versions of EM — MCEM, SEM, SAEM and MCMC-SAEM, with the Lyapunov argument behind SAEM.
• 8. Adaptive MCMC and parallel tempering — learning the proposal, when adaptation is safe, crossing energy barriers with a temperature ladder.
• 9. Approximate Bayesian computation — inference by simulation, tolerance and summary statistics, ABC-MCMC and ABC-SMC.
• 10. Research directions and open problems — gradient-based samplers, Markovian stochastic approximation, theory of EM, simulation-based inference (self-study).

Session plan

40 h in class = 10 sessions of 4 h (2 h lecture + 2 h lab).

The volume is the Université Paris-Saclay programme table (20 h of lectures and 20 h of labs), which matches the master-mva.com description (20 h of lectures and 20 h of training and lab work). The split into sessions, and the placement of the four homework assignments after four of the labs, is a proposal.

The order follows the official description, with one change: the stochastic versions of EM (official part 1) are taught after simulation and MCMC, because MCMC-SAEM needs both. Level 1 = foundations, 2 = core results and algorithms, 3 = advanced and recent methods; the level never decreases.

Self-study, outside the 40 h (about 5 h in total): ch. 10 Research directions and open problems — Langevin and Hamiltonian samplers, stochastic-gradient MCMC, Markovian stochastic approximation, theory of EM, tempering, simulation-based inference, stochastic algorithms in computational anatomy (≈ 5 h, with the open-ended parts of Ex. 16).

How to study this module

• Read ch. 0 before the first session, even if the tools look familiar: the later proofs cite its numbered results without restating them.
• After each lecture, redo the proofs of the boxed theorems without notes; the written-exam proposal asks for exactly these arguments.
• Run every code block of a chapter during its lab, then change one parameter (step size, proposal scale, tolerance) and predict the effect before running again.
• Do the exercises listed in the Practice column right after the matching session; the ★★★ exercises are good preparation for the paper presentation.
• Python environment: NumPy and SciPy only for the course code; every block fixes its random seed and runs in under 30 seconds on a laptop.

0. Review: Monte Carlo, convexity and martingales

Every algorithm of the course is analysed with four M1 tools. The law of large numbers and the central limit theorem give Monte Carlo its error bar. The Kullback–Leibler divergence makes EM work. Smooth and strongly convex functions make gradient methods converge. Martingales make stochastic approximation converge.

This chapter states each tool in the form used later, with short proofs. Background is in the refreshers Probability and Statistics and Linear Algebra, Analysis and Optimization.

0.1 Monte Carlo integration

import numpy as np

rng = np.random.default_rng(0)
n = 100_000
for d in [1, 10, 50]:
    X = rng.standard_normal((n, d))
    h = np.exp(-0.5 * (X ** 2).sum(axis=1))      # h(x) = exp(-|x|^2 / 2)
    est = h.mean()
    se = h.std(ddof=1) / np.sqrt(n)              # CLT standard error
    exact = 2.0 ** (-d / 2)
    print(f"d={d:2d}  estimate={est:.4e}  exact={exact:.4e}  relative s.e.={se / est:.4f}")

0.2 Jensen’s inequality and the Kullback–Leibler divergence

This inequality is the whole reason EM increases the likelihood (Theorem 3.6).

0.3 Exponential families

Chapters 3 and 7 use this structure: when the complete data of a latent-variable model form an exponential family, the EM algorithm only manipulates the expected sufficient statistics.

0.4 Smooth and strongly convex functions

The stochastic version of this argument is Theorem 1.3; the deterministic theory is developed in Convex Optimization and Applications in Machine Learning.

0.5 Martingales

This is the probabilistic half of the step-size condition \sum_n \gamma_n^2 < \infty in stochastic approximation (ch. 4.2).

1. The stochastic gradient algorithm

Fitting a model means minimising a risk that is an expectation, or an average over n data points with n in the millions. One exact gradient then costs a full pass over the data. The stochastic gradient algorithm replaces it with an unbiased estimate computed from a single sample.

This chapter proves when that works and at what rate, and lists what goes wrong in practice. The general theory of such noisy iterations is ch. 4.

1.1 The algorithm

For least squares, f(\theta, (x, y)) = \frac12 (y - x^\top \theta)^2 and the update reads \theta_t = \theta_{t-1} - \gamma_t\, (x_t^\top \theta_{t-1} - y_t)\, x_t . This is the least-mean-squares rule of adaptive signal processing.

The economics are simple. With n = 10^6 points in dimension d = 100, a full gradient costs about 2 \times 10^8 operations and one stochastic gradient about 200. One step of gradient descent therefore buys a million steps of SGD.

1.2 Convergence for strongly convex objectives

1.3 Convex objectives and averaging

1.4 Numerical bottlenecks

The theorems hide constants that dominate practice. Each bottleneck below has a precise source.

• Step size. For least squares, a step larger than 2/\|x_t\|^2 amplifies the error along x_t; a few such steps can make the iterates explode before 1/t steps become small. Practical schedules use \gamma_t = \gamma_0 / (1 + \gamma_0 \mu t), which caps the first steps.
• Conditioning. For least squares, the error along an eigenvector of the Hessian with eigenvalue \lambda_j contracts like (1 - \gamma\lambda_j)^t. Directions with small \lambda_j are slow; standardising the features and preconditioning reduce \kappa = L/\mu.
• Mini-batches. Averaging b stochastic gradients divides their variance by b and multiplies the cost by b. This pays only when the b gradients are computed in parallel.
• Finite sums. For F = \frac1n \sum_i f_i, variance-reduced methods such as SVRG (Johnson and Zhang, 2013) and SAGA (Defazio, Bach and Lacoste-Julien, 2014) store or periodically recompute full gradients. They recover a linear rate for smooth strongly convex problems.
• Monitoring and non-convexity. F itself is too expensive to evaluate, so one tracks a running average of the sampled losses or a held-out loss. Without convexity, SGD only reaches approximate stationary points (Ghadimi and Lan, 2013).

import numpy as np

rng = np.random.default_rng(1)
d, T = 10, 50_000
eig = np.linspace(0.1, 1.0, d)            # Hessian spectrum: mu = 0.1, L = 1
theta_star = np.ones(d)
X = np.sqrt(eig) * rng.standard_normal((T, d))
Y = X @ theta_star + rng.standard_normal(T)   # noise variance sigma^2 = 1

def excess_risk(theta):                    # F(theta) - F(theta*)
    return 0.5 * np.sum(eig * (theta - theta_star) ** 2)

def sgd(step, average=False):
    theta, avg = np.zeros(d), np.zeros(d)
    for t in range(1, T + 1):
        x, y = X[t - 1], Y[t - 1]
        theta = theta - step(t) * (x @ theta - y) * x    # stochastic gradient step
        avg += (theta - avg) / t                          # running Polyak-Ruppert average
    return excess_risk(avg if average else theta)

mu, L = 0.1, 1.0
runs = {
    "constant step 0.05": (lambda t: 0.05, False),
    "1/(mu (t + 2L/mu))": (lambda t: 1 / (mu * (t + 2 * L / mu)), False),
    "c/t with 2 mu c = 0.2": (lambda t: 1.0 / (t + 2), False),
    "0.5/sqrt(t), averaged": (lambda t: 0.5 / np.sqrt(t), True),
}
for name, (step, avg) in runs.items():
    print(f"{name:24s} excess risk = {sgd(step, avg):.2e}")
print(f"{'statistical limit':24s} sigma^2 d / (2T) = {d / (2 * T):.2e}")

2. Bayesian inference and its computational challenges

Bayesian inference treats the parameter as a random variable and conditions on the data. Every Bayesian answer — an estimator, a credible interval, a prediction, a model comparison — is then an integral against the posterior distribution.

This chapter introduces the framework briefly, then pinpoints which integrals are hard. Each obstacle it identifies is removed by one of the later chapters.

2.1 Prior, likelihood and posterior

The posterior is known up to the constant p(y) as soon as the likelihood and the prior can be evaluated. Much of the course is about doing statistics with a density known only up to a constant.

2.2 Bayes estimators

2.3 Large samples: Bernstein–von Mises and Laplace

The posterior forgets the prior and becomes Gaussian around the maximum-likelihood estimator. Credible sets then become asymptotic confidence sets, and the posterior mean, median and mode agree to first order.

import numpy as np
from scipy import stats, special, integrate

a, b, n, k = 2, 2, 10, 7                     # Beta(2, 2) prior, 7 successes out of 10
post = stats.beta(a + k, b + n - k)          # conjugacy: Beta(9, 5)
mode = (a + k - 1) / (a + b + n - 2)
print(f"posterior mean {post.mean():.4f}  median {post.median():.4f}  "
      f"MAP {mode:.4f}  MLE {k / n:.4f}")
print("95% credible interval", np.round(post.interval(0.95), 3))

# evidence p(y): closed form, adaptive quadrature, Laplace approximation
log_ev = np.log(special.comb(n, k)) + special.betaln(a + k, b + n - k) - special.betaln(a, b)
unnorm = lambda p: stats.binom.pmf(k, n, p) * stats.beta.pdf(p, a, b)
quad, _ = integrate.quad(unnorm, 0, 1)
f = lambda p: np.log(unnorm(p))
h = 1e-4
curv = -(f(mode + h) - 2 * f(mode) + f(mode - h)) / h ** 2   # -(log density)'' at the mode
laplace = unnorm(mode) * np.sqrt(2 * np.pi / curv)
print(f"evidence: exact {np.exp(log_ev):.5f}  quadrature {quad:.5f}  Laplace {laplace:.5f}")

These asymptotics fail exactly where computational statistics is needed: small samples, multimodal posteriors, parameters on a boundary, and latent-variable models whose number of hidden variables grows with n.

2.4 Where the computations get hard

Four obstacles make this ratio hard to compute, and each is removed by a later chapter.

1. Dimension. A grid with m points per axis needs m^d evaluations: 10^{20} for m = 10 and d = 20. Monte Carlo (ch. 0.1) avoids the grid but needs samples from the posterior, which is known only up to p(y). Rejection and importance sampling (ch. 5) work in small dimension; MCMC (ch. 6) works in general; adaptive MCMC and tempering (ch. 8) make it robust.
2. Hidden variables. In many models the likelihood is itself an integral over unobserved variables. Even the maximum-likelihood or MAP estimator then requires EM (ch. 3), and stochastic versions of EM when the E-step is an intractable integral (ch. 7).
3. No likelihood at all. Simulators in population genetics or epidemiology can generate data but cannot evaluate p(y \mid \theta). Approximate Bayesian computation (ch. 9) works with simulations only.
4. Large data. When each evaluation of the likelihood costs n operations, stochastic gradients (ch. 1) and stochastic approximation (ch. 4) replace full passes over the data by single samples.

The constant p(y) cancels in the ratio \pi(\theta' \mid y)/\pi(\theta \mid y). This simple remark is what makes the Metropolis–Hastings algorithm (ch. 6.2) possible.

2.5 Latent-variable models

A Bayesian model is itself a latent-variable model in which \theta is the hidden variable. This is why the same simulation tools serve both for computing posteriors and for maximising likelihoods with hidden variables.

3. Majorize–Minimize and the EM algorithm

When the likelihood is an integral over hidden variables, direct maximisation is awkward, while the complete-data likelihood is often simple. The EM algorithm exploits this: at each step it maximises a surrogate that touches the log-likelihood from below.

Seen from further away, EM is one instance of a general principle, Majorize–Minimize, which we state first.

3.1 The Majorize–Minimize principle

Lange, Hunter and Yang (2000) and Hunter and Lange (2004) survey MM algorithms. Mairal (2015) extends them to large finite sums with incremental surrogates.

3.2 EM as a minorize–maximize algorithm

Consider a latent-variable model (Definition 2.8) and the log-likelihood \ell(\theta) = \log p(y; \theta).

With q_k = p(\cdot \mid y; \theta_k), the lower bound \mathcal L(q_k, \theta) equals Q(\theta \mid \theta_k) plus a term that does not depend on \theta.

EM therefore maximises, at each step, a minorizer of \ell that touches it at \theta_k: it is an MM algorithm for -\ell.

Neal and Hinton (1998) read EM as coordinate ascent on \mathcal L(q, \theta): the E-step maximises over q, the M-step over \theta. This view justifies partial E-steps, incremental variants and variational EM, where q is restricted to a tractable family.

3.3 Exponential-family complete data

Assume that the complete data form a curved exponential family (Definition 0.6), p(y, z; \theta) = h(y, z) \exp\big(\langle \varphi(\theta), S(y, z) \rangle - \psi(\theta)\big), with \varphi : \mathbb{R}^p \to \mathbb{R}^q. Then, up to a constant, Q(\theta \mid \theta_k) = \langle \varphi(\theta), s_k \rangle - \psi(\theta), where s_k is the conditional expectation of S(y, Z) under \theta_k.

The E-step is an expectation of a fixed function; the M-step is usually explicit, as for complete-data maximum likelihood. Stochastic versions of EM (ch. 7) keep the M-step and replace the expectation \bar s by simulation.

3.4 Gaussian mixtures

import numpy as np
from scipy.stats import norm

rng = np.random.default_rng(3)
n = 1000
z = rng.random(n) < 0.3                                   # hidden labels, P(z = 1) = 0.3
y = np.where(z, rng.normal(3.0, 1.0, n), rng.normal(0.0, 1.0, n))

def loglik(w, m, s):
    return np.log(w * norm.pdf(y, m[0], s[0]) + (1 - w) * norm.pdf(y, m[1], s[1])).sum()

w, m, s = 0.5, np.array([1.0, -1.0]), np.array([2.0, 2.0])  # w = weight of component 0
trace = [loglik(w, m, s)]
for it in range(200):
    p0 = w * norm.pdf(y, m[0], s[0])
    p1 = (1 - w) * norm.pdf(y, m[1], s[1])
    r = p0 / (p0 + p1)                                     # E-step: responsibilities
    w = r.mean()                                           # M-step: weighted moments
    m = np.array([r @ y / r.sum(), (1 - r) @ y / (1 - r).sum()])
    s = np.sqrt(np.array([r @ (y - m[0]) ** 2 / r.sum(),
                          (1 - r) @ (y - m[1]) ** 2 / (1 - r).sum()]))
    trace.append(loglik(w, m, s))

trace = np.array(trace)
print("log-likelihood at iterations 0, 1, 5, 20, 200:", np.round(trace[[0, 1, 5, 20, 200]], 2))
print("monotone:", bool(np.all(np.diff(trace) > -1e-9)))
print(f"weights ({w:.3f}, {1 - w:.3f})  means {np.round(m, 3)}  std {np.round(s, 3)}")

The likelihood of a Gaussian mixture is unbounded: centring one component on a data point and letting its variance tend to zero sends \ell to +\infty. EM finds local maxima, not this degenerate supremum. In practice one bounds the variances from below or adds a prior and runs MAP-EM, which adds \log \pi(\theta) to Q.

3.5 Convergence and speed

4. Stochastic approximation

SGD (ch. 1), the stochastic versions of EM (ch. 7) and adaptive MCMC (ch. 8) share one structure: a parameter moves along a noisy estimate of a mean field, with decreasing steps. Stochastic approximation, started by Robbins and Monro in 1951, studies this structure once and for all.

This chapter gives the convergence conditions that the syllabus asks for, with a complete proof in the Lyapunov case, then the rates and the effect of averaging.

4.1 The Robbins–Monro algorithm

4.2 Step-size conditions

The first condition lets the iterates travel any distance: the deterministic part moves by \sum_n \gamma_n h(\theta_{n-1}). The second makes the accumulated noise \sum_n \gamma_n \xi_n converge (Corollary 0.12). Steps \gamma_n = \gamma_0 n^{-a} satisfy both when 1/2 < a \le 1.

4.3 The Robbins–Siegmund lemma

4.4 Almost-sure convergence

4.5 The ODE method

When the mean field has several roots, as for non-convex SGD or for EM-like fixed points, no single quadratic Lyapunov function exists. The standard tool compares the iterates with the solutions of the ordinary differential equation \dot\theta = h(\theta): the steps \gamma_n play the role of time increments, and the noise averages out.

For SGD, V = F and \mathcal L is the set of critical points. For SAEM, V is minus the log-likelihood read through the M-step map (Proposition 7.4).

Assumption 3, stability, is the delicate one in practice. A classical remedy is truncation on randomly varying bounds: the iterate is re-initialised whenever it leaves a compact set, and the compact set grows after each re-initialisation (Andrieu, Moulines and Priouret, 2005).

4.6 Rates and averaging

The optimal step 1/(\lambda n) requires the unknown slope \lambda; when \gamma_0\lambda < 1/2 the error decays only like n^{-\gamma_0\lambda}, the phenomenon of Remark 1.4. Averaging removes this dependence.

import numpy as np

rng = np.random.default_rng(4)
alpha, N = 0.9, 100_000
target = np.log(10.0)                      # 0.9-quantile of Exp(1); f(target) = 0.1
X = rng.exponential(size=N)

def robbins_monro(step, average=False):
    th, avg = 0.0, 0.0
    for n in range(1, N + 1):
        th += step(n) * (alpha - (X[n - 1] <= th))   # H(theta, X) = alpha - 1{X <= theta}
        avg += (th - avg) / n
    return avg if average else th

for name, step, av in [("gamma_n = 1/n", lambda n: 1.0 / n, False),
                       ("gamma_n = 10/n", lambda n: 10.0 / n, False),
                       ("n^(-2/3), averaged", lambda n: n ** (-2 / 3), True)]:
    print(f"{name:20s} error = {robbins_monro(step, av) - target:+.4f}")

5. Simulation of random variables

Every stochastic algorithm of this course consumes random variables: the samples of SGD, the hidden variables drawn by stochastic EM, the proposals of MCMC, the synthetic data sets of ABC. This chapter builds non-uniform random variables from uniform ones, by inversion, transformation and rejection.

It then introduces importance sampling, which changes the sampling distribution instead of the samples.

5.1 Inversion

We assume a generator of independent uniform variables on (0, 1). NumPy’s default_rng uses the PCG64 generator; the Mersenne Twister (Matsumoto and Nishimura, 1998) is another standard choice.

Examples follow directly.

• Exponential. -\log(U)/\lambda has law \mathrm{Exp}(\lambda), since 1 - U is uniform too.
• Cauchy. \tan(\pi(U - 1/2)) is standard Cauchy.
• Discrete laws. For probabilities p_1, \dots, p_K, return the first k with p_1 + \dots + p_k \ge U; a binary search on the cumulative sums costs O(\log K) per draw.

Inversion needs F^{-} in closed form or a cheap numerical inverse. The Gaussian distribution function has neither, which motivates the next section.

5.2 Transformations

Two further transformations are used constantly.

• Gaussian vectors. With \Sigma = LL^\top a Cholesky factorisation, \mu + LZ with Z \sim \mathcal N(0, I_d) has law \mathcal N(\mu, \Sigma).
• Mixtures, by composition. Draw the label z from the weights, then y from component z. This is how the data of Example 3.8 were generated.

5.3 Rejection sampling

import numpy as np

rng = np.random.default_rng(5)
M = np.sqrt(2 * np.e / np.pi)                      # sup f/g, attained at x = 1

def half_normal(size):
    out, tried = [], 0
    while sum(len(o) for o in out) < size:
        x = rng.exponential(size=size)             # proposals from g = Exp(1)
        u = rng.random(size)
        out.append(x[u <= np.exp(-0.5 * (x - 1.0) ** 2)])   # accept if U <= f / (M g)
        tried += size
    return np.concatenate(out)[:size], sum(len(o) for o in out) / tried

x, acc = half_normal(200_000)
z = x * rng.choice([-1.0, 1.0], size=x.size)       # random sign -> N(0, 1)
print(f"acceptance {acc:.4f} (1/M = {1 / M:.4f})")
print(f"mean of |Z| {x.mean():.4f} (sqrt(2/pi) = {np.sqrt(2 / np.pi):.4f}),  var of Z {z.var():.4f}")

Rejection is therefore a tool for low dimension. In this course it samples full conditionals inside Gibbs samplers (ch. 6), and it is the ancestor of ABC (ch. 9). For log-concave densities in dimension one, adaptive rejection sampling (Gilks and Wild, 1992) builds the envelope automatically from tangents.

5.4 Importance sampling

When f is known only up to a constant, the self-normalised estimator \hat I_n^{\mathrm{SN}} = \frac{\sum_i w(X_i)\, h(X_i)}{\sum_i w(X_i)} needs only the unnormalised weights. It is consistent, as a ratio of two convergent averages, with a bias of order 1/n.

The effective sample size (\sum_i w_i)^2 / \sum_i w_i^2, between 1 and n, is the standard diagnostic of weight degeneracy. ABC-SMC (ch. 9.3) relies on it.

import numpy as np
from scipy.stats import norm

rng = np.random.default_rng(5)
n, p = 100_000, norm.sf(4.0)
naive = rng.standard_normal(n) > 4.0
y = rng.normal(4.0, 1.0, n)                        # proposal N(4, 1)
terms = (y > 4.0) * np.exp(-4.0 * y + 8.0)         # h(y) w(y)
print(f"exact {p:.4e} | naive {naive.mean():.4e} ({naive.sum()} hits) | "
      f"IS {terms.mean():.4e} +/- {1.96 * terms.std() / np.sqrt(n):.1e}")

6. Markov chain Monte Carlo: Metropolis–Hastings and Gibbs

Rejection and importance sampling fail in high dimension because they draw independent points from a fixed proposal, which ignores where the target mass lies. Markov chain Monte Carlo draws dependent points instead. Each point is a local modification of the previous one, accepted or rejected so that the target distribution stays invariant.

This chapter builds the two classical samplers of the syllabus, Metropolis–Hastings and Gibbs, and the theory that justifies their averages.

6.1 Markov kernels, invariance and ergodicity

A chain of length n is worth n/\tau_h independent samples: this is its effective sample size (ESS). Positive autocorrelations make \tau_h > 1, and every sampler below is judged by this number.

6.2 The Metropolis–Hastings algorithm

Irreducibility holds, for instance, when q(x, y) > 0 for all x and y in the support of \pi, and aperiodicity holds as soon as rejections have positive probability (Tierney, 1994). Two choices of proposal dominate practice.

• Random-walk Metropolis (RWM). Y = x + \sigma \varepsilon with \varepsilon symmetric. The proposal densities cancel, and \alpha = \min(1, \tilde\pi(y)/\tilde\pi(x)).
• Independent Metropolis–Hastings. q(x, y) = g(y), and \alpha = \min(1, w(y)/w(x)) with w = \tilde\pi/g. The chain is uniformly ergodic if and only if w is bounded (Mengersen and Tweedie, 1996).

Always compute \log \alpha rather than \alpha: densities in high dimension underflow.

The practical rule is to tune the proposal scale to an acceptance rate of about 0.23 in high dimension, and about 0.44 in dimension one. The proposal covariance should be proportional to the target covariance, with factor 2.38^2/d (Gelman, Roberts and Gilks, 1996). Adaptive MCMC (ch. 8) automates both rules.

import numpy as np

rng = np.random.default_rng(6)

def ess(x):
    """Effective sample size: n / (1 + 2 sum of autocorrelations up to the first negative one)."""
    x = x - x.mean()
    n = x.size
    f = np.fft.rfft(x, 2 * n)
    acf = np.fft.irfft(f * np.conj(f))[:n] / (x.var() * n)
    neg = np.argmax(acf[1:] < 0) + 1
    return n / (1 + 2 * acf[1:neg].sum())

d, n = 10, 50_000
log_pi = lambda x: -0.5 * x @ x
for scale in [0.1, 2.38 / np.sqrt(d), 3.0]:
    x, lp, acc, chain = np.zeros(d), 0.0, 0, np.empty(n)
    for t in range(n):
        y = x + scale * rng.standard_normal(d)
        lpy = log_pi(y)
        if np.log(rng.random()) < lpy - lp:        # MH acceptance for a symmetric proposal
            x, lp, acc = y, lpy, acc + 1
        chain[t] = x[0]
    print(f"RWM scale {scale:.3f}: acceptance {acc / n:.3f}, ESS of x_1 = {ess(chain):7.0f} / {n}")

6.3 The Gibbs sampler

import numpy as np

rng = np.random.default_rng(6)

def ess(x):
    x = x - x.mean()
    n = x.size
    f = np.fft.rfft(x, 2 * n)
    acf = np.fft.irfft(f * np.conj(f))[:n] / (x.var() * n)
    neg = np.argmax(acf[1:] < 0) + 1
    return n / (1 + 2 * acf[1:neg].sum())

rho, n = 0.99, 50_000
xs, x, y = np.empty(n), 0.0, 0.0
for t in range(n):
    x = rho * y + np.sqrt(1 - rho ** 2) * rng.standard_normal()   # X | Y = y
    y = rho * x + np.sqrt(1 - rho ** 2) * rng.standard_normal()   # Y | X = x
    xs[t] = x
lag1 = np.corrcoef(xs[:-1], xs[1:])[0, 1]
print(f"lag-1 autocorrelation {lag1:.4f} (theory rho^2 = {rho ** 2:.4f})")
print(f"n / ESS = {n / ess(xs):.0f} (theory (1 + rho^2) / (1 - rho^2) = {(1 + rho ** 2) / (1 - rho ** 2):.1f})")

6.4 Practical diagnostics

The ergodic theorem says nothing about how long to wait. Practice relies on the following diagnostics, which can reveal a failure but never prove convergence.

• Burn-in. Discard an initial segment, until traces of the parameters and of \log\tilde\pi(X_n) look stationary.
• Autocorrelation and ESS. Report the effective sample size of every quantity of interest, not the raw chain length.
• Several chains. Run several chains from over-dispersed starting points. The potential scale reduction factor \hat R compares between-chain and within-chain variances (Gelman and Rubin, 1992). The rank-normalised version of Vehtari and coauthors (2021) recommends \hat R < 1.01.
• Known failure mode. A chain trapped in one mode has stationary-looking traces and a reasonable ESS. Only several chains, or the tempering methods of ch. 8, reveal the missing modes.

7. Stochastic versions of EM

The E-step of EM needs the conditional expectation \bar s(\theta) of the complete-data sufficient statistics (Algorithm 3.7). In a mixture it is a finite sum. In a nonlinear mixed-effects model, a deformable template model or a random-effects logistic regression it is an intractable integral.

The stochastic versions of EM replace this expectation by simulation: exact simulation when possible (ch. 5), MCMC otherwise (ch. 6). Stochastic approximation (ch. 4) then controls the resulting noise. Throughout, the complete data form the exponential family of ch. 3.3, with M-step map \bar\theta.

7.1 Monte Carlo EM

The Monte Carlo error of s_k is of order m_k^{-1/2} at every iteration and does not vanish unless m_k \to \infty. Monotonicity of the likelihood is lost, and the total cost \sum_k m_k is high. Booth and Hobert (1999) increase m_k automatically once the change in \theta falls below the Monte Carlo error.

7.2 Stochastic EM

The sequence (\theta_k) is a Markov chain that never settles at a point: at every iteration, a new draw of the hidden variables moves it. Under conditions it is ergodic, with a stationary distribution concentrated near a maximum of the likelihood, and one estimates \theta by averaging the iterates after a burn-in.

The random perturbations help SEM escape saddle points and poor local maxima, which is why it is used to initialise EM. Exercise 14 shows why it cannot converge pointwise.

7.3 Stochastic approximation EM

With \gamma_k = 1 SAEM is SEM. With decreasing steps it averages the simulations over the iterations, so a single simulation per iteration is enough: the average does the work that the large m_k of MCEM did.

The recursion is a Robbins–Monro algorithm (Definition 4.1) on the statistic s, with mean field h(s) = \bar s\big(\bar\theta(s)\big) - s and martingale increments \xi_k = S(y, z_k) - \bar s(\bar\theta(s_{k-1})). The roots of h are the fixed points of EM written in statistics, s^* = \bar s(\bar\theta(s^*)).

Along the mean ODE \dot s = h(s), the likelihood increases, exactly as along EM iterations. The ODE method (Theorem 4.7) with V = -W then gives the convergence of SAEM.

In practice, one takes \gamma_k = 1 during a first phase of K_1 iterations, which explores like SEM. Then \gamma_k = (k - K_1)^{-a} with 1/2 < a \le 1 makes the iterates settle. Stability is enforced by truncation on random bounds (ch. 4.5) when needed.

7.4 MCMC-SAEM

When p(z \mid y; \theta) cannot be sampled exactly, the simulation step becomes a few steps of an MCMC kernel \Pi_\theta whose invariant distribution is p(\cdot \mid y; \theta): z_k \sim \Pi_{\theta_{k-1}}(z_{k-1}, \cdot) . The chain of hidden variables continues from one iteration to the next, while its target moves slowly with \theta_k.

The noise is now Markovian rather than a martingale increment. The proof decomposes it, through the solution of a Poisson equation, into a martingale increment plus a remainder that vanishes thanks to the slow variation of \theta_k (Benveniste, Métivier and Priouret, 1990).

import numpy as np
from scipy.optimize import minimize
from scipy.special import expit, logsumexp

rng = np.random.default_rng(7)
N, m, mu_true, om_true = 200, 10, -0.5, 1.0
phi_true = rng.normal(mu_true, om_true, N)                 # random effects (never observed)
k = (rng.random((N, m)) < expit(phi_true)[:, None]).sum(1)  # successes per subject

def loglik_y(phi, kk):          # log p(y_i | phi_i) up to a constant
    return kk * phi - m * np.logaddexp(0.0, phi)

# MCMC-SAEM: S(z) = (sum phi_i, sum phi_i^2), M-step mu = S1/N, omega^2 = S2/N - mu^2
K_burn, K = 200, 1000
mu, om2, s1, s2 = 0.0, 4.0, 0.0, 0.0
phi = np.zeros(N)
for it in range(1, K + 1):
    for _ in range(5):                                      # simulation: 5 RWM steps per subject
        prop = phi + 0.8 * rng.standard_normal(N)
        log_r = (loglik_y(prop, k) - loglik_y(phi, k)
                 - ((prop - mu) ** 2 - (phi - mu) ** 2) / (2 * om2))
        phi = np.where(np.log(rng.random(N)) < log_r, prop, phi)
    gamma = 1.0 if it <= K_burn else 1.0 / (it - K_burn)   # stochastic approximation step
    s1 += gamma * (phi.sum() - s1)
    s2 += gamma * ((phi ** 2).sum() - s2)
    mu, om2 = s1 / N, s2 / N - (s1 / N) ** 2               # maximisation
print(f"MCMC-SAEM          : mu = {mu:+.3f}, omega = {np.sqrt(om2):.3f}")

# Reference: exact maximum likelihood, integral over phi_i by Gauss-Hermite quadrature
nodes, weights = np.polynomial.hermite_e.hermegauss(60)   # weight exp(-x^2 / 2)
def neg_loglik(p):
    ph = p[0] + np.exp(p[1]) * nodes                        # change of variables phi = mu + omega x
    terms = loglik_y(ph[None, :], k[:, None]) + np.log(weights)[None, :]
    return -(logsumexp(terms, axis=1) - 0.5 * np.log(2 * np.pi)).sum()
res = minimize(neg_loglik, x0=[0.0, 0.0], method="Nelder-Mead", options={"xatol": 1e-6, "fatol": 1e-8})
print(f"Gauss-Hermite MLE  : mu = {res.x[0]:+.3f}, omega = {np.exp(res.x[1]):.3f}")

8. Adaptive MCMC and parallel tempering

The efficiency of Metropolis–Hastings depends on its proposal. Example 6.9 lost a factor of more than ten with a badly scaled random walk, and Example 8.8 below shows a sampler that never leaves the mode it starts in.

Adaptive MCMC learns the proposal from the chain’s own history. Parallel tempering crosses energy barriers through a ladder of flattened targets. Both are recent advances named in the syllabus, and both rely on stochastic approximation (ch. 4).

8.1 Adaptive Metropolis

Theorem 6.8 and its extensions recommend a Gaussian random-walk proposal with covariance \frac{2.38^2}{d}\Sigma_\pi, where \Sigma_\pi is the covariance of the target. The covariance \Sigma_\pi is unknown, but the chain itself produces samples from which it can be estimated.

With \gamma_n = 1/(n+1), (m_n, \Sigma_n) is essentially the empirical mean and covariance of the past states. This is a Robbins–Monro recursion whose mean field vanishes at the mean and covariance of \pi. The term \varepsilon I_d keeps the proposal non-degenerate.

A second, widespread adaptation tunes the global scale \sigma towards a target acceptance rate \alpha^* = 0.234 (Andrieu and Thoms, 2008; Roberts and Rosenthal, 2009): \log\sigma_{n+1} = \log\sigma_n + \gamma_{n+1}\big(\alpha(X_n, Y_{n+1}) - \alpha^*\big) . It is a Robbins–Monro algorithm for the root of \sigma \mapsto \bar\alpha(\sigma) - \alpha^*, where \bar\alpha is the mean acceptance rate.

import numpy as np

rng = np.random.default_rng(8)

def ess(x):
    x = x - x.mean(); n = x.size
    f = np.fft.rfft(x, 2 * n)
    acf = np.fft.irfft(f * np.conj(f))[:n] / (x.var() * n)
    neg = np.argmax(acf[1:] < 0) + 1
    return n / (1 + 2 * acf[1:neg].sum())

# Target N(0, Sigma) in d = 8 with standard deviations from 0.1 to 10, random rotation
d, n = 8, 30_000
Q, _ = np.linalg.qr(rng.standard_normal((d, d)))
Sigma = Q @ np.diag(np.logspace(-1, 1, d) ** 2) @ Q.T
P = np.linalg.inv(Sigma)
log_pi = lambda x: -0.5 * x @ P @ x
sd = 2.38 ** 2 / d

def metropolis(adaptive):
    x, lp = np.zeros(d), 0.0
    mean, cov = np.zeros(d), np.eye(d)
    chain, acc = np.empty((n, d)), 0
    for t in range(1, n + 1):
        C = sd * (cov + 1e-6 * np.eye(d)) if adaptive else sd * np.eye(d)
        y = x + np.linalg.cholesky(C) @ rng.standard_normal(d)
        lpy = log_pi(y)
        if np.log(rng.random()) < lpy - lp:
            x, lp, acc = y, lpy, acc + 1
        chain[t - 1] = x
        if adaptive:                               # Robbins-Monro updates, step 1/(t+1)
            g = 1.0 / (t + 1)
            diff = x - mean
            mean = mean + g * diff
            cov = cov + g * (np.outer(diff, diff) - cov)
    second = chain[n // 2:]
    u = np.linalg.eigh(Sigma)[1][:, -1]           # direction of largest variance
    return acc / n, ess(second @ u)

for name, ad in [("RWM, proposal (2.38^2/d) I", False), ("Adaptive Metropolis", True)]:
    a, e = metropolis(ad)
    print(f"{name:28s} acceptance {a:.3f}  ESS along the widest axis {e:6.0f} / {n // 2}")

8.2 When adaptation is safe

Each kernel used by an adaptive chain may be \pi-invariant while the adaptive chain is not, because the choice of kernel depends on the past.

Adaptive Metropolis satisfies diminishing adaptation because \gamma_n \to 0: the covariance changes by O(1/n) per step. Containment holds on compact state spaces, and on unbounded spaces under tail conditions (Saksman and Vihola, 2010).

In Example 8.3 the kernel keeps switching at full amplitude, so adaptation never diminishes. A simple safe alternative is to adapt during the burn-in only, then freeze the proposal: the remaining chain is an ordinary MCMC chain.

8.3 Multimodality and tempering

Write \tilde\pi = e^{-U}, with U the energy. A random walk that must cross a barrier of height \Delta U between two modes needs a number of steps of order e^{\Delta U}.

For the mixture \frac12\mathcal N(-5, 1) + \frac12\mathcal N(5, 1), the energy at the midpoint exceeds the energy at a mode by 12.5 - \log 2 \approx 11.8. Crossings are so rare that no realistic run observes one.

When a hot chain sits in a mode that the cold chain has never visited, a swap is accepted with high probability, and the cold chain jumps to that mode at once (Exercise 12).

import numpy as np

rng = np.random.default_rng(9)

def log_pi(x):          # 0.5 N(-5, 1) + 0.5 N(5, 1), up to a constant
    return np.logaddexp(-0.5 * (x + 5.0) ** 2, -0.5 * (x - 5.0) ** 2)

n = 20_000
betas = np.array([1.0, 0.3, 0.1, 0.03])        # inverse temperatures, beta_0 = 1
scales = 1.0 / np.sqrt(betas)                  # wider local moves when hotter

def run(tempering):
    K = len(betas) if tempering else 1
    x = np.full(K, -5.0)                       # every chain starts in the left mode
    cold, swaps, tries = np.empty(n), np.zeros(K - 1), np.zeros(K - 1)
    for t in range(n):
        y = x + scales[:K] * rng.standard_normal(K)          # local RWM move at each level
        acc = np.log(rng.random(K)) < betas[:K] * (log_pi(y) - log_pi(x))
        x = np.where(acc, y, x)
        if tempering:                                        # swap a random adjacent pair
            k = rng.integers(K - 1)
            log_a = (betas[k] - betas[k + 1]) * (log_pi(x[k + 1]) - log_pi(x[k]))
            tries[k] += 1
            if np.log(rng.random()) < log_a:
                x[k], x[k + 1] = x[k + 1], x[k]
                swaps[k] += 1
        cold[t] = x[0]
    return cold, swaps / np.maximum(tries, 1)

for name, tempering in [("single RWM chain", False), ("parallel tempering", True)]:
    cold, rates = run(tempering)
    print(f"{name:20s} fraction of time in the right mode {np.mean(cold > 0):.3f} (target 0.5)"
          + (f", swap rates {np.round(rates, 2)}" if tempering else ""))

The ladder is a compromise. Gaps that are too large make swaps fail; gaps that are too small require many chains. For product targets in dimension d, Atchadé, Roberts and Rosenthal (2011) show that the swap acceptance rate between neighbours should be about 0.234, and the number of levels needed grows like \sqrt d.

8.4 Adaptive parallel tempering

Choosing the ladder by hand requires knowing where the barriers are. Miasojedow, Moulines and Vihola (2013) adapt it with stochastic approximation, together with the local proposals.

The step sizes decrease, so the adaptation diminishes, and the authors prove the convergence of the adapted parameters and a law of large numbers for the cold chain. The result is a sampler whose only inputs are the target, the number of levels and the target acceptance rates.

9. Approximate Bayesian computation

Some models can be simulated but not evaluated. The likelihood of a genealogy in population genetics, of a stochastic epidemic, or of a distribution defined by its quantile function has no tractable expression, so Metropolis–Hastings cannot compute its acceptance ratio.

Approximate Bayesian computation (ABC) replaces the evaluation of the likelihood by simulation and comparison. It is the last of the recent advances in the syllabus, and it builds directly on rejection sampling (ch. 5.3) and MCMC (ch. 6).

9.1 Rejection ABC

For discrete data and \varepsilon = 0 with s the identity, the acceptance probability is exactly p(y \mid \theta) and rejection ABC is exact. For continuous data, the tolerance trades accuracy against acceptance rate.

import numpy as np

rng = np.random.default_rng(10)
n_obs, tau = 20, 3.0                          # y_i ~ N(theta, 1), prior theta ~ N(0, 3^2)
y_obs = rng.normal(1.5, 1.0, n_obs)
s_obs = y_obs.mean()                          # sufficient statistic
post_var = 1.0 / (1.0 / tau ** 2 + n_obs)     # exact conjugate posterior
post_mean = post_var * n_obs * s_obs
print(f"exact posterior: mean {post_mean:.3f}, sd {np.sqrt(post_var):.3f}")

n_sim = 500_000
theta = rng.normal(0.0, tau, n_sim)                                      # 1. prior draws
s_sim = (theta[:, None] + rng.standard_normal((n_sim, n_obs))).mean(1)   # 2. simulated data sets
dist = np.abs(s_sim - s_obs)
for eps in [1.0, 0.3, 0.1, 0.03]:
    keep = dist <= eps                                                   # 3. accept if close
    th = theta[keep]
    print(f"eps = {eps:4.2f}: acceptance {keep.mean():.4f}, ABC mean {th.mean():.3f}, sd {th.std():.3f}")

# Regression adjustment (Beaumont, Zhang and Balding, 2002) at the loosest tolerance
keep = dist <= 1.0
beta = np.polyfit(s_sim[keep] - s_obs, theta[keep], 1)[0]
adj = theta[keep] - beta * (s_sim[keep] - s_obs)
print(f"eps = 1.00 + regression adjustment: mean {adj.mean():.3f}, sd {adj.std():.3f}")

9.2 Tolerance, dimension and summary statistics

With q = 1 the mean squared error decreases like N^{-0.8}, close to the Monte Carlo rate N^{-1}. With q = 20 it decreases like N^{-1/6}. Summaries must therefore be few and informative.

• Semi-automatic summaries. Fearnhead and Prangle (2012) run a pilot simulation, regress \theta on many candidate features of y', and use the fitted regression, an estimate of \mathbb{E}[\theta \mid y], as the summary. This gives exactly q = d summaries.
• Regression adjustment. Beaumont, Zhang and Balding (2002) fit a linear regression of \theta on s(y') among the accepted pairs and correct each accepted parameter: \theta_i^{\mathrm{adj}} = \theta_i - \hat\beta^\top\big(s(y'_i) - s(y)\big) . This removes the first-order effect of the tolerance. In Example 9.3 the model is linear and Gaussian, so the correction is exact.

9.3 ABC-MCMC and ABC-SMC

Rejection ABC proposes from the prior, which wastes simulations when the posterior is much narrower than the prior. Two remedies follow the two ideas of ch. 5–6: Markov chains, and sequential importance sampling.

In the tails of the posterior, almost every simulated data set is rejected, so an ABC-MCMC chain can stick for very long periods.

The weights are those of importance sampling (Definition 5.8): the proposal is the mixture of kernels centred on the previous population. Beaumont and coauthors take for K_t a Gaussian with twice the weighted empirical variance of the previous population. The tolerances can be set adaptively, for instance as a quantile of the previous distances, and the effective sample size of the weights monitors degeneracy.

9.4 Caveats

• Information loss. Without sufficient summaries, \pi_\varepsilon approximates \pi(\theta \mid s(y)), which can be much wider than \pi(\theta \mid y).
• Model choice. Bayes factors computed by ABC with summaries that are not sufficient across models can be inconsistent (Robert, Cornuet, Marin and Pillai, 2011).
• Asymptotics. Frazier, Martin, Robert and Rousseau (2018) show that the ABC posterior concentrates on the true parameter when \varepsilon decreases fast enough with the sample size, and characterise its asymptotic shape.

10. Research directions and open problems

This self-study chapter connects the course to current research. Each section names published work, states what is known, and ends with questions that remain open. It assumes ch. 1–9.

10.1 Gradient-based samplers

The Langevin diffusion dX_t = \nabla \log\pi(X_t)\, dt + \sqrt 2\, dW_t leaves \pi invariant. Its Euler discretisation with step h is the unadjusted Langevin algorithm (ULA): X_{k+1} = X_k + h\, \nabla\log\pi(X_k) + \sqrt{2h}\ \xi_{k+1}, \qquad \xi_{k+1} \sim \mathcal N(0, I_d) . ULA is biased for every h > 0. Using its step as a Metropolis–Hastings proposal removes the bias: this is MALA, the Metropolis-adjusted Langevin algorithm (Roberts and Tweedie, 1996).

• Scaling with the dimension. Random-walk Metropolis needs steps of order d^{-1} (Theorem 6.8), MALA steps of order d^{-1/3}, with an optimal acceptance rate of 0.574 (Roberts and Rosenthal, 1998).
• Non-asymptotic theory. For smooth, strongly log-concave targets, Dalalyan (2017) and Durmus and Moulines (2017) bound the distance between the law of ULA after k steps and \pi, with explicit, polynomial dependence on d.
• Hamiltonian Monte Carlo. HMC (Duane, Kennedy, Pendleton and Roweth, 1987; Neal, 2011) adds a momentum variable and simulates Hamiltonian dynamics. The No-U-Turn sampler (Hoffman and Gelman, 2014) tunes its trajectory length automatically and is the default sampler of the Stan language.

import numpy as np

rng = np.random.default_rng(17)
d, n = 100, 5_000
grad_log_pi = lambda x: -x                          # target N(0, I_d)
log_pi = lambda x: -0.5 * x @ x

def run(kind, h):
    x, acc, first = rng.standard_normal(d), 0, np.empty(n)
    for t in range(n):
        if kind == "RWM":
            y = x + np.sqrt(h) * rng.standard_normal(d)
            log_q = 0.0
        else:                                       # MALA: Langevin proposal + MH correction
            mean_x = x + h * grad_log_pi(x)
            y = mean_x + np.sqrt(2 * h) * rng.standard_normal(d)
            mean_y = y + h * grad_log_pi(y)
            log_q = (-np.sum((x - mean_y) ** 2) + np.sum((y - mean_x) ** 2)) / (4 * h)
        if np.log(rng.random()) < log_pi(y) - log_pi(x) + log_q:
            x, acc = y, acc + 1
        first[t] = x[0]
    return acc / n, np.corrcoef(first[:-1], first[1:])[0, 1]

print("RWM , h = 2.38^2/d     : acceptance %.3f, lag-1 autocorrelation of x_1 %.3f" % run("RWM", 2.38 ** 2 / d))
print("MALA, h = 1.3 d^(-1/3) : acceptance %.3f, lag-1 autocorrelation of x_1 %.3f" % run("MALA", 1.3 * d ** (-1 / 3)))

In dimension 100, the optimally scaled random walk has a lag-1 autocorrelation of 0.995; MALA, at an acceptance rate near 0.59, has 0.80. One gradient evaluation per step buys a large gain in mixing.

Open questions: sharp mixing times for targets that are not log-concave; samplers whose cost does not grow with the separation of the modes; the links between Langevin sampling and the score-based generative models of Diffusion and Multiscale Models in Artificial Intelligence and Stochastic Calculus in Machine Learning.

10.2 Stochastic gradients inside samplers

Stochastic gradient Langevin dynamics (Welling and Teh, 2011) runs ULA with mini-batch estimates of \nabla\log\pi and decreasing steps. It merges ch. 1 and ch. 6: the data size no longer limits the cost of an iteration.

The price is a bias that no accept–reject step can remove, since the acceptance ratio would need the full likelihood. Open question: how to trade bias, variance and cost optimally, and how to certify the accuracy of such samplers on large models.

10.3 Stochastic approximation with Markovian noise

SAEM with MCMC simulation (ch. 7.4), adaptive MCMC (ch. 8) and online EM all run Robbins–Monro recursions whose noise is Markovian and biased, not a martingale increment. Asymptotic theory follows the Poisson-equation route of Benveniste, Métivier and Priouret (1990).

Non-asymptotic bounds are more recent: Karimi, Miasojedow, Moulines and Wai (2019) give finite-time rates for biased stochastic approximation with Markovian noise. Open questions: stability without truncation; explicit rates for MCMC-SAEM in high dimension.

10.4 Theory of EM

Wu’s theorem (Theorem 3.9) only promises stationary points. Balakrishnan, Wainwright and Yu (2017) separate a population EM operator from its sample version and prove geometric convergence inside a basin of attraction, with statistical error bounds.

For a balanced mixture of two Gaussians with known covariance, Daskalakis, Tzamos and Zampetakis (2017) show that EM converges from almost any initialisation. With three or more components, Jin, Zhang, Balakrishnan, Wainwright and Jordan (2016) exhibit bad local maxima, with arbitrarily small likelihood relative to the global one, that attract EM from random initialisations with high probability. Open question: which initialisations or stochastic versions of EM provably avoid such maxima.

10.5 Tempering and multimodality

Parallel tempering with reversible swaps moves information along the ladder diffusively. Syed, Bouchard-Côté, Deligiannidis and Doucet (2022) propose a non-reversible swap scheme whose round-trip time scales better with the number of chains, and tune the ladder from the estimated communication barrier.

Woodard, Schmidler and Huber (2009) show, conversely, that tempering can mix torpidly on some multimodal targets, for instance when the modes have different widths. Open question: dimension-robust guarantees for multimodal sampling.

10.6 Simulation-based inference beyond ABC

Neural methods learn a conditional density estimator of \theta given y from simulated pairs, instead of rejecting simulations (Papamakarios and Murray, 2016). Cranmer, Brehmer and Louppe (2020) survey this field of simulation-based inference.

Within ABC, Bernton, Jacob, Gerber and Robert (2019) replace summary statistics by Wasserstein distances between empirical distributions, and Gutmann and Corander (2016) use Bayesian optimisation to reduce the number of simulations. Open questions: calibration of the resulting posteriors, and robustness when the simulator is misspecified.

10.7 Stochastic algorithms in computational anatomy

Allassonnière, Amit and Trouvé (2007) formulate the estimation of a deformable template, an average image together with its typical deformations, as maximum likelihood in a Bayesian latent-variable model. The deformation of each observed image is the hidden variable.

Allassonnière, Kuhn and Trouvé (2010) estimate such models with an MCMC-SAEM algorithm, with a Gibbs-type sampler for the deformations, and prove its convergence. The whole chain of this course appears there: a Bayesian model (ch. 2), hidden variables (ch. 3), stochastic approximation (ch. 4), MCMC (ch. 6) and SAEM (ch. 7). Related courses: Medical Image Analysis Based on Generative, Geometric and Biophysical Models and Geometry and Shape Spaces.

10.8 Removing the burn-in bias

Every MCMC estimator is biased at finite length, because the chain starts out of equilibrium. Jacob, O’Leary and Atchadé (2020) couple two chains so that they meet exactly after a random time, and build unbiased estimators from the coupled path. Independent replicas can then run in parallel, and their average has a valid confidence interval. Open question: couplings that meet quickly in high dimension.

Exercises

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

1. ★ How many Monte Carlo samples? (ch. 0) You estimate a probability p \approx 0.3 by the frequency of an event in n independent trials. How large must n be for a 95 % confidence half-width of 0.01? Same question for p \approx 10^{-4} and a half-width of 10 % of p.

   Solution. The half-width is 1.96\sqrt{p(1-p)/n} (Proposition 0.2).

   For p = 0.3 and half-width 0.01, n \ge 1.96^2 \times 0.21 / 10^{-4} \approx 8068.

   For p = 10^{-4} and half-width 10^{-5}, n \ge \frac{1.96^2 \times p(1 - p)}{10^{-10}} \approx 3.8 \times 10^6 . The relative error of crude Monte Carlo is about 1/\sqrt{np}: rare events call for importance sampling (Example 5.10).

2. ★ SGD computes a mean. (ch. 1) Take f(\theta, Z) = \frac12(\theta - Z)^2 with \mathbb{E}Z = m and \mathrm{Var}\,Z = \sigma^2. Show that SGD with \gamma_t = 1/t returns the running mean of Z_1, \dots, Z_t from any \theta_0. Compare its exact error with the bound of Theorem 1.3.

   Solution. The gradient is \theta - Z and F is 1-strongly convex, so \gamma_t = 1/(\mu t) = 1/t. The update reads \theta_t = (1 - 1/t)\theta_{t-1} + Z_t/t.

   At t = 1, \theta_1 = Z_1 and \theta_0 is forgotten; by induction, t\theta_t = (t-1)\theta_{t-1} + Z_t = \sum_{s \le t} Z_s. The exact error is \mathbb{E}(\theta_t - m)^2 = \sigma^2/t.

   The bound of Theorem 1.3 is B^2/t, where B^2 must dominate the second moment of the gradient: B^2 \ge \mathbb{E}(\theta - Z)^2 = (\theta - m)^2 + \sigma^2 \ge \sigma^2 . The bound is therefore tight up to the size of the domain.

3. ★ Gamma–Poisson. (ch. 2) Let y_1, \dots, y_n be independent Poisson with mean \lambda, and let the prior be \mathrm{Gamma}(a, b) with density proportional to \lambda^{a-1}e^{-b\lambda}. Find the posterior.

   With a = 2, b = 1, n = 5 and \sum_i y_i = 13, compute the posterior mean, mode and standard deviation, and the maximum-likelihood estimate.

   Solution. The posterior is proportional to the likelihood times the prior, \lambda^{\sum_i y_i} e^{-n\lambda} \cdot \lambda^{a-1} e^{-b\lambda}, so it is \mathrm{Gamma}(a + \sum_i y_i,\ b + n), here \mathrm{Gamma}(15, 6).

   Mean 15/6 = 2.5, mode 14/6 \approx 2.333, standard deviation \sqrt{15}/6 \approx 0.645, and maximum-likelihood estimate 13/5 = 2.6.

   The posterior mean is the weighted average \frac16 \times 2 + \frac56 \times 2.6 of the prior mean and the data mean.

4. ★ Inversion. (ch. 5) Give inversion formulas for the Pareto law, with distribution function 1 - (x_m/x)^\alpha for x \ge x_m; for the Weibull law, with distribution function 1 - \exp(-(x/\lambda)^k); and for the geometric law on \{1, 2, \dots\} with parameter p.

   Solution. Using that 1 - U is uniform: Pareto, x_m U^{-1/\alpha}; Weibull, \lambda(-\log U)^{1/k}.

   Geometric: X = \lceil \log U / \log(1-p) \rceil. Indeed X > k if and only if \log U < k\log(1 - p), since \log(1-p) < 0. This has probability (1-p)^k, the geometric tail.

5. ★ Rejection for Beta laws. (ch. 5) Sample \mathrm{Beta}(2, 2), with density 6x(1-x), by rejection from the uniform proposal: find M and the acceptance rate. Then do the same for \mathrm{Beta}(2.7, 6.3), and propose a better proposal.

   Solution. For \mathrm{Beta}(2,2), M = \max 6x(1-x) = 1.5 at x = 1/2, and the acceptance rate is 2/3.

   For \mathrm{Beta}(2.7, 6.3), M is the density at the mode 1.7/7 \approx 0.243, that is M \approx 2.67, and the acceptance rate is 0.375.

   The proposal \mathrm{Beta}(2, 6), easy to simulate with integer parameters, gives a density ratio proportional to x^{0.7}(1-x)^{0.3}, maximal at x = 0.7. Then M \approx 1.67 and the acceptance rate rises to 0.60.

6. ★ Gibbs on a correlated Gaussian (Python). (ch. 6) Run the Gibbs sampler of Example 6.12 for \rho \in \{0, 0.5, 0.9, 0.99\} and compare n/\mathrm{ESS} with (1 + \rho^2)/(1 - \rho^2). Then explain how to make the sampler efficient for every \rho.

   Solution.

   import numpy as np
   
   rng = np.random.default_rng(11)
   
   def ess(x):
       x = x - x.mean(); n = x.size
       f = np.fft.rfft(x, 2 * n)
       acf = np.fft.irfft(f * np.conj(f))[:n] / (x.var() * n)
       neg = np.argmax(acf[1:] < 0) + 1
       return n / (1 + 2 * acf[1:neg].sum())
   
   n = 20_000
   for rho in [0.0, 0.5, 0.9, 0.99]:
       xs, x, y, c = np.empty(n), 0.0, 0.0, np.sqrt(1 - rho ** 2)
       for t in range(n):
           x = rho * y + c * rng.standard_normal()
           y = rho * x + c * rng.standard_normal()
           xs[t] = x
       print(f"rho = {rho:4.2f}: n/ESS = {n / ess(xs):6.1f}   theory {(1 + rho**2) / (1 - rho**2):6.1f}")

   The run gives about 1.0, 1.7, 9.3 and 104, against 1, 1.67, 9.53 and 99.5 in theory.

   The rotated coordinates (X + Y)/\sqrt2 and (X - Y)/\sqrt2 are independent. A Gibbs sampler in those coordinates, or a single block update of (X, Y), therefore has \tau = 1 for every \rho.

7. ★★ A step size that is too small. (ch. 1) Take F(\theta) = \mu\theta^2/2 with exact gradients and steps \gamma_t = c/t, where 2\mu c < 1. Prove that, for a constant C > 0, |\theta_t| \ge C\,|\theta_0|\, t^{-\mu c} . What is the error after 10^6 iterations when \mu c = 0.1?

   Solution. The iterates are \theta_t = \theta_0 \prod_{s=1}^t (1 - \mu c/s), with factors in [1/2, 1) because \mu c < 1/2.

   For x \in [0, 1/2], \log(1 - x) \ge -x - x^2. Hence the logarithm of the product is at least -\mu c\sum_{s \le t} 1/s - (\mu c)^2 \sum_s 1/s^2, which is at least -\mu c(1 + \log t) - (\mu c)^2\pi^2/6.

   Exponentiating gives the bound with C = e^{-\mu c - (\mu c)^2\pi^2/6}. With \mu c = 0.1, t^{-0.1} = 10^{-0.6} \approx 0.25 after 10^6 steps: the squared error is still of order 0.06\,\theta_0^2.

   With c = 1/\mu, by contrast, the first step lands exactly on the minimiser, and with noise Theorem 1.3 gives an error of order 1/t (Remark 1.4).

8. ★★ Majorize–Minimize for least absolute deviations. (ch. 3) Using Example 3.3, derive the MM iteration for minimising F(\theta) = \sum_i |y_i - x_i^\top\theta| and prove that F decreases. Run three iterations by hand for the median of (1, 2, 10) from \theta_0 = 5.

   Solution. With r_i(\theta) = y_i - x_i^\top\theta and weights w_i = 1/|r_i(\theta_k)|, the surrogate is \frac12\sum_i w_i r_i(\theta)^2 plus a constant. Its minimiser is the weighted least-squares solution \theta_{k+1} = (X^\top W_k X)^{-1} X^\top W_k\, y . F decreases by Proposition 3.2. This is iteratively reweighted least squares; in practice |r_i| is replaced by \max(|r_i|, \delta), which majorizes a smoothed objective.

   For the median (x_i = 1), \theta_{k+1} is the average of the y_i with weights 1/|y_i - \theta_k|.

   From \theta_0 = 5, the iterates are 3.723, 2.820 and 2.296. Meanwhile F goes from 12 through 10.72 and 9.82 to 9.30, towards the minimum F(2) = 9.

9. ★★ EM for censored lifetimes. (ch. 3) Lifetimes T_i \sim \mathrm{Exp}(\lambda) are censored at times c_i: one observes y_i = \min(T_i, c_i) and \delta_i = \mathbf 1\{T_i \le c_i\}. Derive EM with the T_i as complete data, find its fixed point, and compute its rate.

   Iterate from \lambda_0 = 1 with n = 10, r = \sum_i \delta_i = 6 and \sum_i y_i = 15.

   Solution. The complete log-likelihood is n\log\lambda - \lambda\sum_i T_i. By the memoryless property, \mathbb{E}[T_i \mid T_i > c_i] = c_i + 1/\lambda_k. Hence \lambda_{k+1} = \frac{n}{\sum_i y_i + (n - r)/\lambda_k} . The fixed point solves \lambda\sum_i y_i + n - r = n, so \lambda^* = r/\sum_i y_i = 0.4, the closed-form maximum-likelihood estimate.

   The derivative of the EM map at \lambda^* is (n - r)/n = 0.4, the fraction of censored data, as Proposition 3.10 predicts.

   The iterates are 0.526, 0.442, 0.416, 0.406 and 0.402: the error shrinks by a factor approaching 0.4.

10. ★★ Averaging in practice (Python). (ch. 4) Estimate by simulation the variance of \sqrt N(\bar\theta_N - \theta^*) for the averaged quantile recursion of Example 4.10 with steps n^{-2/3}, over 200 replications, and compare with Theorem 4.9.

    Solution.

    import numpy as np
    
    rng = np.random.default_rng(12)
    alpha, N, R = 0.9, 100_000, 200
    target = np.log(10.0)                            # true 0.9-quantile of Exp(1)
    th, avg = np.zeros(R), np.zeros(R)               # R independent replications, vectorised
    for n in range(1, N + 1):
        x = rng.exponential(size=R)
        th += n ** (-2 / 3) * (alpha - (x <= th))
        avg += (th - avg) / n
    print(f"N Var(averaged RM) = {N * np.var(avg - target):.2f}   (theory alpha(1-alpha)/f^2 = 9)")

    The run gives about 8.8, against the optimal value \alpha(1-\alpha)/f(\theta^*)^2 = 0.09/0.01 = 9, which is also the asymptotic variance of the empirical quantile. With only N = 10^4 steps the value is noticeably higher: the initial transient of the average has not yet faded.

11. ★★ Other acceptance rules. (ch. 6) (a) Show that Barker’s acceptance probability r/(1 + r), where r is the Metropolis–Hastings ratio, also gives a \pi-reversible kernel. (b) Compare it with \min(1, r). (c) For the standard Gaussian target and an independent Gaussian proposal \mathcal N(0, \sigma^2), for which \sigma is the chain uniformly ergodic?

    Solution. (a) With a = \pi(x)q(x,y) and b = \pi(y)q(y,x), the flow from x to y is a \cdot \frac{b/a}{1 + b/a} = \frac{ab}{a + b}, which is symmetric in (x, y).

    (b) r/(1+r) \le \min(1, r), so Barker’s rule rejects more often. By Peskun’s ordering (Peskun, 1973), Metropolis–Hastings then has a smaller asymptotic variance for every function.

    (c) The weight w(x) = \varphi(x)/g(x) is proportional to \exp(-x^2(1 - \sigma^{-2})/2). It is bounded if and only if \sigma \ge 1, which is the condition for uniform ergodicity (ch. 6.2). For \sigma < 1, a chain in the tails, where w is huge, almost never leaves.

12. ★★ The swap move in numbers. (ch. 8) In Example 8.8, the cold chain sits at x_0 = -5. Compute the swap acceptance probability with the level \beta_1 = 0.3 when x_1 = 5, and when x_1 = 0. Explain why a direct swap between \beta = 1 and \beta = 0.03 would rarely be accepted.

    Solution. Up to a constant, the log-target is \log\tilde\pi(x) = \log\big(e^{-(x+5)^2/2} + e^{-(x-5)^2/2}\big) . Two cases arise.

    • for x_1 = 5, the two values of \tilde\pi are equal, the ratio is 1, and the swap is always accepted: the cold chain jumps to the other mode;
    • for x_1 = 0, the log-ratio is \log\tilde\pi(0) - \log\tilde\pi(-5) \approx -12.5 + \log 2 = -11.81 , and the acceptance probability is e^{-0.7 \times 11.81} \approx 2.6 \times 10^{-4}.

    A chain at \beta = 0.03 wanders far from both modes, where \log\tilde\pi is very negative. Exchanging it with the cold chain multiplies a large energy difference by 0.97, so the swap is almost never accepted: the intermediate levels are necessary.

13. ★★ ABC for the g-and-k distribution (Python). (ch. 9) The g-and-k distribution is defined by its quantile function Q(u) = A + B\big[1 + 0.8\tanh(g z_u/2)\big](1 + z_u^2)^k\, z_u, \qquad z_u = \Phi^{-1}(u), and has no closed-form density (Allingham, King and Mengersen, 2009). With B = 1 and k = 0.5 known, estimate (A, g) from 200 observations by rejection ABC, using the median and two quantile-based skewness measures as summaries.

    Solution. Sampling is inversion (ch. 5.1) with z_u replaced by a standard Gaussian draw.

    import numpy as np
    
    rng = np.random.default_rng(13)
    
    def gk_sample(A, B, g, k, size):
        z = rng.standard_normal(size)                # z = Phi^{-1}(U): inversion of the quantile function
        return A + B * (1 + 0.8 * np.tanh(g * z / 2)) * (1 + z ** 2) ** k * z
    
    def summaries(y):                                # median and two quantile skewness measures, row-wise
        q = np.sort(y, axis=-1)
        m = y.shape[-1]
        qq = lambda p: q[..., int(p * (m - 1))]
        q25, q40, q50, q60, q75 = (qq(p) for p in (0.25, 0.4, 0.5, 0.6, 0.75))
        return np.stack([q50, (q75 + q25 - 2 * q50) / (q75 - q25),
                         (q60 + q40 - 2 * q50) / (q60 - q40)], axis=-1)
    
    n, B, k = 200, 1.0, 0.5
    y_obs = gk_sample(3.0, B, 2.0, k, n)              # true (A, g) = (3, 2)
    s_obs = summaries(y_obs)
    
    N = 40_000
    A = rng.uniform(0, 10, N); g = rng.uniform(0, 5, N)       # uniform priors
    s = summaries(gk_sample(A[:, None], B, g[:, None], k, (N, n)))
    scale = s.std(axis=0)
    d = np.sqrt((((s - s_obs) / scale) ** 2).sum(axis=1))
    keep = d <= np.quantile(d, 0.01)                           # keep the closest 1 %
    print(f"A: mean {A[keep].mean():.2f}, sd {A[keep].std():.2f} | "
          f"g: mean {g[keep].mean():.2f}, sd {g[keep].std():.2f}  (truth A = 3, g = 2)")

    The ABC posterior has mean about 3.1 for A and 2.4 for g, with standard deviations of about 0.5 and 0.8. The truth lies within one standard deviation.

    The skewness of the quartiles saturates for large g, because \tanh does, which is why an inner skewness measure helps. Regression adjustment (ch. 9.2) sharpens both marginals.

14. ★★★ SEM, SAEM and the Lyapunov function on a mixture. (ch. 7) Consider the mixture, with equal weights, of two unit-variance Gaussians with unknown means \mu_1 and \mu_2. (a) Write S, \bar\theta and \bar s. (b) Compute \langle \nabla W(s), h(s)\rangle of Proposition 7.4 explicitly. (c) Prove that SEM cannot converge almost surely.

    Solution. (a) The complete log-likelihood is \sum_i \sum_j \mathbf 1\{z_i = j\}(y_i\mu_j - \mu_j^2/2) plus a constant. So S = (n_1, \sum_{z_i = 1} y_i, n_2, \sum_{z_i = 2} y_i) and \bar\theta(s) = (s_2/s_1,\ s_4/s_3). The expectation \bar s replaces the indicators by the responsibilities r_i and 1 - r_i.

    (b) Here \varphi(\theta) = (-\mu_1^2/2,\ \mu_1,\ -\mu_2^2/2,\ \mu_2), \qquad H(s) = \mathrm{diag}(-s_1, -s_3) . A direct computation gives \langle \nabla W(s), h(s)\rangle = \frac{\bar s_1^2}{s_1}\big(\mu_1^{\mathrm{EM}} - \mu_1\big)^2 + \frac{\bar s_3^2}{s_3}\big(\mu_2^{\mathrm{EM}} - \mu_2\big)^2 , where \mu^{\mathrm{EM}} is one EM update from \bar\theta(s). It vanishes exactly at the fixed points of EM.

    (c) Suppose that (\theta_k) converges with positive probability. On that event the iterates eventually stay in a compact set, where the responsibilities lie in [\eta, 1 - \eta] for some \eta > 0. Every labelling then has probability at least \eta^n, and two labellings give means at least 2\delta apart, for some \delta > 0. Hence, eventually, P\big(|\theta_k - \theta_{k-1}| > \delta \mid \mathcal F_{k-1}\big) \ge \eta^n . By Lévy’s conditional Borel–Cantelli lemma, jumps larger than \delta occur infinitely often on that event, which contradicts convergence.

15. ★★★ Making a bad adaptation safe. (ch. 8) (a) In Example 8.3, check that each kernel is \pi-reversible, and compute the stationary law of the adaptive chain for \varepsilon = 0.1. (b) Modify the rule so that at time n the state-dependent choice is used only with probability 1/n, and P_b otherwise. Show that Theorem 8.4 applies. (c) Discuss what plays the role of \varepsilon I_d in Algorithm 8.1.

    Solution. (a) Both matrices are symmetric with rows summing to 1, hence reversible for the uniform law. The adaptive chain has \pi'(0) = 1/(1 + 2\varepsilon) = 0.833 instead of 0.5.

    (b) At time n, the chain now uses the state-dependent kernel with probability 1/n, and P_b otherwise. Every kernel used is still one of the two kernels of (a), hence \pi-invariant. Two consecutive kernels differ only if one of them comes from the state-dependent choice, which has probability at most 2/n: adaptation diminishes in probability.

    Containment holds because only two kernels are ever used, and both are uniformly ergodic: P_b reaches \pi in one step, and P_a comes within any fixed distance of \pi after a fixed number of steps, since its second eigenvalue is 1 - 2\varepsilon. Theorem 8.4 then gives convergence to \pi.

    (c) Without \varepsilon I_d, a chain that stays still for a while makes \Sigma_n shrink, smaller proposals make it move less, and the proposal can collapse. The regularisation keeps a minimal proposal scale, which is the containment condition in disguise (Saksman and Vihola, 2010).

16. ★★★ Optimal scaling, empirically (Python). (ch. 6, 8, 10) For random-walk Metropolis on \mathcal N(0, I_{20}) with proposal \mathcal N(x, \frac{\ell^2}{d} I_d), plot the efficiency \mathrm{ESS}/n of x_1 against the acceptance rate as \ell varies. Is the optimum near 0.234? Open part: repeat in dimension 1 and 5, on a correlated Gaussian, and on the mixture of Example 8.8, and compare with MALA (ch. 10.1).

    Hints and expected result.

    import numpy as np
    
    rng = np.random.default_rng(16)
    
    def ess(x):
        x = x - x.mean(); n = x.size
        f = np.fft.rfft(x, 2 * n)
        acf = np.fft.irfft(f * np.conj(f))[:n] / (x.var() * n)
        neg = np.argmax(acf[1:] < 0) + 1
        return n / (1 + 2 * acf[1:neg].sum())
    
    d, n = 20, 20_000
    for ell in [0.5, 1.0, 1.7, 2.38, 3.5, 5.0]:
        x, lp, acc, chain = np.zeros(d), 0.0, 0, np.empty(n)
        z = rng.standard_normal((n, d)) * ell / np.sqrt(d)
        logu = np.log(rng.random(n))
        for t in range(n):
            y = x + z[t]
            lpy = -0.5 * y @ y
            if logu[t] < lpy - lp:
                x, lp, acc = y, lpy, acc + 1
            chain[t] = x[0]
        print(f"ell = {ell:4.2f}: acceptance {acc / n:.3f}, ESS/n = {ess(chain) / n:.4f}")

    The acceptance rates go from about 0.81 to 0.02 as \ell grows. The efficiency peaks near \ell = 2.38, with acceptance about 0.25, on a flat plateau between \ell = 1.7 and 3.5.

    In dimension 1 the optimum moves to an acceptance rate near 0.44. For correlated targets the rule holds only after rescaling by the target covariance, which is what adaptive Metropolis provides. For the mixture no scale works, and tempering is needed.

Projects and assessment

The official assessment, on master-mva.com, has two parts. Homework following each training session, four assignments in total (“(4)” on the page), counts for 2/3 of the grade. The oral presentation of one scientific paper, prepared as team work with no written report, counts for 1/3. The Université Paris-Saclay course detail describes the assessment as “Evaluation Continue non Intégrale”.

The solo subjects below are sized like an extended homework. The group subjects are paper studies for a team of two that end, as in the official format, with an oral presentation and no report. The written exam is a self-assessment proposal: no written exam is announced.

Solo projects

1. [Solo · ≈ 20 h] Step sizes and averaging in SGD. How sensitive is SGD to its step-size constant, and does averaging remove that sensitivity? A paper study of Nemirovski, Juditsky, Lan and Shapiro (2009), with Polyak and Juditsky (1992).
   • Implement: SGD for \ell_2-regularised logistic regression with steps c/t over a grid of c, a constant step, and averaged steps t^{-1/2} and t^{-2/3}.
   • Data: the breast-cancer data set bundled with scikit-learn (569 samples, 30 features, standardised before use); reference optimum by SciPy’s L-BFGS.
   • Deliverables: code, 4–6 page report, 10-minute talk. Chapters: ch. 1, 4.
   • Success criterion: reproduce the t^{-2\mu c} slowdown of Remark 1.4, and show averaged SGD within a factor 2 of the best tuned c/t.
2. [Solo · ≈ 20 h] Optimal scaling of random-walk Metropolis. Does the 0.234 rule hold in moderate dimension and beyond product targets? A paper study of Roberts, Gelman and Gilks (1997).
   • Implement: random-walk Metropolis with an ESS estimator, on Gaussian product targets in dimensions 1, 5, 20 and 100, then on a correlated Gaussian and a two-component mixture.
   • Deliverables: code, 4–6 page report with efficiency-versus-acceptance curves, 10-minute talk. Chapters: ch. 6, 8.
   • Success criterion: the efficiency peak moves from about 0.44 in dimension 1 to about 0.23 in dimension 100, within the Monte Carlo error.
3. [Solo · ≈ 20 h] EM, SEM and SAEM on a hard mixture. Which algorithm finds the global maximum of a mixture likelihood most often, and how fast? A paper study of Celeux and Diebolt (1985) and Delyon, Lavielle and Moulines (1999).
   • Implement: the three algorithms for a one-dimensional three-component Gaussian mixture with overlapping components (for instance equal weights, means 0, 1.5 and 3, unit variances, 500 synthetic points), each from the same 100 random initialisations.
   • Deliverables: code, 4–6 page report, 10-minute talk. Chapters: ch. 3, 7.
   • Success criterion: a table giving, for each algorithm, how often it reaches the best likelihood found and after how many iterations; the empirical linear rate of EM near its limit matches the spectral radius of I_c^{-1} I_m (Proposition 3.10) within 10 %.

Group projects

1. [Group of 2 · ≈ 25 h each] MCMC-SAEM for a pharmacokinetic mixed-effects model. Can MCMC-SAEM recover population parameters and their standard errors in a nonlinear mixed-effects model? A paper study of Kuhn and Lavielle (2004) and Delyon, Lavielle and Moulines (1999).
   • Model and data: a one-compartment model with first-order absorption and log-normal random effects on absorption rate, volume and elimination rate; 20 synthetic data sets of 50 subjects with 10 samples each.
   • Split: member A writes the model, the simulator and the MCMC kernel; member B writes the stochastic approximation, the M-step and the Fisher information by Louis’ formula.
   • Deliverables: code and a 20-minute joint presentation, no report. Chapters: ch. 3, 4, 6, 7.
   • Success criterion: the 95 % intervals cover the true values in at least 17 of the 20 data sets, the range expected for a 95 % coverage.
2. [Group of 2 · ≈ 25 h each] Adaptive parallel tempering. Does an adaptive temperature ladder match a hand-tuned one on strongly multimodal targets? A paper study of Miasojedow, Moulines and Vihola (2013).
   • Implement: parallel tempering with the ladder adaptation of Algorithm 8.9 and adaptive Metropolis at every level.
   • Targets: a two-dimensional mixture of 20 well-separated Gaussians, and a bimodal Gaussian mixture in dimension 10.
   • Split: member A writes the sampler and the local adaptation; member B writes the ladder adaptation and the diagnostics (mode-visit frequencies, swap rates, ESS).
   • Deliverables: code and a 20-minute joint presentation, no report. Chapters: ch. 6, 8. Success criterion: every mode visited with its correct frequency within 10 %.
3. [Group of 2 · ≈ 25 h each] Likelihood-free inference for the g-and-k distribution. How many simulations do rejection ABC, ABC-MCMC and ABC-SMC need for the same accuracy? A paper study of Beaumont, Cornuet, Marin and Robert (2009), with Allingham, King and Mengersen (2009) for the model.
   • Implement: the three samplers with regression adjustment, and an exact-likelihood MCMC baseline that evaluates the density by numerical inversion of the quantile function.
   • Split: member A writes rejection ABC and ABC-MCMC; member B writes ABC-SMC and the exact baseline.
   • Deliverables: code and a 20-minute joint presentation, no report. Chapters: ch. 5, 6, 9.
   • Success criterion: posterior means within 0.1 posterior standard deviation of the exact ones, reported against the number of simulations.

Written exam (proposal)

Two hours, one handwritten A4 sheet allowed, no calculator needed; four independent questions.

1. Question (5 points). State the Metropolis–Hastings algorithm for a target known up to a constant. Prove that its kernel satisfies detailed balance, and explain why the normalising constant is not needed. What goes wrong with a random-walk proposal that is too small, or too large?

   Expected: Algorithm 6.6; the flow \min(\pi(x)q(x,y), \pi(y)q(y,x)) is symmetric and rejections sit on the diagonal (Theorem 6.7); the constant cancels in the ratio. Too small: high acceptance, tiny moves, high autocorrelation. Too large: almost all rejected. Rule of thumb: acceptance near 0.234.

2. Question (6 points). Observations y_1, \dots, y_n come from a mixture of two Poisson laws: \mathrm{Poisson}(\lambda_1) with weight w and \mathrm{Poisson}(\lambda_2) with weight 1 - w. Write the complete-data log-likelihood, the E-step and the M-step of EM, and prove in three lines that the likelihood does not decrease.

   Expected: responsibilities r_i proportional to w\, e^{-\lambda_1}\lambda_1^{y_i} against (1 - w)\, e^{-\lambda_2}\lambda_2^{y_i}; updates w = \frac1n\sum_i r_i and \lambda_1 = \sum_i r_i y_i / \sum_i r_i, and symmetrically for \lambda_2.

   Monotonicity follows from the lower bound \mathcal L(q, \theta) and Gibbs’ inequality (Theorem 3.6).

3. Question (5 points). Let F be \mu-strongly convex on a compact convex set \Theta, with stochastic gradients of second moment at most B^2. Prove that projected SGD with steps 1/(\mu t) satisfies \mathbb{E}\|\theta_t - \theta^*\|^2 \le B^2/(\mu^2 t).

   Expected: non-expansiveness of the projection; the one-step recursion \delta_t \le (1 - 2\mu\gamma_t)\,\delta_{t-1} + \gamma_t^2 B^2 from strong monotonicity and optimality; the induction of Theorem 1.3, including the case t = 1.

4. Question (4 points). (a) Sample \mathrm{Beta}(2, 2) by rejection from the uniform law: give M, the acceptance rate and the mean number of proposals. (b) Show that rejection ABC with discrete data, s the identity and \varepsilon = 0 samples the exact posterior.

   Expected: (a) M = 3/2, acceptance 2/3, mean 3/2 proposals (Theorem 5.5). (b) A draw \theta is accepted with probability P_\theta(Y' = y) = p(y \mid \theta), so accepted values have density proportional to \pi(\theta)p(y \mid \theta) (Proposition 9.2).

Total: 5 + 6 + 5 + 4 = 20 points.

Grading rubric (proposal)

References

Official references. Neither the master-mva.com course page nor the Université Paris-Saclay course detail lists a bibliography for 2026–2027. The works below are standard sources for the topics of the official description.

Textbooks

• C. P. Robert, G. Casella, Monte Carlo Statistical Methods, 2nd ed., Springer, 2004 — simulation, MCMC, stochastic EM (ch. 5–8).
• G. J. McLachlan, T. Krishnan, The EM Algorithm and Extensions, 2nd ed., Wiley, 2008 — ch. 3, 7.
• A. Benveniste, M. Métivier, P. Priouret, Adaptive Algorithms and Stochastic Approximations, Springer, 1990 — ch. 4, 7.
• H. J. Kushner, G. G. Yin, Stochastic Approximation and Recursive Algorithms and Applications, 2nd ed., Springer, 2003 — ch. 4.
• V. S. Borkar, Stochastic Approximation: A Dynamical Systems Viewpoint, Cambridge University Press, 2008 — ch. 4.
• R. Douc, E. Moulines, P. Priouret, P. Soulier, Markov Chains, Springer, 2018 — ch. 6.
• S. Meyn, R. L. Tweedie, Markov Chains and Stochastic Stability, 2nd ed., Cambridge University Press, 2009 — ch. 6.
• C. P. Robert, The Bayesian Choice, 2nd ed., Springer, 2007 — ch. 2.
• A. Gelman, J. B. Carlin, H. S. Stern, D. B. Dunson, A. Vehtari, D. B. Rubin, Bayesian Data Analysis, 3rd ed., CRC Press, 2013 — ch. 2, 6.
• A. W. van der Vaart, Asymptotic Statistics, Cambridge University Press, 1998 — Bernstein–von Mises theorem (ch. 2).
• L. Devroye, Non-Uniform Random Variate Generation, Springer, 1986 — ch. 5.
• D. Williams, Probability with Martingales, Cambridge University Press, 1991 — ch. 0.
• A. S. Nemirovski, D. B. Yudin, Problem Complexity and Method Efficiency in Optimization, Wiley, 1983 — lower bounds (ch. 1).

Stochastic gradient and stochastic approximation (ch. 1, 4)

• H. Robbins, S. Monro, “A stochastic approximation method”, Annals of Mathematical Statistics 22(3), 400–407, 1951.
• H. Robbins, D. Siegmund, “A convergence theorem for non negative almost supermartingales and some applications”, in Optimizing Methods in Statistics, Academic Press, 233–257, 1971.
• B. T. Polyak, A. B. Juditsky, “Acceleration of stochastic approximation by averaging”, SIAM Journal on Control and Optimization 30(4), 838–855, 1992.
• A. Nemirovski, A. Juditsky, G. Lan, A. Shapiro, “Robust stochastic approximation approach to stochastic programming”, SIAM Journal on Optimization 19(4), 1574–1609, 2009.
• E. Moulines, F. Bach, “Non-asymptotic analysis of stochastic approximation algorithms for machine learning”, Advances in Neural Information Processing Systems 24, 2011.
• L. Bottou, F. E. Curtis, J. Nocedal, “Optimization methods for large-scale machine learning”, SIAM Review 60(2), 223–311, 2018.
• R. Johnson, T. Zhang, “Accelerating stochastic gradient descent using predictive variance reduction”, Advances in Neural Information Processing Systems 26, 2013.
• A. Defazio, F. Bach, S. Lacoste-Julien, “SAGA: A fast incremental gradient method with support for non-strongly convex composite objectives”, Advances in Neural Information Processing Systems 27, 2014.
• S. Ghadimi, G. Lan, “Stochastic first- and zeroth-order methods for nonconvex stochastic programming”, SIAM Journal on Optimization 23(4), 2341–2368, 2013.
• M. Benaïm, “Dynamics of stochastic approximation algorithms”, Séminaire de Probabilités XXXIII, Lecture Notes in Mathematics 1709, Springer, 1–68, 1999.
• C. Andrieu, É. Moulines, P. Priouret, “Stability of stochastic approximation under verifiable conditions”, SIAM Journal on Control and Optimization 44(1), 283–312, 2005.

EM, MM and their stochastic versions (ch. 3, 7)

• A. P. Dempster, N. M. Laird, D. B. Rubin, “Maximum likelihood from incomplete data via the EM algorithm”, Journal of the Royal Statistical Society B 39(1), 1–38, 1977.
• C. F. J. Wu, “On the convergence properties of the EM algorithm”, Annals of Statistics 11(1), 95–103, 1983.
• T. A. Louis, “Finding the observed information matrix when using the EM algorithm”, Journal of the Royal Statistical Society B 44(2), 226–233, 1982.
• X.-L. Meng, D. B. Rubin, “Maximum likelihood estimation via the ECM algorithm: a general framework”, Biometrika 80(2), 267–278, 1993.
• R. M. Neal, G. E. Hinton, “A view of the EM algorithm that justifies incremental, sparse, and other variants”, in Learning in Graphical Models, Kluwer, 355–368, 1998.
• R. Varadhan, C. Roland, “Simple and globally convergent methods for accelerating the convergence of any EM algorithm”, Scandinavian Journal of Statistics 35(2), 335–353, 2008.
• K. Lange, D. R. Hunter, I. Yang, “Optimization transfer using surrogate objective functions”, Journal of Computational and Graphical Statistics 9(1), 1–20, 2000.
• D. R. Hunter, K. Lange, “A tutorial on MM algorithms”, The American Statistician 58(1), 30–37, 2004.
• J. Mairal, “Incremental majorization-minimization optimization with application to large-scale machine learning”, SIAM Journal on Optimization 25(2), 829–855, 2015.
• G. C. G. Wei, M. A. Tanner, “A Monte Carlo implementation of the EM algorithm and the poor man’s data augmentation algorithms”, Journal of the American Statistical Association 85(411), 699–704, 1990.
• J. G. Booth, J. P. Hobert, “Maximizing generalized linear mixed model likelihoods with an automated Monte Carlo EM algorithm”, Journal of the Royal Statistical Society B 61(1), 265–285, 1999.
• G. Celeux, J. Diebolt, “The SEM algorithm: a probabilistic teacher algorithm derived from the EM algorithm for the mixture problem”, Computational Statistics Quarterly 2, 73–82, 1985.
• B. Delyon, M. Lavielle, E. Moulines, “Convergence of a stochastic approximation version of the EM algorithm”, Annals of Statistics 27(1), 94–128, 1999.
• E. Kuhn, M. Lavielle, “Coupling a stochastic approximation version of EM with an MCMC procedure”, ESAIM: Probability and Statistics 8, 115–131, 2004.
• O. Cappé, E. Moulines, “On-line expectation–maximization algorithm for latent data models”, Journal of the Royal Statistical Society B 71(3), 593–613, 2009.

Simulation and MCMC (ch. 5, 6)

• M. Matsumoto, T. Nishimura, “Mersenne Twister: a 623-dimensionally equidistributed uniform pseudo-random number generator”, ACM Transactions on Modeling and Computer Simulation 8(1), 3–30, 1998.
• W. R. Gilks, P. Wild, “Adaptive rejection sampling for Gibbs sampling”, Journal of the Royal Statistical Society C (Applied Statistics) 41(2), 337–348, 1992.
• N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, E. Teller, “Equation of state calculations by fast computing machines”, Journal of Chemical Physics 21(6), 1087–1092, 1953.
• W. K. Hastings, “Monte Carlo sampling methods using Markov chains and their applications”, Biometrika 57(1), 97–109, 1970.
• S. Geman, D. Geman, “Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images”, IEEE Transactions on Pattern Analysis and Machine Intelligence 6(6), 721–741, 1984.
• A. E. Gelfand, A. F. M. Smith, “Sampling-based approaches to calculating marginal densities”, Journal of the American Statistical Association 85(410), 398–409, 1990.
• L. Tierney, “Markov chains for exploring posterior distributions”, Annals of Statistics 22(4), 1701–1728, 1994.
• K. L. Mengersen, R. L. Tweedie, “Rates of convergence of the Hastings and Metropolis algorithms”, Annals of Statistics 24(1), 101–121, 1996.
• P. H. Peskun, “Optimum Monte-Carlo sampling using Markov chains”, Biometrika 60(3), 607–612, 1973.
• G. O. Roberts, A. Gelman, W. R. Gilks, “Weak convergence and optimal scaling of random walk Metropolis algorithms”, Annals of Applied Probability 7(1), 110–120, 1997.
• A. Gelman, G. O. Roberts, W. R. Gilks, “Efficient Metropolis jumping rules”, in Bayesian Statistics 5, Oxford University Press, 599–607, 1996.
• J. Diebolt, C. P. Robert, “Estimation of finite mixture distributions through Bayesian sampling”, Journal of the Royal Statistical Society B 56(2), 363–375, 1994.
• A. Gelman, D. B. Rubin, “Inference from iterative simulation using multiple sequences”, Statistical Science 7(4), 457–472, 1992.
• A. Vehtari, A. Gelman, D. Simpson, B. Carpenter, P.-C. Bürkner, “Rank-normalization, folding, and localization: an improved \hat R for assessing convergence of MCMC”, Bayesian Analysis 16(2), 667–718, 2021.

Adaptive MCMC and tempering (ch. 8)

• H. Haario, E. Saksman, J. Tamminen, “An adaptive Metropolis algorithm”, Bernoulli 7(2), 223–242, 2001.
• G. O. Roberts, J. S. Rosenthal, “Coupling and ergodicity of adaptive Markov chain Monte Carlo algorithms”, Journal of Applied Probability 44(2), 458–475, 2007.
• G. O. Roberts, J. S. Rosenthal, “Examples of adaptive MCMC”, Journal of Computational and Graphical Statistics 18(2), 349–367, 2009.
• C. Andrieu, J. Thoms, “A tutorial on adaptive MCMC”, Statistics and Computing 18(4), 343–373, 2008.
• E. Saksman, M. Vihola, “On the ergodicity of the adaptive Metropolis algorithm on unbounded domains”, Annals of Applied Probability 20(6), 2178–2203, 2010.
• R. H. Swendsen, J.-S. Wang, “Replica Monte Carlo simulation of spin-glasses”, Physical Review Letters 57(21), 2607–2609, 1986.
• C. J. Geyer, “Markov chain Monte Carlo maximum likelihood”, in Computing Science and Statistics: Proceedings of the 23rd Symposium on the Interface, 156–163, 1991.
• Y. F. Atchadé, G. O. Roberts, J. S. Rosenthal, “Towards optimal scaling of Metropolis-coupled Markov chain Monte Carlo”, Statistics and Computing 21(4), 555–568, 2011.
• B. Miasojedow, E. Moulines, M. Vihola, “An adaptive parallel tempering algorithm”, Journal of Computational and Graphical Statistics 22(3), 649–664, 2013.

Approximate Bayesian computation (ch. 9)

• S. Tavaré, D. J. Balding, R. C. Griffiths, P. Donnelly, “Inferring coalescence times from DNA sequence data”, Genetics 145(2), 505–518, 1997.
• J. K. Pritchard, M. T. Seielstad, A. Perez-Lezaun, M. W. Feldman, “Population growth of human Y chromosomes: a study of Y chromosome microsatellites”, Molecular Biology and Evolution 16(12), 1791–1798, 1999.
• M. A. Beaumont, W. Zhang, D. J. Balding, “Approximate Bayesian computation in population genetics”, Genetics 162(4), 2025–2035, 2002.
• P. Marjoram, J. Molitor, V. Plagnol, S. Tavaré, “Markov chain Monte Carlo without likelihoods”, Proceedings of the National Academy of Sciences 100(26), 15324–15328, 2003.
• S. A. Sisson, Y. Fan, M. M. Tanaka, “Sequential Monte Carlo without likelihoods”, Proceedings of the National Academy of Sciences 104(6), 1760–1765, 2007.
• M. A. Beaumont, J.-M. Cornuet, J.-M. Marin, C. P. Robert, “Adaptive approximate Bayesian computation”, Biometrika 96(4), 983–990, 2009.
• M. G. B. Blum, “Approximate Bayesian computation: a nonparametric perspective”, Journal of the American Statistical Association 105(491), 1178–1187, 2010.
• P. Fearnhead, D. Prangle, “Constructing summary statistics for approximate Bayesian computation: semi-automatic approximate Bayesian computation”, Journal of the Royal Statistical Society B 74(3), 419–474, 2012.
• J.-M. Marin, P. Pudlo, C. P. Robert, R. J. Ryder, “Approximate Bayesian computational methods”, Statistics and Computing 22(6), 1167–1180, 2012.
• C. P. Robert, J.-M. Cornuet, J.-M. Marin, N. S. Pillai, “Lack of confidence in approximate Bayesian computation model choice”, Proceedings of the National Academy of Sciences 108(37), 15112–15117, 2011.
• D. T. Frazier, G. M. Martin, C. P. Robert, J. Rousseau, “Asymptotic properties of approximate Bayesian computation”, Biometrika 105(3), 593–607, 2018.
• C. Andrieu, G. O. Roberts, “The pseudo-marginal approach for efficient Monte Carlo computations”, Annals of Statistics 37(2), 697–725, 2009.
• D. Allingham, R. A. R. King, K. L. Mengersen, “Bayesian estimation of quantile distributions”, Statistics and Computing 19(2), 189–201, 2009.

Research directions (ch. 10)

• G. O. Roberts, R. L. Tweedie, “Exponential convergence of Langevin distributions and their discrete approximations”, Bernoulli 2(4), 341–363, 1996.
• G. O. Roberts, J. S. Rosenthal, “Optimal scaling of discrete approximations to Langevin diffusions”, Journal of the Royal Statistical Society B 60(1), 255–268, 1998.
• A. S. Dalalyan, “Theoretical guarantees for approximate sampling from smooth and log-concave densities”, Journal of the Royal Statistical Society B 79(3), 651–676, 2017.
• A. Durmus, É. Moulines, “Nonasymptotic convergence analysis for the unadjusted Langevin algorithm”, Annals of Applied Probability 27(3), 1551–1587, 2017.
• S. Duane, A. D. Kennedy, B. J. Pendleton, D. Roweth, “Hybrid Monte Carlo”, Physics Letters B 195(2), 216–222, 1987.
• R. M. Neal, “MCMC using Hamiltonian dynamics”, in Handbook of Markov Chain Monte Carlo, CRC Press, 113–162, 2011.
• M. D. Hoffman, A. Gelman, “The No-U-Turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo”, Journal of Machine Learning Research 15, 1593–1623, 2014.
• M. Welling, Y. W. Teh, “Bayesian learning via stochastic gradient Langevin dynamics”, Proceedings of the 28th International Conference on Machine Learning (ICML), 681–688, 2011.
• B. Karimi, B. Miasojedow, E. Moulines, H.-T. Wai, “Non-asymptotic analysis of biased stochastic approximation scheme”, Proceedings of the 32nd Conference on Learning Theory (COLT), PMLR 99, 1944–1974, 2019.
• S. Balakrishnan, M. J. Wainwright, B. Yu, “Statistical guarantees for the EM algorithm: from population to sample-based analysis”, Annals of Statistics 45(1), 77–120, 2017.
• C. Daskalakis, C. Tzamos, M. Zampetakis, “Ten steps of EM suffice for mixtures of two Gaussians”, Proceedings of the 30th Conference on Learning Theory (COLT), PMLR 65, 704–710, 2017.
• C. Jin, Y. Zhang, S. Balakrishnan, M. J. Wainwright, M. I. Jordan, “Local maxima in the likelihood of Gaussian mixture models: structural results and algorithmic consequences”, Advances in Neural Information Processing Systems 29, 2016.
• S. Syed, A. Bouchard-Côté, G. Deligiannidis, A. Doucet, “Non-reversible parallel tempering: a scalable highly parallel MCMC scheme”, Journal of the Royal Statistical Society B 84(2), 321–350, 2022.
• D. B. Woodard, S. C. Schmidler, M. Huber, “Sufficient conditions for torpid mixing of parallel and simulated tempering”, Electronic Journal of Probability 14, 780–804, 2009.
• G. Papamakarios, I. Murray, “Fast ε-free inference of simulation models with Bayesian conditional density estimation”, Advances in Neural Information Processing Systems 29, 2016.
• K. Cranmer, J. Brehmer, G. Louppe, “The frontier of simulation-based inference”, Proceedings of the National Academy of Sciences 117(48), 30055–30062, 2020.
• E. Bernton, P. E. Jacob, M. Gerber, C. P. Robert, “Approximate Bayesian computation with the Wasserstein distance”, Journal of the Royal Statistical Society B 81(2), 235–269, 2019.
• M. U. Gutmann, J. Corander, “Bayesian optimization for likelihood-free inference of simulator-based statistical models”, Journal of Machine Learning Research 17(125), 1–47, 2016.
• S. Allassonnière, Y. Amit, A. Trouvé, “Towards a coherent statistical framework for dense deformable template estimation”, Journal of the Royal Statistical Society B 69(1), 3–29, 2007.
• S. Allassonnière, E. Kuhn, A. Trouvé, “Construction of Bayesian deformable models via a stochastic approximation algorithm: a convergence study”, Bernoulli 16(3), 641–678, 2010.
• P. E. Jacob, J. O’Leary, Y. F. Atchadé, “Unbiased Markov chain Monte Carlo methods with couplings”, Journal of the Royal Statistical Society B 82(3), 543–600, 2020.