Probability Theory for ML: Applications to Monte Carlo Methods and Generative Models

Contents

• Introduction and overview
• 0. Review: measure, random variables, convergence
• 1. Measure and integration: expectation, Tonelli-Fubini
• 2. Convergence: LLN, CLT, delta method
• 3. Simulating: inverse, rejection, Box-Muller
• 4. Monte Carlo and variance reduction
• 5. Importance sampling
• 6. Conditional expectation
• 7. Markov chains: stationarity, reversibility, ergodicity
• 8. MCMC: Metropolis-Hastings, Gibbs, diagnostics, annealing
• 9. Langevin dynamics and generative models
• 10. Research (Durmus): sampling in high dimension
• Exercises
• Projects and assessment
• References

Introduction and overview

Audience: an engineer who has already drawn Gaussians and estimated a mean by sample average — no need to open on sigma-algebras.

Empirical first: a Monte Carlo estimator that is too noisy, a Metropolis chain that will not mix, a sampler that “works” in dim 2 and dies in dim 100. Only then do we name measure, Markov, Langevin.

Intractable expectations, complex laws, plausible data: from measure to Markov chains, Monte Carlo, MCMC, Langevin. Behind diffusion / score models: samplers and stochastic dynamics. This is the program’s mathematical depth on the randomness side.

What you will be able to do

• Handle expectation and integration in the measure-theoretic sense (Tonelli-Fubini, conditional expectation) without loss of rigor.
• State the convergence results (LLN, CLT, delta method) and use them to build confidence intervals by simulation.
• Analyze a Markov chain: stationarity, reversibility, ergodicity, speed of convergence to the stationary distribution.
• Implement the basic sampling methods: inverse, rejection, Box-Muller.
• Build a Monte Carlo estimator and make it efficient: variance reduction, importance sampling.
• Implement MCMC algorithms — Metropolis-Hastings, Gibbs — with the diagnostics that certify convergence, and use simulated annealing.
• Connect Langevin dynamics to generative models: sample a high-dimensional density by following a stochastic process.

Official objectives and syllabus (2026-2027 catalog)

Official objectives (paraphrased from the description published on SynapseS): the course sets itself three objectives.

First, present the foundations of probability on the basis of abstract measure theory: construction of probability spaces, of probabilities on measured spaces and of integration over abstract spaces, with the essential properties of the integral.

Next, present and analyze Monte Carlo methods and their “Markov chain” version (MCMC). Finally, conclude with a brief introduction to generative models.

By the end of the course, a student can:

• apply the fundamental concepts of probability theory and situate it within modern statistics, machine learning and applied contexts;
• apply the basic Monte Carlo methods;
• solve basic problems in machine learning and computational statistics that involve probability;
• solve complex problems involving stochastic processes.

Official syllabus “Plan for the Course” (9 lectures) and mapping to our chapters:

All official themes are covered. Two chapters go beyond the official syllabus: ch. 9 treats Langevin dynamics in much more depth than the announced “brief introduction”, and ch. 10 (research of the course coordinator: sampling in high dimension) is an addition of these notes, as is the emphasis on MCMC diagnostics and simulated annealing in ch. 8.

Where this course fits in the program

Notes kept for reference: this course unit is not in the LLGA 2026–2027 lists (it is still scheduled in Period 1 for other tracks of the École, 36 h), so no LLGA module requires it. For an LLGA student it is a mathematical deepening rather than an application course.

The official description requires an undergraduate-level background in statistics, probability, linear algebra and calculus; for students of the X engineering cycle, the sheet additionally requires having taken at least one course among MAP432, MAP433.

In the master’s program, its natural prerequisite is the Refresher in Statistics module, whose measure theory, random variables and convergence it takes up to raise them to the level required by simulation: the review chapter ensures the transition.

No Python prerequisite is required: the official sheet specifies that a mini-training and a tutorial are provided at the beginning of the course.

Three LLGA modules use material from these notes; none of them assumes the course.

Computational Optimal Transport for ML and Fairness defines the Wasserstein distance (its ch. 2) in which Theorem 9.5 measures the error of ULA, and studies the sampling error of empirical measures (its ch. 4); Multimodal Generative AI learns, in its ch. 3, the scores that our ch. 9 samples with Langevin dynamics; Reinforcement Learning and Autonomous Agents builds on Markov chains and sampling (Monte Carlo methods for control and evaluation).

For students aiming at a PhD in applied probability or in generative modeling, these notes are a strong mathematical complement to the Year 1 courses.

Course outline

• Review: measure, random variables, convergence: the foundation of the Refresher in Statistics, taken up at the theoretical level required afterwards.
• Measure and integration: expectation, Tonelli-Fubini: constructing expectation properly and swapping limits and integrals safely.
• Convergence: LLN, CLT, delta method: what simulation guarantees as the number of samples grows.
• Simulating: inverse, rejection, Box-Muller: generating standard distributions, with their proofs and their limits.
• Monte Carlo and variance reduction: estimating an integral by simulation, then making it converge faster.
• Importance sampling: targeting the regions that matter to reduce the variance of an estimator.
• Conditional expectation: the best possible prediction in the least-squares sense — a central object for generative modeling and for RL alike.
• Markov chains: stationarity, reversibility, ergodicity: the theory underlying all the samplers in the course.
• MCMC: Metropolis-Hastings, Gibbs, diagnostics, annealing: sampling intractable distributions, and knowing when the algorithm has converged.
• Langevin dynamics and generative models: from the stochastic process to sampling high-dimensional densities — the gateway to score-based models.
• Research (Durmus): sampling in high dimension: Alain Durmus’s research on MCMC, Langevin, and generative modeling.

Session plan

36 h in class = 9 lectures of 2 h + 9 exercise sessions of 2 h. The volume and its split are those of the SynapseS 2026–2027 sheet (official teaching methods: 9 lectures of 2 h and 9 exercise sessions of 2 h, some of them on laptop); the sheet gives the theme of each lecture but no chapter-level programme, so the assignment of chapters below is our proposal.

Each row is one week: a 2 h lecture followed by its 2 h exercise session.

Level 1 = foundations (reviews, definitions, first worked examples), 2 = core methods (Monte Carlo, conditioning, Markov chains, MCMC), 3 = advanced material (Langevin dynamics and its link with generative models); the level never goes down from one week to the next.

Self-study, outside the 36 h (about 3 h in total): ch. 10 Research (Durmus) — three papers on the non-asymptotic analysis of Langevin algorithms (≈ 3 h).

How to study this module

• First secure the review material: the opening chapter states precisely what is taken up from the Refresher in Statistics module and at what level.
• Do the proofs on scratch paper (Tonelli-Fubini, ergodicity): in this module, one proof understood saves ten computational exercises.
• Implement inverse, rejection, and Metropolis-Hastings yourself before using any library: these algorithms fit in a few dozen lines and lend themselves well to code.
• Practice MCMC diagnostics (trajectories, autocorrelation, multiple runs): being able to say “this sampler has not converged” is an exam and research skill.
• Play with dimension: redo a Langevin experiment in increasing dimension to get a feel for the difficulty of high-dimensional sampling, the heart of the final chapter’s work.
• The module is demanding upstream (measure) and abstract downstream (Langevin): schedule review time rather than moving ahead superficially.

0. Review: measure, random variables, convergence

Every chapter of this module computes, bounds or estimates an expectation, so it needs a precise notion of probability space, random variable and distribution, together with two elementary inequalities. This review fixes this vocabulary at the level used in the proofs.

1. Measure and integration: expectation, Tonelli-Fubini

Concrete scene: you estimate \mathbb{E}[f(X)] by the sample mean of N=10^4 draws. In dim 1 it converges; in dim 50 almost all points fall outside the useful region and the estimator freezes. This is not a Python bug: it is high-dimensional integration.

This chapter builds expectation properly (Lebesgue, Tonelli-Fubini) so you know when you may swap limit and integral — the reflex before Monte Carlo / MCMC.

2. Convergence: LLN, CLT, delta method

A Monte Carlo estimate is an average of random draws; to trust it, one must know in what sense it approaches the true value, and how large its error is for a given number of draws. This chapter defines the modes of convergence, proves the law of large numbers and the central limit theorem, and derives the delta method, which transfers the CLT to smooth functions of an average.

3. Simulating: inverse, rejection, Box-Muller

A computer only provides uniform random numbers; every simulation algorithm must turn them into draws from a prescribed distribution. This chapter gives three exact ways of doing so — inverting the distribution function, rejection, and the Box-Muller change of variables — with their proofs and their costs.

import numpy as np
rng = np.random.default_rng(42)

def box_muller(n):                      # n standard Gaussians (n even)
    u = rng.random(n // 2); v = rng.random(n // 2)
    r = np.sqrt(-2.0 * np.log(u)); t = 2.0 * np.pi * v
    return np.concatenate([r * np.cos(t), r * np.sin(t)])

def gaussian_rejection(n, m=1.3155):    # N(0,1) by rejection with Laplace
    out = []
    while len(out) < n:
        x = rng.laplace(0.0, 1.0)       # density (1/2) exp(-|x|)
        u = rng.random()
        if u < np.exp(-x*x/2.0) / (m * 0.5 * np.exp(-abs(x))):
            out.append(x)
    return np.array(out)

z1, z2 = box_muller(200000), gaussian_rejection(200000)
print(z1.mean(), z1.var())              # ~0.000 , ~1.000
print(z2.mean(), z2.var())              # ~0.001 , ~1.000

def mixture(n):                          # 0.4 N(-2,0.5^2) + 0.6 N(2.5,0.8^2)
    j = rng.random(n) < 0.6
    return np.where(j, 2.5 + 0.8*rng.standard_normal(n),
                       -2.0 + 0.5*rng.standard_normal(n))

4. Monte Carlo and variance reduction

Many quantities of interest are expectations that cannot be computed in closed form: an integral in high dimension, a probability of failure, the mean of a posterior distribution. Monte Carlo replaces the expectation by an average of simulated draws.

This chapter computes its error, explains why that error does not depend directly on the dimension, and presents two ways of reducing it at almost no cost: antithetic variables and control variates.

import numpy as np
rng = np.random.default_rng(7)

n = 500_000
x, y = rng.uniform(-1, 1, (2, n))
pi_hat = 4.0 * np.mean(x*x + y*y <= 1.0)
se = 4.0 * np.std(x*x + y*y <= 1.0, ddof=1) / np.sqrt(n)
print(f"pi = {pi_hat:.5f} +/- {1.96*se:.5f}")

# antithetic for ln 2
u = rng.random(n // 2)
est = 0.5 * ((1/(1+u)) + (1/(2-u)))
print(f"ln2 = {est.mean():.6f} (variance per draw {est.var(ddof=1):.2e})")

# control variate for E[e^U]
u = rng.random(n)
e = np.exp(u); beta = np.cov(e, u, ddof=1)[0, 1] / u.var(ddof=1)
mc = e - beta * (u - 0.5)
print(f"E[e^U] = {mc.mean():.6f}, e-1 = {np.e-1:.6f}")

5. Importance sampling

When the event of interest is rare, or when one can only sample from a distribution other than the target, plain Monte Carlo wastes almost all its draws. Importance sampling draws from a better-chosen distribution and corrects with weights.

This chapter computes its variance and its optimal proposal, and quantifies the weight degeneracy that limits it in high dimension.

import numpy as np
from scipy import stats
rng = np.random.default_rng(2024)

n = 200_000
p_exact = 1 - stats.norm.cdf(4.0)               # 3.167e-5

mc = rng.standard_normal(n) > 4.0               # crude Monte Carlo
g = rng.normal(4.0, 1.0, n)                     # exponential tilt
w = np.exp(-4.0*g + 8.0)                        # f/g, N(0,1)/N(4,1)
isw = w * (g > 4.0)

for name, est in [("crude", mc.mean()), ("IS", isw.mean())]:
    print(f"{name:5s} {est:.3e}  (exact {p_exact:.3e})")

ess = n / (1.0 + w.var(ddof=1))                 # ESS on the overall weights
print(f"ESS = {ess:.1f} / {n}")                 # weights highly degenerate outside the tail

6. Conditional expectation

A prediction can only use the information available, and an estimator can often be improved by averaging out what is irrelevant to it. Conditional expectation formalizes both ideas: it is the best prediction of a variable given some information.

This chapter defines it in full generality and proves the tower property and the law of total variance, the two rules used by variance reduction and by Markov chains.

7. Markov chains: stationarity, reversibility, ergodicity

When independent draws from a distribution are out of reach, one can still build a random sequence whose distribution converges to it: a Markov chain. This chapter sets up Markov kernels and invariant and reversible distributions, and states the ergodic theorem, which turns averages along a single trajectory into estimates of expectations.

8. MCMC: Metropolis-Hastings, Gibbs, diagnostics, annealing

This chapter puts the two previous ideas together: a Markov chain (ch. 7) built so that its invariant distribution is a prescribed target known only up to a constant, whose trajectory averages replace the independent averages of Monte Carlo (ch. 4).

It presents the Metropolis-Hastings and Gibbs algorithms, proves that they leave the target invariant, and explains how to diagnose their convergence.

import numpy as np
rng = np.random.default_rng(0)

# target: bimodal mixture 0.5 N(-2.5, 0.7^2) + 0.5 N(2.5, 0.7^2)
logpi = lambda x: np.logaddexp(-((x+2.5)**2)/(2*0.49),
                               -((x-2.5)**2)/(2*0.49))

def rw_mh(x0, n, s):
    x = np.empty(n); x[0] = x0; acc = 0
    for k in range(n-1):
        y = x[k] + s * rng.standard_normal()
        if np.log(rng.random()) < logpi(y) - logpi(x[k]):
            x[k+1] = y; acc += 1
        else:
            x[k+1] = x[k]
    return x, acc/n

x, rate = rw_mh(0.0, 200_000, 4.0)     # large step to jump between modes
print(f"acceptance {rate:.3f}, mean {x[1000:].mean():+.3f} (exact 0)")

def ess(z):                            # truncated positive autocorrelations
    z = z - z.mean(); n = len(z)
    acf = [1.0] + [np.corrcoef(z[:-k], z[k:])[0,1] for k in range(1, 200)]
    tau, s = 1.0, 0.0
    for k in range(1, len(acf), 2):    # Geyer pairs
        pair = acf[k] + acf[k+1] if k+1 < len(acf) else acf[k]
        if pair <= 0: break
        s += pair
    return n / (1 + 2*s)

print(f"ESS = {ess(x[1000:]):.0f} / 199000")

9. Langevin dynamics and generative models

Random-walk proposals explore a high-dimensional space slowly because they ignore the shape of the target. Langevin dynamics follows the gradient of the log-density toward high-probability regions, while noise keeps it exploring.

This chapter shows that the target is invariant for the continuous dynamics, computes the bias of its discretization (ULA), and explains why the same equation drives score-based generative models.

import numpy as np
rng = np.random.default_rng(3)

def ula(gamma, n=200_000):             # target N(0,1): U(x) = x^2/2
    x = 3.0
    for k in range(n):
        x = x - gamma*x + np.sqrt(2*gamma)*rng.standard_normal()
        if k >= n//2:                  # ignore burn-in
            yield x

for g in (0.5, 0.1, 0.01):
    v = np.fromiter(ula(g), float)
    print(f"gamma={g:4.2f}  stationary var {v.var():.4f}"
          f"  theory {2/(2-g):.4f}")
# at gamma=0.01 the chain decorrelates in ~1/gamma steps: expect a few % of noise

def mala(gamma, n=200_000):            # same proposal + MH acceptance
    x, out = 3.0, []
    U, dU = lambda z: z*z/2, lambda z: z
    q = lambda a, b: -(b - a + gamma*dU(a))**2 / (4*gamma)   # log q(a -> b) + const
    for k in range(n):
        y = x - gamma*dU(x) + np.sqrt(2*gamma)*rng.standard_normal()
        loga = U(x) - U(y) + q(y, x) - q(x, y)   # log[pi(y) q(y->x)] - log[pi(x) q(x->y)]
        if np.log(rng.random()) < loga: x = y
        if k >= n//2: out.append(x)
    return np.array(out)

print(f"MALA gamma=0.5 : var {mala(0.5).var():.4f} (exact 1)")   # ~0.995

10. Research (Durmus): sampling in high dimension

This self-study chapter presents three papers of the course coordinator that give non-asymptotic guarantees for the Langevin algorithms of chapter 9.

Exercises

1. ★ Tonelli and exact computation. Let U_1, U_2 be i.i.d. \mathcal{U}[0,1]. Writing \mathbb{E}[Y] = \int_0^1 \mathbb{P}(Y > t)\,dt for Y \geq 0 (layer-cake formula), compute \mathbb{E}[\min(U_1,U_2)] then \mathbb{E}[|U_1 - U_2|] via a double integral (Fubini, nonnegative or integrable function).

   Solution. \mathbb{P}(\min > t) = (1-t)^2 so \mathbb{E}[\min] = \int_0^1 (1-t)^2 dt = \tfrac13. For the difference: \mathbb{E}|U_1-U_2| = 2\int_0^1\!\!\int_0^{u_1}(u_1-u_2)\, du_2\,du_1 = 2\int_0^1 \tfrac{u_1^2}{2}du_1 = \tfrac13 (the factor 2 by symmetry of the two half-triangles); check via layers: \mathbb{P}(|U_1-U_2| > t) = (1-t)^2, same integral.

2. ★ Walk on the triangle, then on the square. A random walk jumps uniformly to one of the two other vertices of a triangle, so its matrix is P = (J - I)/2 with J the all-ones matrix. Show that the chain is irreducible and aperiodic, find its invariant distribution and the speed of convergence, and compare with the lazy version (stay in place with probability 1/2). Then do the same for the walk on the square (the 4-cycle).

   Solution. Triangle: the graph is connected, so the chain is irreducible. It returns to its start in 2 steps (A→B→A) and in 3 steps (A→B→C→A), and \gcd(2,3) = 1: aperiodic. P is symmetric and doubly stochastic, so the uniform distribution is invariant and reversible.

   The spectrum of P is \{1, -\tfrac12, -\tfrac12\}, and P^n = J/3 + (-\tfrac12)^n (I - J/3). From any vertex, \|P^n(x,\cdot) - \pi\|_{\mathrm{TV}} = \tfrac23\,(\tfrac12)^n . The lazy chain (I+P)/2 has spectrum \{1, \tfrac14, \tfrac14\} and distance \tfrac23(\tfrac14)^n: here laziness is not needed, and it even speeds convergence up.

   Square: each step changes the parity of the vertex, so returns take an even number of steps: period 2, eigenvalue -1, and P^n(x,x) = 0 for every odd n, so P^n(x,\cdot) does not converge. The lazy walk has spectrum \{1, \tfrac12, \tfrac12, 0\} and converges to the uniform distribution (Thm 7.3), at distance (\tfrac12)^{n+1} from any vertex.

3. ★ Optimal rejection (full computation). We want to simulate the distribution with density f(x) \propto e^{-x^2/2}\mathbf{1}_{x>0} (truncated positive Gaussian) by rejection with the standard exponential proposal g(x) = e^{-x}\mathbf{1}_{x>0}. Compute the explicit M, the acceptance rate, and the expected number of proposals.

   Solution. f(x) = \sqrt{2/\pi}\,e^{-x^2/2}\mathbf{1}_{x>0}; \frac{f}{g} = \sqrt{2/\pi}\,e^{-x^2/2+x}, maximal at x=1 where the exponent equals 1/2.

   Hence M = \sqrt{2e/\pi} \approx 1.3155; rate = 1/M \approx 76\% — the same constants as Example 3.3 (no coincidence: by symmetry, drawing the absolute value of a truncated N(0,1) is equivalent to the symmetric Gaussian/Laplace problem); expected number of trials = M \approx 1.32.

4. ★ Gaussian degeneracy of IS. Redo the computation of Prop. 5.2: show \mathbb{E}_g[w^2] = e^{\mu^2} for f=\mathcal{N}(0,1), g=\mathcal{N}(\mu,1), deduce \mathrm{ESS}/n = e^{-\mu^2}, then comment on the d-dimensional case with a proposal centered at \mu\mathbf{1}.

   Solution. w^2 = e^{-2\mu Y + \mu^2} and \mathbb{E}_g[e^{tY}] = e^{t\mu + t^2/2} with t = -2\mu gives \mathbb{E}_g[e^{-2\mu Y}] = e^{-2\mu^2+2\mu^2} = 1, hence \mathbb{E}_g w^2 = e^{\mu^2}.

   In d product dimensions, \mathbb{E}w^2 = e^{d\mu^2}: \mathrm{ESS}/n = e^{-d\mu^2} — for \mu = 0.5 and d = 50, e^{-12.5} \approx 3.7\times10^{-6}: one million draws yield ~4 effective points.

5. ★ ULA on a Gaussian (proof to complete). For \pi = \mathcal{N}(0,1), show that the stationary distribution of the AR(1) X_{k+1} = (1-\gamma)X_k + \sqrt{2\gamma}Z_{k+1} is \mathcal{N}(0, 2/(2-\gamma)) (stationary variance computation), deduce the relative bias of the variance and the step needed for a bias < 1%.

   Solution. The stationary variance satisfies (independence of the innovations) \sigma^2_\infty = (1-\gamma)^2\sigma^2_\infty + 2\gamma, \quad \text{so} \quad \sigma^2_\infty = 2\gamma/(\gamma(2-\gamma)) = 2/(2-\gamma). Relative bias = \gamma/(2-\gamma) \approx \gamma/2; < 1% requires \gamma < 0.02. Generalization (bonus question): in d isotropic dimensions, same computation per axis; see paper [3] of ch. 10 (tamed ULA) for drifts that are not globally Lipschitz.

6. ★★ Control tower. Prove the tower property in the case \mathcal{H} = \sigma(Y), \mathcal{G} = \sigma(Y,Z) discrete: verify point by point the tower identity by summing over the values of Z: \mathbb{E}[\mathbb{E}[X \mid Y,Z] \mid Y] = \mathbb{E}[X \mid Y]. Recover the law of total probability. Hint. Use Definition 6.1 with A = \{Y = y\} and the partition \{Y=y\} = \bigcup_z \{Y=y, Z=z\}.

   \mathbb{E}[X \mid Y=y] = \sum_z \mathbb{E}[X \mid Y=y, Z=z]\,\mathbb{P}(Z=z \mid Y=y): the average of conditional averages, weighted by the distribution of Z given Y. Trivial \mathcal{H} case: \mathbb{E}[\mathbb{E}[X\mid Y]] = \mathbb{E}[X] = total probability.

7. ★★ MH-independence. Establish the acceptance probability of the independence sampler (q(x,\cdot) = g) and verify detailed balance directly. What happens when g has lighter tails than \pi (e.g. g Gaussian, \pi Cauchy)? Hint. Redo the proof of Thm 8.2 with \rho(x,y) = \pi(x)g(y).

   \alpha(x,y) = \min\big(1, \tfrac{\pi(y)g(x)}{\pi(x)g(y)}\big). Detail: \pi(x)g(y)\,\alpha(x,y) = \min(\pi(x)g(y), \pi(y)g(x)) is symmetric.

   If g decays faster than \pi, the ratio \pi/g is unbounded: rejections dominate in the tails, the chain sticks to its states and is no longer geometrically ergodic — golden rule: heavier tails for the proposal.

8. ★★ ESS and weights. Show that for normalized weights w_1,\dots,w_n > 0, \bar w_i = w_i/\sum_j w_j, the quantity 1/\sum_i \bar w_i^2 \in [1,n] (participation index) coincides with the definition \mathrm{ESS} = n/(1+\mathrm{cv}^2) in expectation, and interpret the two extremes.

   Solution (sketch). Expanding \sum_i \bar w_i^2 to first order around the mean weights, \mathbb{E}\big[\sum_i \bar w_i^2\big] \approx \frac1n\big(1 + \mathrm{cv}^2\big), hence the identification with the definition \mathrm{ESS} \approx 1/\sum_i \bar w_i^2.

   Extremes: equal weights (\bar w_i = 1/n) give \sum \bar w_i^2 = 1/n, i.e. \mathrm{ESS} = n (all draws count); if a single particle carries all the weight, \sum \bar w_i^2 = 1 and \mathrm{ESS} = 1. This index is the diagnostic displayed by all SMC software.

Projects and assessment

Official assessment: “Exam and Lab” — the SynapseS 2026–2027 sheet of APM_51056_EP gives a numerical grade out of 20, a resit (the better of the two grades counts) and acquisition at 10/20, for the degrees it lists. The LLGA lists do not include this unit, so no LLGA rule applies. Everything below is a proposal sized for a 36 h course: solo subjects of about 20 h, group subjects of 25–30 h per member, a three-hour written exam and a rubric.

The projects exercise two master objectives: O1 (foundations: Monte Carlo error bars, the LLN and the CLT, convergence of Markov chains) and O2 (generative modelling theory: MCMC, Langevin dynamics, the score). The score-based subject continues in Multimodal Generative AI, and the molecule subject adds O4 (graph-structured data) through Graph ML. All of them train O7 (hands-on practice): every sampler is first checked on a target with a known answer, and every number comes with its error bar.

Solo projects

1. [Solo · ≈ 20 h] Diagnosing MCMC on multimodal targets. Which diagnostics catch a sampler that has not converged, which ones are fooled, and how does the answer change with the dimension?
   • Implement: random-walk MH, MALA and parallel tempering on mixtures of two well-separated Gaussians in dimension 2, 10 and 50; acceptance rate, ESS by Geyer’s initial sequence and by batch means, \hat R over 4 overdispersed chains, mode-switching rate.
   • Paper study: Vehtari, Gelman, Simpson, Carpenter and Bürkner, “Rank-normalization, folding, and localization: an improved \hat R for assessing convergence of MCMC”, Bayesian Analysis 16(2), 2021 — compare the classic and the rank-normalized \hat R on the same chains.
   • Reproduce: Example 8.7 from its code, and the MALA check of ch. 9 (stationary variance close to 1 at \gamma = 0.5, against 1.33 for ULA).
   • Deliverables: code, 4–6 page report ending with a one-page diagnostic checklist, 10-minute talk. Chapters: ch. 7, 8, 9.
   • Success criterion: over 50 replications, the 95 % MCMC intervals cover the exact mean at least 90 % of the time whenever the checklist says “converged”; the report shows one case where ESS looks fine and \hat R does not.
2. [Solo · ≈ 20 h] Rare events and the curse of dimension for importance sampling. How small a probability can crude Monte Carlo, exponential tilting and the cross-entropy method estimate, and how does a mismatched proposal break importance sampling as the dimension grows?
   • Implement: the three estimators for \mathbb{P}(Z > a) with a = 3 to 6, then for the loss of a portfolio of d correlated Gaussian factors exceeding a high quantile (closed form, for checking); the ESS of the proposal \mathcal{N}(\mu\mathbf{1}, I_d) against \mathcal{N}(0, I_d) for d up to 50.
   • Prove: \mathrm{ESS}/n = e^{-d\mu^2} (Prop. 5.2), and the second moment of the tilted estimator with g = \mathcal{N}(\mu, 1), which is e^{\mu^2}\,\Phi(-(a+\mu)).
   • Paper: de Boer, Kroese, Mannor and Rubinstein, “A tutorial on the cross-entropy method”, Annals of Operations Research 134, 2005.
   • Deliverables: code, 4-page report, 10-minute talk. Chapters: ch. 2, 4, 5.
   • Success criterion: every estimate comes with a CLT interval, and over 100 replications the 95 % intervals of the methods that work cover the exact value at least 90 % of the time. With n = 10^5 draws, the measured ESS follows e^{-d\mu^2} within 20 % while e^{-d\mu^2} \geq 0.15; beyond, the report shows that it becomes over-optimistic and explains why (the error on \sum_i w_i^2 involves \mathbb{E}_g[w^4] = e^{6d\mu^2}). The report states where each method fails.

Group projects

1. [Group of 3 · ≈ 25 h each] From ULA to score-based generation on toy data. Does the step-size bias of ULA (ch. 9) predict the quality of samples generated with a learned score, and what do the number of noise levels and of steps per level change?
   • Paper study: Song and Ermon, “Generative modeling by estimating gradients of the data distribution”, NeurIPS 2019 — a noise-conditional score network sampled by annealed Langevin (Prop. 9.6).
   • Implement: a small MLP score network trained by denoising score matching on 2-D data (Swiss roll and two moons from scikit-learn, a mixture of 8 Gaussians on a circle), then annealed Langevin; ULA, MALA and tamed ULA (paper 3 of ch. 10) on known targets, including the steep double well U(x) = x^4/4 - x^2/2 in dimension 10 to 100.
   • Question to answer: why can MALA not correct a learned score? (A score network gives no log-density, so the Metropolis-Hastings ratio of Def. 8.1 cannot be computed.)
   • Split of work: A trains the score networks; B writes the samplers and checks them on known targets; C owns the metrics (sliced Wasserstein distance to held-out data, mode weights), the seeds and the report.
   • Link: Multimodal Generative AI, ch. 3 (denoising score matching, Theorem 3.4, annealed Langevin).
   • Deliverables: repository, 6-page report, 15-minute talk. Chapters: ch. 7, 8, 9, 10.
   • Success criterion: on the Gaussian target the samplers reproduce Example 9.4 within error bars; ULA diverges on the double well for large steps while tamed ULA does not; on the 8-Gaussian mixture, annealed Langevin recovers every mode with weight 1/8 \pm 0.03 over 5 seeds; on the two-mode mixture with weights 1/5 and 4/5 of Song and Ermon, Langevin at a single noise level gets the weights wrong and annealed Langevin does not.
2. [Group of 3 · ≈ 30 h each] High-dimensional Bayesian logistic regression: ULA, MALA and Gibbs. At equal computing time, which sampler gives the most accurate posterior means, variances and predictive probabilities as the dimension grows? A paper study of Durmus and Moulines, “High-dimensional Bayesian inference via the Unadjusted Langevin Algorithm”, Bernoulli 25(4A), 2019 (paper 2 of ch. 10).
   • Data: German Credit (Statlog, OpenML id 31: 1,000 applicants, 20 attributes, one-hot encoded), and a synthetic design with d from 10 to 1,000 and known coefficients.
   • Implement: ULA and MALA with a Gaussian prior and a step scaled with d; the Pólya-Gamma Gibbs sampler of Polson, Scott and Windle (JASA, 2013), for example with the polyagamma package; a long MALA run as the reference posterior.
   • Measure: errors on posterior means and variances against the reference, ESS per second, held-out log-loss and calibration, and the ULA bias as \gamma varies (Theorem 9.5).
   • Split of work: A: ULA and step tuning; B: MALA and the reference chains; C: the Gibbs sampler and the predictive evaluation.
   • Deliverables: repository, 6-page report, 15-minute talk. Chapters: ch. 8, 9, 10.
   • Success criterion: every comparison is at equal wall-clock time, with error bars over 5 seeds; the report tests whether the ULA bias on posterior variances is linear in \gamma, as Example 9.4 predicts for a Gaussian, and says which sampler wins at d = 1{,}000 and why.
3. [Group of 3 · ≈ 30 h each] Metropolis-Hastings over graphs for molecule design. Can a Markov chain on molecular graphs, cooled by simulated annealing, find molecules with a high property score that stay valid, new and diverse?
   • Paper study: Xie, Shi, Zhou, Yang, Zhang, Yu and Li, “MARS: Markov molecular sampling for multi-objective drug discovery”, ICLR 2021.
   • Data: ZINC250k, about 250,000 drug-like molecules selected by Gómez-Bombarelli et al. (ACS Central Science, 2018), as SMILES; RDKit for validity, QED and logP.
   • Implement: an energy U(G) equal to minus a property (QED, or a GNN predictor trained on ZINC) plus validity and size penalties; graph edits (add or delete an atom or a bond, change an atom type) with their reverse probabilities, so that the ratio of Def. 8.1 is exact; temperature schedules (Prop. 8.6).
   • Check: detailed balance of the edit kernel (Theorem 8.2) on a space of small graphs small enough to enumerate.
   • Split of work: A: graph edits and the proposal ratio; B: the property predictor, with Graph ML ch. 9 (graph generation) as background; C: schedules, diagnostics (trace and ESS of the energy) and evaluation.
   • Deliverables: repository, 6-page report, 15-minute talk. Chapters: ch. 7, 8.
   • Success criterion: on the enumerable space, visit frequencies match \pi (chi-squared test at 5 %); on ZINC, three schedules are compared on the mean QED, validity, novelty and diversity of their top 100 molecules, with error bars over 5 seeds.

Written exam (proposal)

Three hours, one handwritten A4 sheet and a basic calculator allowed; four independent questions. Use \Phi(-3) \approx 1.35\times10^{-3}, \Phi(-6) \approx 9.87\times10^{-10} and e^9 \approx 8\,103.

1. Question (5 points). Two runs of random-walk Metropolis on the bimodal target of Example 8.7 report: run A, step 0.7, acceptance 0.71; run B, step 4, acceptance 0.32. A colleague prefers A “because it accepts more”. (a) What does each of the trace, the autocorrelation (or ESS) and \hat R detect? (b) Which run do you trust, and how do you check it? (c) Why is the 0.234 rule not the right guide here?

   Expected: (a) the trace shows mode trapping and drift; the autocorrelation gives \tau and the ESS, the variance inflation of the chain CLT (Prop. 8.4); \hat R compares overdispersed chains and reveals a mode that one chain has not visited. (b) B: it switches mode about every 10 iterations (ESS ≈ 20,000 against ≈ 200); check with at least 4 chains started in both modes, \hat R \approx 1 and mode weights near 1/2. (c) The rule comes from product targets in growing dimension; acceptance measures local moves, not crossings between modes.

2. Question (6 points). Let p = \mathbb{P}(Z > 3) with Z \sim \mathcal{N}(0,1). (a) How many crude Monte Carlo draws give a 10 % relative standard error? (b) With the proposal g = \mathcal{N}(3,1), write the weight w and show that \mathbb{E}_g\big[w(Y)^2\,\mathbf{1}_{Y>3}\big] = e^{9}\,\Phi(-6).

   3. Deduce the number of draws for a 10 % relative standard error, and the 95 % interval returned by n draws.

   Expected:

   • 1. The relative standard error is \sqrt{(1-p)/(np)}, so n \approx 100(1-p)/p \approx 7.4\times10^4.
   • 2. w(y) = e^{-3y + 9/2}; then w^2\,\varphi(y-3) = e^{9}\varphi(y+3), and integrating over y > 3 gives e^9\,\Phi(-6).
   • 3. e^9\Phi(-6) \approx 8.0\times10^{-6}; subtracting p^2 \approx 1.8\times10^{-6} leaves a per-draw standard deviation \approx 2.5\times10^{-3}, that is 1.84 times p, so n \approx 100 \times 1.84^2 \approx 340 draws — about 220 times fewer.
   • The interval is \hat p \pm 1.96\,\hat\sigma/\sqrt n, with \hat\sigma the empirical standard deviation of the weighted indicators.

3. Question (5 points). (a) On a finite space, prove that if \pi(x)P(x,y) = \pi(y)P(y,x) for all states, then \pi P = \pi. (b) Give an irreducible, aperiodic chain on three states whose invariant distribution is uniform but which is not reversible. (c) Prove that the Gibbs update of one coordinate is a Metropolis-Hastings step with acceptance probability 1.

   Expected:

   • 1. Sum the identity over x and use \sum_x P(y,x) = 1.
   • 2. The biased rotation P(x, x+1) = 2/3, P(x, x-1) = 1/3 (indices mod 3) is doubly stochastic, so the uniform law is invariant; it returns in 2 and 3 steps, so it is aperiodic; but \pi_1 P_{12} = 2/9 \neq 1/9 = \pi_2 P_{21}.
   • 3. As in Prop. 8.3: with the full conditional as proposal and y_{-i} = x_{-i}, the factorization \pi(x) = \pi(x_{-i})\,\pi(x_i \mid x_{-i}) makes the ratio equal to 1.

4. Question (4 points). Let \pi = \mathcal{N}(0, \sigma^2), so U(x) = x^2/(2\sigma^2). (a) Write the ULA recursion and find the steps for which it is stable. (b) Prove that its stationary variance is \sigma^2/(1 - \gamma/(2\sigma^2)). (c) In \mathbb{R}^d with independent coordinates of variances \sigma_1^2 \geq \dots \geq \sigma_d^2, explain why the number of iterations grows like \sigma_1^2/\sigma_d^2.

   Expected:

   • 1. The recursion is an AR(1), X_{k+1} = \big(1 - \gamma/\sigma^2\big)X_k + \sqrt{2\gamma}\,Z_{k+1}, stable if and only if 0 < \gamma < 2\sigma^2.
   • 2. The stationary variance v solves v = (1-\gamma/\sigma^2)^2 v + 2\gamma, which gives the formula (Example 9.4 when \sigma = 1).
   • 3. Stability and a small bias on the narrowest coordinate force \gamma \lesssim \sigma_d^2; the widest coordinate forgets its start like (1-\gamma/\sigma_1^2)^k, so it needs about \sigma_1^2/\gamma iterations — the ratio L/m of Theorem 9.5.

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

Grading rubric (proposal)

References

• A. Durmus, É. Moulines, Non-asymptotic convergence analysis for the Unadjusted Langevin Algorithm, Annals of Applied Probability 27(3), 2017 — Theorem 9.5 and chapter 10 [1].
• A. Durmus, É. Moulines, High-dimensional Bayesian inference via the Unadjusted Langevin Algorithm, Bernoulli 25(4A), 2019 — dimension, W_2 bounds, smoother potentials, empirical averages: chapter 10 [2].
• N. Brosse, A. Durmus, É. Moulines, S. Sabanis, The Tamed Unadjusted Langevin Algorithm, Stochastic Processes and their Applications 129(10), 2019 — chapter 10 [3].
• C. Robert, G. Casella, Monte Carlo Statistical Methods, 2nd ed., Springer, 2004 — the pedagogical reference for chapters 3, 5 and 8 of this module (rejection, IS, MH, Gibbs, diagnostics).
• D. Levin, Y. Peres, E. Wilmer, Markov Chains and Mixing Times, AMS, 2009 — free online version: complete proofs of Theorem 7.3 for finite chains, mixing times, spectrum (Example 7.4).
• A. Gelman, J. B. Carlin, H. S. Stern, D. B. Dunson, A. Vehtari, D. B. Rubin, Bayesian Data Analysis, 3rd ed., CRC Press, 2013 — book website and data: \hat R diagnostics, ESS, MCMC practice of ch. 8.
• J. R. Norris, Markov Chains, Cambridge University Press, 1997 — free pdf version: discrete-time chains (Norris’s ch. 1): invariant distributions, convergence to equilibrium, reversibility, ergodic theorem — Theorem 7.3 on countable spaces.
• P. Billingsley, Probability and Measure, 3rd ed., Wiley, 1995 — chapters 1–3: probability measures, Lebesgue measure and integral, product measure and Fubini (§18); chapter 5: characteristic functions and the CLT (§26–27, Thm 2.3); chapter 6: conditional expectation (§34).
• J. S. Liu, Monte Carlo Strategies in Scientific Computing, Springer, 2001 — importance sampling, weights and degeneracy, sequential methods — background for ch. 5.
• A. E. Gelfand, A. F. M. Smith, “Sampling-Based Approaches to Calculating Marginal Densities”, JASA 85(410), 1990 — the paper that launched the modern Gibbs sampler (Prop. 8.3).
• L. Tierney, “Markov Chains for Exploring Posterior Distributions”, Annals of Statistics 22(4), 1994 — convergence of MCMC chains; Metropolis-Hastings chains are Harris recurrent (Theorem 7.3).
• G. O. Roberts, J. S. Rosenthal, “General state space Markov chains and MCMC algorithms”, Probability Surveys 1, 2004 — Theorem 7.3 in general spaces, with proofs.
• G. O. Roberts, A. Gelman, W. R. Gilks, “Weak convergence and optimal scaling of random walk Metropolis algorithms”, Annals of Applied Probability 7(1), 1997 — the 0.234 rule of Definition 8.5 and its product-target assumption (Example 8.7).
• A. Gelman, G. O. Roberts, W. R. Gilks, “Efficient Metropolis jumping rules”, Bayesian Statistics 5, Oxford University Press, 1996 — the 0.44 acceptance rate for a one-dimensional Gaussian target (Definition 8.5).
• A. S. Dalalyan, Theoretical guarantees for approximate sampling from smooth and log-concave densities, Journal of the Royal Statistical Society B 79(3), 2017 — the first non-asymptotic bounds for ULA (chapter 10, paper 1).
• J. Ho, A. Jain, P. Abbeel, Denoising Diffusion Probabilistic Models, NeurIPS 2020 — Y. Song, S. Ermon, Generative Modeling by Estimating Gradients of the Data Distribution, NeurIPS 2019 (annealed Langevin, Prop. 9.6) — and Y. Song et al., Score-Based Generative Modeling through Stochastic Differential Equations, ICLR 2021 — the “generative models” continuation of Prop. 9.6 (see the Multimodal GenAI module, chapter 3).