Bayesian Machine Learning

Contents

• Introduction and overview
• 0. Review: probability tools for Bayesian inference
• 1. Decision theory
• 2. Fifty shades of Bayes: subjective and objective interpretations
• 3. Bayesian supervised learning
• 4. Bayesian unsupervised learning
• 5. Bayesian computation I: advanced Monte Carlo
• 6. Bayesian computation II: variational methods
• 7. Bayesian nonparametrics
• 8. Bayesian methods for deep learning
• 9. Complements for self-study
• 10. Research directions and open problems
• Exercises
• Projects and assessment
• References

Introduction and overview

Audience: a student with a good M1 in mathematics or computer science who knows maximum likelihood and has trained a classifier, and who wants models that say how much they do not know, and decisions that use it.

Empirical first: fit polynomials of degree 0 to 9 to 20 noisy points of a sine curve. The training error falls from 0.73 to 0.17 as the degree grows, so it cannot choose the degree. The Bayesian marginal likelihood of the same models, computed in closed form from the training data alone, peaks at degree 3 (ch. 3).

Five neural networks trained on 60 points agree to within 0.02 inside the data and spread to 0.24 far outside it, where all five miss the truth by 3 to 3.5 (ch. 8). Averaging models helps, but it is not yet a calibrated posterior.

This course treats Bayesian inference as a way of deciding under uncertainty. Choose the loss that matters, put a probability distribution on everything unknown, condition on the data, and act to minimize the posterior expected loss. The course then asks where priors come from, how to compute posteriors at the scale of machine learning, how to put priors on infinite-dimensional objects, and what survives of all this in deep learning.

It continues Introduction to Probabilistic Graphical Models and Deep Generative Models and Computational Statistics, and shares tools with Kernel Methods for Machine Learning and Learning Physical Phenomena: Metamodelling and Uncertainty Quantification.

What you will be able to do

• Formulate a learning task as a statistical decision problem and derive its Bayes rule for a given loss.
• Compare decision rules by frequentist risk, Bayes risk, admissibility and minimaxity, and explain the James–Stein phenomenon.
• Choose and defend a prior: exchangeability and elicitation, conjugate priors, Jeffreys and reference priors, empirical Bayes.
• Check the frequentist behaviour of a Bayesian procedure with the Bernstein–von Mises theorem and coverage simulations.
• Compute posterior and predictive distributions for Bayesian linear and logistic regression and for Gaussian processes, and select models by their marginal likelihood.
• Build and sample Bayesian latent-variable models: mixtures, probabilistic PCA, topic models.
• Implement and diagnose Metropolis–Hastings, Gibbs, Langevin, Hamiltonian and sequential Monte Carlo samplers.
• Derive mean-field, stochastic and black-box variational approximations, and say what they get wrong.
• Use the Dirichlet process and its mixtures, and state what posterior consistency and contraction rates mean.
• Apply Laplace, variational, dropout, stochastic-gradient MCMC and ensemble approximations to neural networks, and evaluate them with proper scoring rules and calibration.
• Read the abstract of a Bayesian machine learning paper, place it on the map of this course, and find a weak point worth a student seminar.

Official objectives and syllabus (2026–2027)

The course description on master-mva.com states that by the end of the course students should have a high-level view of the main approaches to making decisions under uncertainty, and be able to detect when being Bayesian helps and why. They should be able to design and run a Bayesian machine learning pipeline for standard supervised or unsupervised learning.

They should also have a global view of the current limitations of Bayesian approaches and of the research landscape, and be able to understand the abstract of most Bayesian machine learning papers.

The description lists six themes. It announces 8 × 3 hours of lectures and practical sessions, in Python and R, and 4 hours of student seminars for the evaluation. The Université Paris-Saclay course detail gives the volume (16 h of lectures, 12 h of labs) and the type of assessment, with no programme of its own.

Every official theme has a section in these notes. The extensions beyond the syllabus are the review of ch. 0 and two self-study chapters: the complements of ch. 9 (consistency and contraction rates, mixtures of finite mixtures, expectation propagation, accuracy of variational Bayes) and the research directions of ch. 10.

Where this course fits in the programme

Official prerequisites: an undergraduate course in probability, and a basic knowledge of Python and R. The description recommends the semester-1 course of P. Latouche and N. Chopin on probabilistic graphical models or that of S. Allassonière on computational statistics. Ch. 0 reviews the probability used throughout.

• Refreshers — Probability and Statistics (conditional laws, Gaussian vectors, Markov chains), Linear Algebra, Analysis and Optimization (Schur complements, Hessians) and Scientific Python, NumPy and PyTorch (the code).
• Semester 1 courses that prepare it — Introduction to Probabilistic Graphical Models and Deep Generative Models and Computational Statistics, the two recommended courses, give the first encounter with EM, mixtures, MCMC and variational methods. Stochastic Calculus in Machine Learning treats the Langevin diffusions behind ch. 5 and ch. 8. Deep Learning and Theoretical Foundations of Deep Learning prepare ch. 8.
• Related courses of semester 2 — Kernel Methods for Machine Learning (the Gaussian-process posterior mean of ch. 3 is kernel ridge regression), Learning Physical Phenomena (uncertainty quantification), Inverse Problems and Imaging: Statistical and Stochastic Approaches and Diffusion and Multiscale Models in Artificial Intelligence (posterior sampling with Langevin dynamics), Sequential Learning (online decisions under uncertainty) and Responsible Machine Learning.
• Afterwards — the research internship: the paper-study format of the evaluation is a rehearsal for it.

Course outline

• 0. Review — Bayes’ formula, Gaussian conditioning, exponential families and conjugacy, the Kullback–Leibler divergence, Markov chains.
• 1. Decision theory — losses and risks, Bayes rules, admissibility and minimaxity, complete classes, Stein’s paradox, cost-sensitive classification.
• 2. Fifty shades of Bayes — subjective foundations (coherence, exchangeability), conjugate and objective priors, empirical Bayes, Bernstein–von Mises, the likelihood principle.
• 3. Bayesian supervised learning — linear regression, marginal likelihood, the Laplace approximation, logistic regression, Gaussian processes.
• 4. Bayesian unsupervised learning — mixtures and the Gibbs sampler, label switching, probabilistic PCA, topic models.
• 5. Advanced Monte Carlo — importance sampling, Metropolis–Hastings, Langevin and Hamiltonian samplers, diagnostics, sequential Monte Carlo.
• 6. Variational methods — the ELBO, mean field and coordinate ascent, stochastic and black-box VI, what the reverse KL misses.
• 7. Bayesian nonparametrics — the Dirichlet process, its urn and stick-breaking representations, DP mixtures and their collapsed Gibbs sampler, Pitman–Yor and hierarchical extensions.
• 8. Bayesian deep learning — priors on weights, the Gaussian-process limit, Laplace, VI, dropout, SG-MCMC and ensembles, calibration, cold posteriors.
• 9. Complements for self-study — posterior consistency and contraction rates, the number of clusters, the Pitman–Yor urn, expectation propagation, the accuracy of variational Bayes.
• 10. Research directions and open problems — a self-study tour of current work, and how to read a paper.

Session plan

28 h in class = 9 sessions (8 × 3 h + 1 × 4 h). The volume of 28 h comes from the Université Paris-Saclay programme table (16 h of lectures and 12 h of labs). master-mva.com describes the same total as 8 × 3 h of lectures and practicals plus 4 h of student seminars for the evaluation. The content of each session and its split between lecture and practical are a proposal.

Level 1 = foundations, 2 = core results and algorithms, 3 = advanced and research-level material. The level never goes down from one session to the next.

Self-study, outside the 28 h (about 6 h in total): ch. 9 Complements for self-study (≈ 3 h) · ch. 10 Research directions and open problems (≈ 3 h).

How to study this module

• Read ch. 0 before session 1 if Gaussian conditioning or Markov chains feel distant; every later chapter uses them without warning.
• After each session, do the exercises listed in its row, then run the code of the chapter and change one ingredient (the prior, the step size, the number of samples) to see its effect.
• Write each sampler yourself once (Gibbs, Metropolis–Hastings, HMC, collapsed Gibbs for the DP mixture) before using a probabilistic programming language; each one is a few lines of NumPy.
• The code uses NumPy, SciPy and PyTorch (CPU) with synthetic data; every block runs in under 30 seconds on a laptop.
• Read ch. 9 after session 7 if you want the theory behind nonparametric posteriors and variational approximations.
• Choose the paper for the seminar after session 5, and use ch. 10 to place it in the landscape.

0. Review: probability tools for Bayesian inference

Bayesian machine learning manipulates joint, marginal and conditional distributions, Gaussian vectors, exponential families, the Kullback–Leibler divergence and Markov chains. This chapter states the M1-level facts used later. The Probability and Statistics refresher covers the same ground at a slower pace.

Bayes’ formula for densities

Two remarks are used constantly. First, the posterior is known up to a constant as soon as the product of likelihood and prior is, and the constant p(x) is usually the hard part. Second, for independent observations the posterior can be updated one observation at a time: yesterday’s posterior is today’s prior.

Gaussian vectors

Exponential families and conjugacy

The Kullback–Leibler divergence

The divergence is not symmetric. Variational inference (ch. 6) minimizes \mathrm{KL}(q \,\|\, p) over an approximation q, and expectation propagation works with the other direction; the two give different answers.

Markov chains

This theorem is the reason Markov chain Monte Carlo works (ch. 4 and 5): build a chain whose invariant distribution is the posterior, run it, and average along the path.

1. Decision theory

Learning is acting: predicting a label, choosing a dose, flagging a transaction. Decision theory compares ways of acting when the state of the world is unknown, and shows that the Bayes rule answers a precise question.

Losses, risks and Bayes rules

Estimation takes \mathcal{A} = \Theta and, for instance, the squared loss. Testing takes \mathcal{A} = \{0, 1\} and the 0–1 loss. Probabilistic prediction takes for \mathcal{A} a set of distributions and the log loss, the setting of ch. 8.

The risk is a function of \theta, and two rules are rarely ordered for every \theta. There are two classical ways to summarize a risk function by a number: average it over \theta (Bayes), or take its worst case (minimax).

The theorem is the operational content of Bayesian decision theory. There is no need to optimize over functions: after seeing x, minimize a function of the action. The average-over-data question and the given-the-data question have the same answer.

For a continuous parameter, the MAP is in general not the Bayes estimator of a fixed loss. Under regularity conditions, it is only a limit of Bayes estimators for the losses \mathbb{1}\{|\theta - a| > \varepsilon\} as \varepsilon \to 0. It also depends on the parameterization, through the Jacobian of the change of variables (Exercise 13). It is best seen as a cheap summary of the posterior.

Admissibility and minimaxity

The code computes the exact risk curves of three estimators. It reports their worst-case risk and their average risk under the uniform prior:

import numpy as np
from scipy import stats

n = 16
theta = np.linspace(0.001, 0.999, 999)
k = np.arange(n + 1)
pmf = stats.binom.pmf(k[None, :], n, theta[:, None])     # P(K = k | theta)

def risk(est):                                             # squared-error risk R(theta, delta)
    return (pmf * (est[None, :] - theta[:, None]) ** 2).sum(axis=1)

rules = {
    "MLE k/n": k / n,
    "Bayes, uniform prior": (k + 1) / (n + 2),
    "minimax, Beta(2,2) prior": (k + 2) / (n + 4),          # a = sqrt(16)/2 = 2
}
for name, est in rules.items():
    r = risk(est)
    print(f"{name:26s} max risk {r.max():.5f}   Bayes risk (uniform) {r.mean():.5f}")

The maximum-likelihood estimator has worst-case risk 0.0156 and average risk 0.0104. The uniform-prior Bayes rule has the smallest average, 0.0093, as it must, and a worst case of 0.0123. The minimax rule has the smallest worst case, 0.0100, flat in \theta. No estimator wins for every \theta.

Complete classes and Stein’s paradox

This is one precise sense in which being Bayesian helps: any sensible rule, in the sense of admissibility, is Bayes for some prior, so restricting attention to Bayes rules loses nothing. The next result shows that the most natural non-Bayesian estimator is not admissible.

At \theta = 0, \|X\|^2 follows a \chi^2_d law with \mathbb{E}[1/\chi^2_d] = 1/(d-2), so the risk drops from d to 2. For d = 10, that is from 10 to 2. Shrinking all coordinates towards a common point helps, although they are independent: ch. 2 explains why, as empirical Bayes (Example 2.10).

Decisions in classification and the foundations

The threshold is applied to the predictive probability. If that probability is miscalibrated, the decision is wrong even with the right costs; this is why ch. 3 and ch. 8 care about predictive distributions rather than point predictions. An option to abstain at cost c_r adds a third action, chosen when both expected losses exceed c_r.

2. Fifty shades of Bayes: subjective and objective interpretations

Theorem 1.3 needs a prior. Where does it come from, and what does a posterior probability mean? Two families of answers coexist. For the subjective view, probability is a coherent degree of belief. For the objective view, the prior is a convention derived from the model. Two pragmatic positions complete the picture: estimate the prior from the data (empirical Bayes), or judge Bayesian procedures by their frequentist properties.

Subjective Bayes: coherence and exchangeability

Coherence gives probabilities, not a specific prior. De Finetti’s second theorem explains where parameters and priors come from in the first place.

The interpretation matters. If your beliefs about a sequence do not depend on the order of the observations, you behave as if there were a parameter \theta, with a prior \mu, and independent data given \theta. The parameter and the prior are consequences of a judgement on observable quantities.

Conjugate priors as an elicitation tool

Proposition 0.5 reads a conjugate prior as pseudo-data. An expert can express a prior mean and an equivalent sample size, “as sure as if I had seen m observations”, and the two numbers fix the prior.

A responsible analysis reports how the conclusions move when the prior moves within a plausible family. This sensitivity analysis is the subjective answer to the charge of arbitrariness.

Objective Bayes: Jeffreys and reference priors

The point of the proposition: “ignorance” cannot be expressed by a flat prior, since a flat prior on \theta is not flat on \theta^2. Jeffreys’ rule gives the same answer whatever parameterization one starts from.

The last three priors are improper: they do not integrate to 1. Such a prior is acceptable only if the posterior is proper, that is if \int p(x \mid \theta)\, \pi(\theta)\, d\theta is finite. Improper priors leave the evidence defined only up to an arbitrary constant, so they cannot be used to compare models (ch. 3).

In several dimensions the Jeffreys prior can behave badly. Reference priors (Bernardo, 1979; Berger, Bernardo and Sun, 2009) fix this: they maximize the expected Kullback–Leibler divergence from prior to posterior, in the limit of many replications of the experiment. In regular one-dimensional models they coincide with the Jeffreys prior.

Empirical Bayes

Frequentist properties of Bayesian procedures

Three consequences follow. In fixed dimension, the prior is forgotten at rate 1/\sqrt{n}. A 95 % credible set is asymptotically a 95 % confidence set. And the theorem fails exactly where Bayesian methods are most interesting: in high dimension, in nonparametric models (ch. 7), and under misspecification, where the posterior concentrates at the best approximation within the model but with the wrong spread (Kleijn and van der Vaart, 2012).

The code computes the exact frequentist coverage of three 95 % intervals for a proportion with n = 30: the Wald interval, and the equal-tailed credible intervals under the Jeffreys and the uniform priors.

import numpy as np
from scipy import stats

n = 30
k = np.arange(n + 1)
z = stats.norm.ppf(0.975)

def coverage(lo, hi, theta):
    inside = (lo <= theta) & (theta <= hi)
    return stats.binom.pmf(k, n, theta)[inside].sum()

phat = k / n
wald = (phat - z * np.sqrt(phat * (1 - phat) / n), phat + z * np.sqrt(phat * (1 - phat) / n))
jeff = (stats.beta.ppf(0.025, k + 0.5, n - k + 0.5), stats.beta.ppf(0.975, k + 0.5, n - k + 0.5))
unif = (stats.beta.ppf(0.025, k + 1, n - k + 1), stats.beta.ppf(0.975, k + 1, n - k + 1))

thetas = np.linspace(0.01, 0.99, 197)
for name, (lo, hi) in {"Wald": wald, "Jeffreys": jeff, "uniform": unif}.items():
    c = np.array([coverage(lo, hi, t) for t in thetas])
    print(f"{name:9s} mean coverage {c.mean():.3f}   min coverage {c.min():.3f}"
          f"   at theta=0.1: {coverage(lo, hi, 0.1):.3f}")

The Wald interval has mean coverage 0.887 and drops to 0.26 near the boundary. The Jeffreys interval has mean coverage 0.950 and minimum 0.889; the uniform-prior interval has 0.952 and 0.879. At \theta = 0.1 the three coverages are 0.809, 0.932 and 0.974. Brown, Cai and DasGupta (2001) recommend the Jeffreys interval on these grounds: an objective prior gives a good frequentist procedure.

The likelihood principle

Bayesian inference depends on the data only through the likelihood function, up to a constant factor. Frequentist inference also depends on the experiments that could have happened but did not. Birnbaum (1962) showed that the sufficiency and conditionality principles together imply this likelihood principle.

from scipy import stats
# 9 heads and 3 tails: P(at least 9 heads) under each design, for a fair coin
print("design A, binomial n = 12:        ", stats.binom.sf(8, 12, 0.5))
print("design B, stop at the third tail: ", stats.nbinom.sf(8, 3, 0.5))

3. Bayesian supervised learning

Supervised learning the Bayesian way: put a prior on the function, through its parameters or directly, condition on the training pairs, and predict by averaging over the posterior. Linear-Gaussian models give closed forms; other models need approximations, the first of which is Laplace’s.

Bayesian linear regression

Take features \phi(x) \in \mathbb{R}^p, the design matrix \Phi with rows \phi(x_i)^\top, Gaussian noise of precision \beta and the prior w \sim \mathcal{N}(0, \alpha^{-1} I): y_i = w^\top \phi(x_i) + \varepsilon_i, \qquad \varepsilon_i \sim \mathcal{N}(0, \beta^{-1}) .

The predictive variance has two parts: the noise \beta^{-1}, which no amount of data removes (aleatoric uncertainty), and the parameter uncertainty \phi^{\star\top} S_N \phi^\star, which shrinks as data accumulate near x^\star (epistemic uncertainty).

The marginal likelihood and model selection

The evidence implements Occam’s razor. The first term rewards fit. The log-determinant penalizes models that spread their prior predictive mass over many possible data sets: since p(y \mid \mathcal{M}) integrates to 1 over y, a flexible model cannot give high density to everything.

To compare two models, Bayes’ formula multiplies the prior odds by the Bayes factor p(y \mid \mathcal{M}_1) / p(y \mid \mathcal{M}_2) (Kass and Raftery, 1995). The hyperparameters \alpha, \beta can be fitted by maximizing the evidence (Definition 2.9), as MacKay (1992) proposed.

The code fits polynomials of degree 0 to 9 to 20 noisy points of \sin(\pi x) and prints the training error of least squares and the log evidence:

import numpy as np

rng = np.random.default_rng(1)
n, alpha, beta = 20, 1.0, 1 / 0.2**2          # prior precision, noise precision (sigma = 0.2)
x = np.sort(rng.uniform(-1, 1, n))
y = np.sin(np.pi * x) + 0.2 * rng.normal(size=n)

def log_evidence(Phi):
    """log p(y) for y = Phi w + noise, w ~ N(0, I/alpha), noise ~ N(0, I/beta)."""
    C = np.eye(n) / beta + Phi @ Phi.T / alpha   # marginal covariance of y
    _, logdet = np.linalg.slogdet(C)
    return -0.5 * (y @ np.linalg.solve(C, y) + logdet + n * np.log(2 * np.pi))

for deg in range(10):
    Phi = np.vander(x, deg + 1, increasing=True)
    w_ml = np.linalg.lstsq(Phi, y, rcond=None)[0]
    rmse = np.sqrt(np.mean((Phi @ w_ml - y) ** 2))
    print(f"degree {deg}: training RMSE {rmse:.3f}   log evidence {log_evidence(Phi):7.2f}")

The training error decreases from 0.726 to 0.165. The log evidence jumps from −55.3 at degree 2 to −19.2 at degree 3, a Bayes factor of e^{36}, then slowly decreases to −20.3 at degree 9. The evidence picks degree 3 without any validation set, although its preference over degree 9 is weak (a factor of about 3).

A warning, known as Lindley’s paradox (Lindley, 1957): as the prior on the extra coefficients becomes vaguer, \alpha \to 0, the log-determinant grows without bound and the evidence of the larger model goes to 0. Evidence-based comparisons need proper, carefully chosen priors.

The Laplace approximation

Both come from the second-order Taylor expansion of \ell at \hat\theta and a Gaussian integral. For regular models with n observations, the relative error on the evidence is O(1/n) (Tierney and Kadane, 1986).

Bayesian logistic regression

With labels y_i \in \{0, 1\}, P(y_i = 1 \mid w) = \sigma(w^\top x_i) and the prior \mathcal{N}(0, \alpha^{-1} I), the posterior has no closed form. The log posterior is concave, so the MAP is found by Newton’s method, and the Hessian at the mode is H = \alpha I + \sum_{i=1}^n \sigma_i (1 - \sigma_i)\, x_i x_i^\top, \qquad \sigma_i = \sigma(\hat w^\top x_i).

Gaussian processes

Bayesian linear regression is a Gaussian process with k(x, x') = \alpha^{-1} \phi(x)^\top \phi(x'). Kernels with infinitely many features, such as the squared-exponential kernel, give priors on infinite-dimensional spaces of functions: a first Bayesian nonparametric model (ch. 7).

The posterior mean is the kernel ridge regression predictor with n\lambda = \sigma^2, studied in Kernel Methods for Machine Learning. The posterior variance does not depend on y. Exact inference costs O(n^3), through a Cholesky factorization. The code fits a GP to 30 points with a gap in (-1, 1) and compares length-scales by marginal likelihood:

import numpy as np

rng = np.random.default_rng(0)
x = np.concatenate([rng.uniform(-3, -1, 15), rng.uniform(1, 3, 15)])
y = np.sin(2 * x) + 0.1 * rng.normal(size=x.size)
noise = 0.1**2

def k(a, b, ell):
    return np.exp(-(a[:, None] - b[None, :]) ** 2 / (2 * ell**2))

def gp(ell, xs):
    K = k(x, x, ell) + noise * np.eye(x.size)
    L = np.linalg.cholesky(K)
    a = np.linalg.solve(L.T, np.linalg.solve(L, y))
    lml = -0.5 * y @ a - np.log(np.diag(L)).sum() - 0.5 * x.size * np.log(2 * np.pi)
    Ks = k(xs, x, ell)
    v = np.linalg.solve(L, Ks.T)
    sd = np.sqrt(1 - (v**2).sum(0) + noise)          # predictive sd of a new observation
    return lml, Ks @ a, sd

xs = np.array([-2.0, 0.0, 6.0])
for ell in [0.1, 0.3, 0.6, 1.0, 3.0]:
    lml, mean, sd = gp(ell, xs)
    print(f"ell={ell}: log marginal lik. {lml:7.2f}  sd at x=-2, 0, 6: {np.round(sd, 2)}")

The marginal likelihood is largest at \ell = 0.6 (8.41, against 7.24 at \ell = 1 and −271.6 at \ell = 3). With \ell = 0.6, the predictive standard deviation is 0.13 among the data, close to the noise level 0.1, then 0.88 in the gap and 1.00 far away, where the prior takes over. The over-smooth \ell = 3 gives 0.11 in the gap: confidently wrong.

Before trusting a model, check it: simulate replicated data sets from the posterior predictive and compare a statistic of interest with its observed value (Box, 1980; Gelman, Meng and Stern, 1996).

4. Bayesian unsupervised learning

Unsupervised learning looks for structure that is not observed: clusters, low-dimensional factors, topics. The Bayesian version writes a generative model with latent variables and computes their posterior. Conjugacy makes the conditional distributions simple, and the Gibbs sampler exploits it.

Mixtures and the Gibbs sampler

The Probabilistic Graphical Models course fits this model by EM, which returns one maximizer of the likelihood. The Bayesian treatment keeps the whole posterior, and the simplest way to explore it is to update one block of variables at a time.

Invariance is not enough: the chain must also be irreducible, which holds here because all conditionals have full support. The ergodic theorem (Theorem 0.10) then justifies the averages.

The code runs this sampler on 300 points from three components with means -4, 0, 3 and weights 0.3, 0.5, 0.2:

import numpy as np

rng = np.random.default_rng(0)
true_mu, true_w = np.array([-4.0, 0.0, 3.0]), np.array([0.3, 0.5, 0.2])
z_true = rng.choice(3, size=300, p=true_w)
x = rng.normal(true_mu[z_true], 1.0)                 # known unit variance

K, s0, a0 = 3, 10.0, 1.0                              # mu_k ~ N(0, s0^2), w ~ Dir(a0, ..., a0)
mu, w = rng.normal(0, 1, K), np.full(K, 1 / K)
draws = []
for it in range(2000):
    # 1. allocations: p(z_i = k | ...) proportional to w_k N(x_i; mu_k, 1)
    logp = np.log(w) - 0.5 * (x[:, None] - mu[None, :]) ** 2
    p = np.exp(logp - logp.max(1, keepdims=True))
    p /= p.sum(1, keepdims=True)
    z = (p.cumsum(1) > rng.uniform(size=(x.size, 1))).argmax(1)
    n_k = np.bincount(z, minlength=K)
    # 2. weights: Dirichlet(a0 + n_k)
    w = rng.dirichlet(a0 + n_k)
    # 3. means: Gaussian, precision 1/s0^2 + n_k, mean (sum of x in k) / precision
    prec = 1 / s0**2 + n_k
    mu = rng.normal(np.bincount(z, weights=x, minlength=K) / prec, 1 / np.sqrt(prec))
    if it >= 500:
        order = np.argsort(mu)                         # relabel by sorted means
        draws.append(np.concatenate([mu[order], w[order]]))
d = np.array(draws)
print("posterior mean of (mu_1, mu_2, mu_3):", d[:, :3].mean(0).round(2))
print("posterior sd   of (mu_1, mu_2, mu_3):", d[:, :3].std(0).round(2))
print("posterior mean of (w_1, w_2, w_3):  ", d[:, 3:].mean(0).round(2))

The posterior means of the component means are (-4.08, 0.02, 2.88), with standard deviations 0.12 to 0.16. The weights come out as (0.25, 0.48, 0.27): they match the proportions actually drawn in this sample, (0.253, 0.483, 0.263), rather than the generating weights.

Label switching

So a sampler that explores the whole posterior gives useless component-wise averages: each average is the mean of all components. The code sorts the means at each iteration, an identifiability constraint that is harmless here because the components are well separated. In general one relabels the draws (Stephens, 2000) or reports label-free quantities: the predictive density, or the co-clustering probabilities P(z_i = z_j \mid x).

The Gibbs sampler above never switches labels, because crossing between the symmetric modes requires moving many allocations at once. This is a failure to mix that happens to be convenient (Celeux, Hurn and Robert, 2000).

The number of components K can be chosen by the evidence, estimated with the methods of ch. 5. Alternatively, take K large with a small a_0, below half the dimension of a component’s parameter (here a_0 < 1/2): Rousseau and Mengersen (2011) show that the posterior then empties the superfluous components. Ch. 7 removes the choice altogether, with infinitely many components.

Probabilistic and Bayesian PCA

Bayesian PCA (Bishop, 1999) puts a prior \mathcal{N}(0, \alpha_j^{-1} I_p) on each column of W and fits the precisions \alpha_j by evidence maximization. Columns not supported by the data get \alpha_j \to \infty and are switched off. This automatic relevance determination selects the latent dimension.

Topic models

5. Bayesian computation I: advanced Monte Carlo

Posterior expectations are integrals in high dimension against a density known up to a constant. Monte Carlo methods replace them by averages over samples. This chapter goes beyond the Gibbs sampler: importance sampling and its curse of dimensionality, Metropolis–Hastings, gradient-based samplers, diagnostics, and sequential Monte Carlo, which also estimates the evidence.

Throughout, the target is \pi = \gamma / Z, where \gamma(\theta) = p(x \mid \theta)\, p(\theta) can be evaluated and Z = p(x) is unknown.

Importance sampling

The effective sample size (\sum_i w_i)^2 / \sum_i w_i^2, between 1 and N, measures how many draws effectively carry the estimate.

Independent proposals fail in high dimension. Markov chains move locally instead.

Metropolis–Hastings

For the random-walk proposal \theta' = \theta + h\xi, with \xi standard Gaussian, on product targets in dimension d, the optimal step is of order d^{-1/2} and the optimal acceptance rate is 0.234 (Roberts, Gelman and Gilks, 1997). The chain then needs of order d iterations to move across the target.

Langevin and Hamiltonian samplers

The Langevin diffusion d\theta_t = \frac12 \nabla \log \pi(\theta_t)\, dt + dB_t has \pi as invariant distribution: for smooth compactly supported f, an integration by parts gives \int \Big(\tfrac12 \nabla\log\pi \cdot \nabla f + \tfrac12 \Delta f\Big)\, \pi \;=\; \tfrac12 \int \nabla\pi \cdot \nabla f - \tfrac12 \int \nabla \pi \cdot \nabla f = 0 . Its Euler discretization with step h^2 is the unadjusted Langevin algorithm (ULA): \theta_{k+1} = \theta_k + \frac{h^2}{2} \nabla \log\gamma(\theta_k) + h\, \xi_k .

ULA is biased: for a fixed step, its invariant law is not \pi. Using the same move as a Metropolis–Hastings proposal removes the bias: this is MALA. Its optimal step is of order d^{-1/6}, with acceptance rate 0.574 (Roberts and Rosenthal, 1998). For strongly log-concave targets, ULA has explicit non-asymptotic error bounds (Dalalyan, 2017; Durmus and Moulines, 2017). The diffusions themselves are studied in Stochastic Calculus in Machine Learning.

Hamiltonian Monte Carlo (HMC) adds a momentum p \sim \mathcal{N}(0, I) and the energy H(\theta, p) = -\log\gamma(\theta) + \frac12 \|p\|^2. The Hamiltonian flow preserves H and volume, so it can travel far along a level set of the joint density. It is integrated by the leapfrog scheme with step \varepsilon: p_{1/2} = p + \tfrac{\varepsilon}{2}\nabla\log\gamma(\theta), \qquad \theta' = \theta + \varepsilon\, p_{1/2}, \qquad p' = p_{1/2} + \tfrac{\varepsilon}{2}\nabla\log\gamma(\theta').

The leapfrog energy error is O(\varepsilon^2) over a trajectory of fixed length, so long moves are accepted with high probability; the optimal step scales as d^{-1/4} (Beskos et al., 2013). The No-U-Turn sampler (Hoffman and Gelman, 2014), used by Stan, chooses the trajectory length adaptively.

The code compares the three samplers on a 50-dimensional Gaussian whose standard deviations range from 0.1 to 1, with steps tuned to the usual acceptance rates. It reports the effective sample size of the widest coordinate.

import numpy as np

rng = np.random.default_rng(0)
d = 50
sd = np.linspace(0.1, 1.0, d)                  # target N(0, diag(sd^2)), condition number 100
logpi = lambda x: -0.5 * np.sum((x / sd) ** 2)
grad = lambda x: -x / sd**2

def ess(chain):                                # sum autocorrelations until the first negative one
    c = chain - chain.mean()
    f = np.fft.rfft(c, 2 * c.size)
    acf = np.fft.irfft(f * np.conj(f))[:c.size] / (c.var() * c.size)
    tau = 1 + 2 * np.sum(acf[1:np.argmax(acf < 0)])
    return c.size / tau

def rwmh(h, n):
    x, out, acc = np.zeros(d), [], 0
    for _ in range(n):
        y = x + h * rng.normal(size=d)
        if np.log(rng.uniform()) < logpi(y) - logpi(x):
            x, acc = y, acc + 1
        out.append(x[-1])
    return np.array(out), acc / n

def mala(h, n):
    x, out, acc = np.zeros(d), [], 0
    logq = lambda b, a: -np.sum((b - a - 0.5 * h**2 * grad(a)) ** 2) / (2 * h**2)
    for _ in range(n):
        y = x + 0.5 * h**2 * grad(x) + h * rng.normal(size=d)
        if np.log(rng.uniform()) < logpi(y) + logq(x, y) - logpi(x) - logq(y, x):
            x, acc = y, acc + 1
        out.append(x[-1])
    return np.array(out), acc / n

def hmc(eps, L, n):
    x, out, acc = np.zeros(d), [], 0
    for _ in range(n):
        p0 = rng.normal(size=d)
        q, p = x.copy(), p0 + 0.5 * eps * grad(x)          # leapfrog integrator
        for i in range(L):
            q = q + eps * p
            p = p + (eps if i < L - 1 else 0.5 * eps) * grad(q)
        if np.log(rng.uniform()) < logpi(q) - 0.5 * p @ p - logpi(x) + 0.5 * p0 @ p0:
            x, acc = q, acc + 1
        out.append(x[-1])
    return np.array(out), acc / n

n = 10000
for name, (ch, acc), grads in [("RWMH h=0.1", rwmh(0.1, n), 0),
                               ("MALA h=0.16", mala(0.16, n), 1),
                               ("HMC eps=0.08, L=20", hmc(0.08, 20, n), 20)]:
    e = ess(ch[n // 10:])
    print(f"{name:20s} acceptance {acc:.2f}   ESS of x_50 per 1000 iterations {1000 * e / (0.9 * n):6.1f}"
          f"   gradients per iteration {grads}")

The acceptance rates are 0.27, 0.56 and 0.93. Per 1000 iterations, the effective sample sizes are about 1 for the random walk, 3 for MALA and 900 for HMC. Even counted per gradient evaluation, HMC gives about 45 effective draws per 1000 gradients, against 2.6 for MALA. The small scales force small steps, and only HMC turns them into long moves.

This comparison flatters HMC a little. Its trajectory length \varepsilon L = 1.6 is close to \pi/2, a quarter of the oscillation period of the unit-scale coordinate, where the exact Hamiltonian flow returns an independent draw. With L between 10 and 40, HMC gives about 17 to 45 effective draws per 1000 gradients. The tuning matters, but HMC stays far ahead of MALA; the No-U-Turn sampler chooses L automatically.

Diagnostics

The \hat R statistic (Gelman and Rubin, 1992; Vehtari et al., 2021) runs several chains from dispersed starting points and compares between-chain and within-chain variances; values above 1.01 signal trouble. No diagnostic proves convergence: they can only reveal its failure, as the label-switching example of ch. 4 shows.

Sequential Monte Carlo and the evidence

The product \hat r_1 \cdots \hat r_T is an unbiased estimator of the evidence Z = p(x), and it is consistent as N \to \infty (Del Moral, Doucet and Jasra, 2006). Temperatures are usually chosen adaptively, so that each step keeps a fixed fraction of effective particles. This gives the evidence that ch. 3 and ch. 4 needed for model choice. Annealed importance sampling (Neal, 2001) is the same scheme without resampling; Chopin and Papaspiliopoulos (2020) give a full treatment.

When n is large, each evaluation of \gamma costs O(n). Bardenet, Doucet and Holmes (2017) review subsampling MCMC: naive subsampling destroys the invariant distribution, and exact methods need good control variates. Stochastic-gradient Langevin dynamics, in ch. 8, trades exactness for speed.

6. Bayesian computation II: variational methods

Monte Carlo is exact in the limit but can be slow. Variational inference (VI) turns inference into optimization: choose a family \mathcal{Q} of tractable distributions and find its member closest to the posterior. VI scales to large data and models, at the price of a bias that this chapter makes precise.

The evidence lower bound

Three consequences follow. The ELBO is a lower bound on the log evidence (Proposition 0.8). Maximizing it over q is the same as minimizing the divergence to the posterior. And neither requires p(x).

Mean field and coordinate ascent

The mean-field family contains the product densities q(\theta) = \prod_j q_j(\theta_j) over blocks of variables.

Cycling through the blocks never decreases the ELBO, and converges to a local optimum. The update has the same shape as a Gibbs conditional, with random neighbours replaced by expectations.

With unit variances and correlation \rho = 0.9, the mean-field variance is 1 - \rho^2 = 0.19: a standard deviation of 0.44 instead of 1. The reverse KL penalizes mass of q where the posterior has little, so it prefers to sit inside the posterior. Minimizing \mathrm{KL}(p \,\|\, q) over the same family gives the exact marginals instead (Exercise 10).

In conditionally conjugate exponential-family models, the CAVI update of a global parameter adds the expected sufficient statistics of all data points to the prior’s natural parameter. Stochastic VI (Hoffman, Blei, Wang and Paisley, 2013) estimates this sum from a minibatch, scaled by n over the batch size, and moves part of the way with a decreasing step. The resulting step is a stochastic natural-gradient step on the ELBO, which scales VI to millions of documents.

Black-box and reparameterized VI

For non-conjugate models, choose a parametric family q_\phi and maximize \mathrm{ELBO}(\phi) by stochastic gradient ascent. Two unbiased gradient estimators exist.

The reparameterization estimator usually has a much smaller variance, but it needs continuous parameters and a differentiable model. For q = \mathcal{N}(m, \mathrm{diag}(e^{2s})), take T(\varepsilon) = m + e^{s} \odot \varepsilon; the entropy of q is \sum_j s_j plus a constant.

The code fits this mean-field Gaussian to the posterior of a two-parameter Bayesian logistic regression with strongly correlated features, and compares it with the exact posterior computed on a grid:

import numpy as np
import torch

torch.set_num_threads(1)
torch.manual_seed(0)
rng = np.random.default_rng(0)
n = 40
X = rng.normal(size=(n, 2))
X[:, 1] = 0.9 * X[:, 0] + 0.45 * X[:, 1]              # correlated features
y = (rng.uniform(size=n) < 1 / (1 + np.exp(-X @ np.array([1.5, -1.0])))).astype(float)
Xt, yt = torch.tensor(X), torch.tensor(y)

def log_joint(w):                                       # w: (S, 2); prior N(0, 4 I)
    logits = w @ Xt.T
    loglik = (yt * logits - torch.nn.functional.softplus(logits)).sum(1)
    return loglik - (w**2).sum(1) / 8

# exact posterior moments on a grid (feasible in two dimensions)
g = torch.linspace(-8, 8, 401, dtype=torch.float64)
W = torch.cartesian_prod(g, g)
p = torch.softmax(log_joint(W), 0)
m_exact = (p[:, None] * W).sum(0)
sd_exact = ((p[:, None] * (W - m_exact) ** 2).sum(0)).sqrt()

# mean-field Gaussian q(w) = N(m, diag(exp(2 s))), reparameterization gradient of the ELBO
m = torch.zeros(2, dtype=torch.float64, requires_grad=True)
s = torch.zeros(2, dtype=torch.float64, requires_grad=True)
opt = torch.optim.Adam([m, s], lr=0.05)
for step in range(3000):
    eps = torch.randn(64, 2, dtype=torch.float64)
    w = m + eps * s.exp()
    elbo = log_joint(w).mean() + s.sum()               # entropy of q up to a constant
    opt.zero_grad(); (-elbo).backward(); opt.step()
print("posterior mean  exact", m_exact.numpy().round(2), "  VI", m.detach().numpy().round(2))
print("posterior sd    exact", sd_exact.numpy().round(2), "  VI", s.exp().detach().numpy().round(2))
corr = (p * (W[:, 0] - m_exact[0]) * (W[:, 1] - m_exact[1])).sum() / sd_exact.prod()
print("exact posterior correlation", round(corr.item(), 2))

The means agree: (0.30, 0.01) exactly, (0.31, 0.05) by VI. The exact posterior correlation is −0.89, and the VI standard deviations (0.37, 0.33) are half the exact ones, (0.74, 0.68). For a Gaussian with this correlation, Proposition 6.4 predicts the ratio \sqrt{1 - 0.89^2} = 0.46; the posterior is not exactly Gaussian, and the observed ratios are 0.50 and 0.49.

Amortized inference learns a map from each data point to its variational parameters, as in variational auto-encoders, treated in the Probabilistic Graphical Models course. Whatever the algorithm, check a variational fit before trusting it. The Pareto diagnostic of Yao, Vehtari, Simpson and Gelman (2018) fits a generalized Pareto tail to the importance weights p(x, \theta)/q(\theta), \theta \sim q; a shape estimate above 0.7 flags an unreliable approximation. Ch. 9 discusses expectation propagation and the theory of variational Bayes.

7. Bayesian nonparametrics

Parametric models fix their complexity before seeing the data. Bayesian nonparametric models put priors on infinite-dimensional objects, such as distributions or functions, and let the posterior decide how much complexity the data support. The Gaussian process of ch. 3 is one cornerstone; the Dirichlet process is the other.

The Dirichlet process

Proposition 7.1 makes these finite-dimensional laws consistent with each other under merging of cells; a Kolmogorov-type argument, or the explicit construction of Theorem 7.7, gives existence. The base measure G_0 is the mean, and the concentration \alpha is a precision.

New draws repeat earlier values with positive probability, so G is almost surely discrete. Grouping equal values gives the Chinese restaurant process: customer n + 1 joins an existing table c of size n_c with probability n_c / (\alpha + n), and opens a new table with probability \alpha / (\alpha + n).

The mass left after K sticks is \prod_{k \leq K}(1 - V_k), with expectation (\alpha/(1+\alpha))^K by independence. For \alpha = 2, 18 sticks leave an expected mass below 10^{-3}. Truncated stick-breaking is the basis of variational algorithms for DP mixtures (Blei and Jordan, 2006).

Dirichlet process mixtures

The code runs this sampler on 150 points from three unit-variance Gaussians centred at -4, 0, 4, with \alpha = 1 and the base measure \mathcal{N}(0, 25) for the means:

import numpy as np

rng = np.random.default_rng(0)
x = np.concatenate([rng.normal(-4, 1, 50), rng.normal(0, 1, 70), rng.normal(4, 1, 30)])
n, alpha, tau2 = x.size, 1.0, 25.0          # components N(mu, 1), base measure mu ~ N(0, tau2)

def log_pred(xi, s, m):
    """log density of xi given m points with sum s already in the cluster (mu integrated out)."""
    prec = 1 / tau2 + m
    mean, var = s / prec, 1 + 1 / prec
    return -0.5 * (np.log(2 * np.pi * var) + (xi - mean) ** 2 / var)

z = np.zeros(n, dtype=int)                    # start with every point in one cluster
counts, big = {}, {}
for it in range(200):
    for i in range(n):
        z[i] = -1                             # remove point i from its cluster
        labels, sizes = np.unique(z[z >= 0], return_counts=True)
        sums = np.array([x[z == k].sum() for k in labels])
        logw = np.log(np.append(sizes, alpha))                    # CRP prior weights
        logw += np.append(log_pred(x[i], sums, sizes), log_pred(x[i], 0.0, 0))
        w = np.exp(logw - logw.max())
        j = rng.choice(w.size, p=w / w.sum())
        z[i] = labels[j] if j < labels.size else (labels.max() + 1 if labels.size else 0)
    if it >= 50:
        sizes = np.unique(z, return_counts=True)[1]
        counts[sizes.size] = counts.get(sizes.size, 0) + 1
        b = int((sizes >= 5).sum())                    # clusters with at least 5 points
        big[b] = big.get(b, 0) + 1
total = sum(counts.values())
print("posterior of the number of clusters:", {k: round(v / total, 2) for k, v in sorted(counts.items())})
print("clusters with at least 5 points:    ", {k: round(v / total, 2) for k, v in sorted(big.items())})
print("prior mean number of clusters, sum alpha/(alpha+i-1):", round(sum(alpha / (alpha + i) for i in range(n)), 2))

The posterior of the number of clusters is spread from 3 to 11, with mode 6 (probability 0.26) and only 0.08 on the true value 3. Counting only clusters of at least 5 points, the mode moves to 3 (probability 0.43). The extra clusters are small: the DP mixture fits the density well, but its number of clusters is not a consistent estimator of a number of components, as Miller and Harrison (2013) proved. Ch. 9 returns to this point.

Two extensions are used widely. The Pitman–Yor process (Pitman and Yor, 1997) adds a discount parameter and makes the number of clusters grow like a power of n, as word frequencies do. The hierarchical Dirichlet process (Teh, Jordan, Beal and Blei, 2006) shares atoms across groups of data, which gives topic models with an unbounded number of topics. Posterior consistency and contraction rates, the frequentist guarantees of these priors, are the subject of ch. 9.

8. Bayesian methods for deep learning

Deep networks are trained to a single parameter vector and are often overconfident away from their training data. Bayesian deep learning puts a prior on the weights and approximates the posterior predictive. Exact inference is out of reach, so this chapter is about approximations, and about how to tell whether they work.

Training with weight decay is MAP estimation under a Gaussian prior, as ridge regression is in Proposition 3.1. The Bayesian question is what the average adds to the mode.

What a prior on weights means

For the ReLU with u \sim \mathcal{N}(0, I) and no bias, the expectation is half the arc-cosine kernel of degree 1 of Cho and Saul (2009), where \vartheta is the angle between x and x': \mathbb{E}\big[\varphi(u^\top x)\, \varphi(u^\top x')\big] = \frac{\|x\|\, \|x'\|}{2\pi} \big(\sin\vartheta + (\pi - \vartheta)\cos\vartheta\big).

Deep networks have the same limit, layer by layer (Lee et al., 2018; Matthews et al., 2018). In this limit the hidden units do not learn features. At finite width, Vladimirova, Verbeek, Mesejo and Arbel (2019) show that Gaussian weight priors with ReLU-like activations give units of layer \ell a sub-Weibull prior with tail parameter \ell/2: tails get heavier with depth, far from Gaussian.

Approximate posteriors for networks

• Laplace (MacKay, 1992; Daxberger et al., 2021): a Gaussian at the MAP, with the Hessian replaced by a generalized Gauss–Newton or Kronecker-factored approximation, often for the last layer only.
• Variational: a mean-field Gaussian over the weights trained by reparameterization, “Bayes by Backprop” (Blundell et al., 2015). It inherits the small variances of Proposition 6.4.
• Monte Carlo dropout (Gal and Ghahramani, 2016): keep dropout at test time and average stochastic forward passes, which the authors interpret as variational inference with a Bernoulli family.
• Deep ensembles (Lakshminarayanan, Pritzel and Blundell, 2017): train M networks from different initializations and average them. SWAG (Maddox et al., 2019) fits a Gaussian to the iterates of stochastic gradient descent around one solution, and MultiSWAG (Wilson and Izmailov, 2020) averages several such Gaussians.

For classification, combine it with the moderation of Proposition 3.6. ReLU networks trained to the MAP become arbitrarily confident far from the data (Hein, Andriushchenko and Bitterwolf, 2019). Kristiadi, Hein and Hennig (2020) prove that even a Laplace approximation on the last layer keeps the confidence bounded away from 1 there (Exercise 15).

SGLD is the unadjusted Langevin algorithm of ch. 5 with a minibatch gradient. With steps satisfying \sum_t \varepsilon_t = \infty and \sum_t \varepsilon_t^2 < \infty, step-weighted averages are consistent (Teh, Thiery and Vollmer, 2016); with a constant step the chain is biased.

The code trains a five-member deep ensemble on 60 points of \sin x with a gap in (-1, 1), and evaluates it in the data, in the gap and far away:

import numpy as np
import torch

torch.set_num_threads(1)
torch.manual_seed(0)
rng = np.random.default_rng(0)
x = np.concatenate([rng.uniform(-4, -1, 30), rng.uniform(1, 4, 30)])
y = np.sin(x) + 0.1 * rng.normal(size=x.size)
X, Y = torch.tensor(x[:, None], dtype=torch.float32), torch.tensor(y[:, None], dtype=torch.float32)

def train(seed):
    torch.manual_seed(seed)
    net = torch.nn.Sequential(torch.nn.Linear(1, 50), torch.nn.Tanh(),
                              torch.nn.Linear(50, 50), torch.nn.Tanh(), torch.nn.Linear(50, 1))
    opt = torch.optim.Adam(net.parameters(), lr=0.01, weight_decay=1e-4)   # MAP with a Gaussian prior
    for _ in range(1000):
        opt.zero_grad()
        loss = ((net(X) - Y) ** 2).mean()
        loss.backward()
        opt.step()
    return net

nets = [train(s) for s in range(5)]
xs = torch.tensor([[-2.5], [0.0], [8.0]])
with torch.no_grad():
    preds = torch.stack([net(xs) for net in nets]).squeeze(-1)    # (members, points)
print("member predictions at x = -2.5, 0, 8:\n", preds.numpy().round(2))
print("ensemble mean:", preds.mean(0).numpy().round(2), "  true:", np.sin([-2.5, 0, 8]).round(2))
print("ensemble sd:  ", preds.std(0).numpy().round(2))

At x = -2.5 the members agree, with spread 0.02 and an error of 0.04. In the gap, at x = 0, the spread grows to 0.05 while the error is 0.20. Far away, at x = 8, the spread is 0.24 but every member predicts between -2.54 and -1.93, while \sin 8 = 0.99. Diversity among a few MAP fits signals extrapolation, but its size is not calibrated; a Gaussian process with a stationary kernel would return to its prior there (ch. 3).

Evaluating predictive uncertainty

This closes the loop with ch. 1. Forecasting a distribution is a decision problem whose actions are distributions; proper scoring rules are exactly the losses under which the Bayes act is to report one’s actual belief (Gneiting and Raftery, 2007).

Guo, Pleiss, Sun and Weinberger (2017) found modern networks overconfident, and corrected them with a single temperature fitted on validation data. Under dataset shift, Ovadia et al. (2019) found deep ensembles the most robust of the methods above.

The cold posterior puzzle

Wenzel et al. (2020) report that SG-MCMC networks for image classification predict better with T well below 1 than with the Bayes posterior. Three explanations have been put forward: priors on weights are misspecified (Fortuin et al., 2022); data augmentation changes the likelihood, and with full-batch HMC Izmailov et al. (2021) see the effect mostly disappear without augmentation; and curated benchmark labels make the effective likelihood sharper (Aitchison, 2021). The puzzle is a reminder that Bayesian deep learning is approximate Bayes with an imperfect model.

9. Complements for self-study

This chapter collects results that the sessions only mention: the frequentist guarantees of nonparametric posteriors, the number of clusters of DP mixtures, expectation propagation, and what is known about the accuracy of variational Bayes. It assumes ch. 2, 6 and 7.

Posterior consistency

Doob’s theorem is reassuring and weak at the same time. The exceptional set has prior probability zero, but in infinite dimension it can contain natural truths: Diaconis and Freedman (1986) build a reasonable prior for a location problem whose posterior is inconsistent at a simple true distribution. One needs conditions that hold at a given P_0.

Dirichlet process location mixtures of Gaussians satisfy the KL condition at many smooth densities, which is why they are consistent density estimators even when their number of clusters is not consistent.

The rate solves “prior mass near the truth ≈ size of the model”, the same balance as for minimax rates. Gaussian process priors reach the minimax rate when their smoothness matches that of the truth, and lose rate otherwise (van der Vaart and van Zanten, 2008). Dirichlet mixtures of normals reach the parametric rate up to a logarithmic factor for the supersmooth densities they are built for (Ghosal and van der Vaart, 2001). The monograph of Ghosal and van der Vaart (2017) develops all of this.

The number of clusters, revisited

The code of ch. 7 put little posterior mass on the true number of components. This is not a sampling artefact. Miller and Harrison (2013, 2014) prove that, for data from a finite mixture, the posterior probability that a DP or Pitman–Yor mixture has the true number of clusters does not tend to 1. The prior keeps creating small clusters at rate \alpha/(\alpha + n) per point.

Two remedies exist. Count only clusters above a size threshold, as the code did. Or use a mixture of finite mixtures: a prior on the number of components K, then a finite mixture given K. Miller and Harrison (2018) show that this model keeps most of the computational conveniences of the DP, including a restaurant-type sampler, while its posterior on K is consistent under the model.

Heuristically, the expected increment of K_n is about \sigma\, \mathbb{E}[K_n]/n for large n, and the solution of this recursion grows like n^\sigma. Power laws of this kind are observed for word types in text, which made the Pitman–Yor process a natural prior for language models (Teh, 2006).

Expectation propagation

Write the posterior as p(\theta \mid x) \propto p(\theta) \prod_i t_i(\theta) with likelihood factors t_i. Expectation propagation (EP; Minka, 2001) replaces each factor by a “site” \tilde t_i in an exponential family, for instance an unnormalized Gaussian, and sets q \propto p \prod_i \tilde t_i.

Step 3 is therefore a local forward-KL projection, the opposite direction from VI. EP has no convergence guarantee in general, but it is often more accurate than Laplace or mean-field VI for Gaussian process classification (Kuss and Rasmussen, 2005).

How accurate is variational Bayes?

Combine Proposition 6.4 with the Bernstein–von Mises theorem 2.11. In a regular parametric model, the posterior is close to \mathcal{N}(\hat\theta, (nI)^{-1}) with I = I(\theta_0). The mean-field approximation of that Gaussian has the right centre and the variances 1/(n I_{jj}), smaller than the exact marginal variances (I^{-1})_{jj}/n unless I is diagonal.

Wang and Blei (2019) make this rigorous: under regularity conditions the variational posterior concentrates at \theta_0 and its mean is asymptotically normal, but its covariance is too small. Credible sets from mean-field VI are thus too narrow even with infinite data. In nonparametric problems, Alquier and Ridgway (2020) show that variational approximations of tempered posteriors contract at the same rate as the posterior, under a prior-mass condition of the type of Theorem 9.4.

10. Research directions and open problems

This self-study chapter surveys active questions and gives a method for reading papers, the skill the student seminar evaluates. It cites only published work; the selection is partial.

Reading a Bayesian machine learning paper

Most abstracts in the field answer four questions, and locating the answers places the paper on the map of this course.

• Model: what is uncertain, and what are the prior and the likelihood (ch. 2, 3, 4, 7, 8)? Is the prior chosen for its meaning or for convenience?
• Inference: is the posterior computed exactly, sampled (ch. 4, 5), or approximated (ch. 3, 6, 8)? What is the cost per iteration, and what bias remains?
• Decision: which loss or scoring rule is the output judged by (ch. 1, 8)? Is the uncertainty used for a decision, or only reported?
• Guarantee: is there a theorem (exactness, consistency, a rate, a bound) or only experiments? Under which assumptions, and are they checked?

A weak point worth a seminar usually hides in the gap between two answers: an exact sampler evaluated with a misspecified model, or an approximation judged on accuracy rather than calibration.

Monte Carlo at scale

Exact MCMC costs a pass over the data per iteration. Bardenet, Doucet and Holmes (2017) analyse subsampling approaches and show that exact methods need control variates whose quality governs the speed-up. A different line replaces independent draws by repulsive ones: Bardenet and Hardy (2020) use determinantal point processes to build Monte Carlo estimators whose error decreases faster than the usual rate in N^{-1/2}. Non-asymptotic bounds for Langevin algorithms in high dimension (Dalalyan, 2017; Durmus and Moulines, 2017) have made sampling a subject of complexity theory, in parallel with optimization.

Priors for deep networks

Weight-space priors are chosen for convenience, and the cold posterior effect suggests they are wrong (Wenzel et al., 2020). Work by Vladimirova and coauthors (2019) characterizes the induced priors on hidden units, and Fortuin and coauthors (2022) compare heavier-tailed and correlated alternatives. Fortuin (2022) reviews the field. The open question is how to specify priors in function space, where knowledge actually lives, and to transfer them to weights.

Misspecification and generalized posteriors

Bayesian guarantees assume a true model. Grünwald and van Ommen (2017) exhibit inconsistency of Bayesian inference for misspecified linear models and repair it by tempering the likelihood with a learning rate. Bissiri, Holmes and Walker (2016) justify posteriors built from a loss instead of a likelihood, and PAC-Bayes theory bounds the risk of such posteriors without any model assumption (Alquier, 2024). Choosing the temperature in a principled way remains open.

Diffusion models as priors

Score-based generative models give powerful priors for images and signals. Sampling the posterior of an inverse problem with such a prior is hard, because the likelihood enters every step of the denoising process. Chung and coauthors (2023) propose an approximate guidance of the reverse diffusion. Cardoso, Janati El Idrissi, Le Corff and Moulines (2024) instead use sequential Monte Carlo to target the exact posterior of linear inverse problems as the number of particles grows. This connects ch. 5 with Inverse Problems and Imaging and Diffusion and Multiscale Models in AI.

Amortized and simulation-based inference

When the likelihood can be simulated but not evaluated, neural networks trained on simulated pairs (\theta, x) learn the posterior map directly (Cranmer, Brehmer and Louppe, 2020). The same idea, with priors over whole data sets, gives prior-data fitted networks: transformers trained to output posterior predictive distributions (Müller et al., 2022), such as TabPFN for small tabular problems (Hollmann et al., 2023). Xie, Raghunathan, Liang and Ma (2022) read in-context learning of language models as implicit Bayesian inference. Whether such learned inference is calibrated outside its training distribution is open.

Uncertainty without a model: conformal prediction

Conformal prediction builds prediction sets with finite-sample coverage guarantees under exchangeability alone (Angelopoulos and Bates, 2023). It complements Bayesian predictive distributions: Bayesian sets are optimal when the model is right, conformal sets are valid when it is wrong. Combining the two, for instance by conformalizing posterior predictive sets, is an active topic.

Open problems

• Diagnosing an approximate posterior in high dimension without a ground-truth sampler.
• Explaining when deep ensembles approximate the Bayesian model average and when they do not (Wilson and Izmailov, 2020).
• Bernstein–von Mises theorems for high-dimensional and nonparametric functionals, and their failure modes.
• Bayesian inference at the scale of large language models, argued for by Papamarkou and coauthors (2024): which parameters to treat as uncertain, and at what cost.
• Inference on the number of components, clusters or topics that is robust to small departures from the model.

Exercises

Levels: ★ apply the course, ★★ combine results or prove a short statement, ★★★ open-ended or research-flavoured. Exercises 7, 9, 14, 16 and 17 are computational.

1. ★ Three Bayes estimators. Two successes in 10 trials, uniform prior. Give the Bayes estimators under the squared, absolute and 0–1 losses, and the posterior expected squared loss of the first.

   Solution: the posterior is Beta(3, 9). Its mean is 3/12 = 0.25, its median 0.236 and its mode 2/10 = 0.2 (Proposition 1.4).

   The posterior expected squared loss of the mean is the posterior variance, ab/((a+b)^2(a+b+1)) = 27/1872 = 0.0144, a standard deviation of 0.12.

2. ★ A cost-sensitive decision. A screening model gives a predictive probability of disease of 0.15. Missing a disease costs 20, a false alarm costs 1. Should the patient be flagged? Where does the decision flip? A review by a specialist, who then settles the case without error, costs 0.5: is it worth it?

   Solution: by Proposition 1.12 the threshold is 1/(1 + 20) = 0.048, and 0.15 > 0.048: flag.

   The expected losses are 1 \times 0.85 = 0.85 for flagging and 20 \times 0.15 = 3 for not flagging. Since the review settles the case with no further loss, its expected loss is 0.5, below both, so referring the case is the Bayes action.

3. ★ Gamma–Poisson update. Daily counts (2, 0, 3, 1, 4) follow a Poisson law with rate \lambda, and the prior is Gamma(2, 1) (shape, rate). Give the posterior, its mean and standard deviation, and the predictive probability of a day with no event.

   Solution: the posterior is Gamma(2 + 10, 1 + 5), that is Gamma(12, 6), with mean 2 and standard deviation \sqrt{12}/6 = 0.577.

   The predictive law is negative binomial, and P(X^\star = 0 \mid x) = \mathbb{E}[e^{-\lambda} \mid x] = (6/7)^{12} = 0.157 .

4. ★ Jeffreys priors and proper posteriors. Compute the Jeffreys prior of the Poisson rate, of the exponential rate and of a Gaussian scale \sigma with known mean 0. After a single observation, when is the posterior proper?

   Solution: Poisson: I(\lambda) = 1/\lambda, prior \lambda^{-1/2}, posterior Gamma(x + 1/2, 1), always proper.

   Exponential: \log p = \log\lambda - \lambda x, so I = 1/\lambda^2, prior 1/\lambda, posterior proportional to e^{-\lambda x}, proper when x > 0.

   Gaussian scale: I = 2/\sigma^2, prior 1/\sigma, posterior proportional to \sigma^{-2} e^{-x^2/2\sigma^2}. The substitution u = 1/\sigma shows it is proper when x \neq 0.

5. ★ Bayesian regression by hand. Redo Example 3.2 with a third point (3, 3.2). Give the posterior and the predictive distribution at x^\star = 4.

   Solution: \sum x_i^2 = 14 and \sum x_i y_i = 14.6. The posterior precision is 1 + 4 \times 14 = 57, the posterior mean 4 \times 14.6/57 = 1.025, the standard deviation 0.132.

   At x^\star = 4: predictive mean 4.098, variance 0.25 + 16/57 = 0.531, standard deviation 0.728. The parameter part of the predictive variance is now x^{\star 2}/57 instead of x^{\star 2}/21: at x^\star = 3 it falls from 0.43 to 0.16.

6. ★ Chinese restaurant arithmetic. For n = 10 customers, compute the expected number of tables for \alpha \in \{0.1, 1, 10\}, and the probability that all customers sit at one table when \alpha = 1.

   Solution: by Example 7.6, the expectations are 1.27, 2.93 and 7.19.

   The probability of a single table is the product over customers 2 to 10 of (i - 1)/(\alpha + i - 1), which telescopes to 1/10 for \alpha = 1.

7. ★ Coverage at n = 100 (Python). Rerun the coverage code of ch. 2 with n = 100. What changes for each interval?

   Hints and expected result: mean coverages become 0.931 (Wald), 0.950 (Jeffreys) and 0.950 (uniform); minimum coverages 0.633, 0.881 and 0.921.

   The Wald interval still fails near the boundaries, where the normal approximation of the binomial is poor. The two credible intervals are close to nominal on average, as Theorem 2.11 predicts, and oscillate because the binomial is discrete.

8. ★★ The Bayes risk of the posterior mean. Show that under the squared loss the Bayes risk is \mathbb{E}[\mathrm{Var}(\theta \mid X)]. For the binomial model with a uniform prior, show that it equals 1/(6(n+2)), and compare with the average risk 1/(6n) of K/n.

   Solution: by Theorem 1.3 and Proposition 1.4, the Bayes risk is the average over X of the posterior expected loss of the posterior mean, which is the posterior variance.

   With a uniform prior, K is uniform on \{0, \dots, n\}, and the posterior Beta(K + 1, n - K + 1) has variance \frac{(K+1)(n-K+1)}{(n+2)^2(n+3)} . The sum of j(n + 2 - j) over j = 1, \dots, n+1 is (n+1)(n+2)(n+3)/6, which gives 1/(6(n+2)).

   For n = 16 this is 0.00926, against 1/96 = 0.0104 for K/n: the values printed by the code of ch. 1.

9. ★★ James–Stein risk (Python). For d = 10, estimate by simulation the risk of \delta^{\mathrm{JS}} and of its positive part, which replaces the factor by its positive part, at \|\theta\| \in \{0, 2, 5, 10\}. Check the formula of Theorem 1.11.

   Hints and expected result: with 2 \times 10^5 replications, the risk of \delta^{\mathrm{JS}} is about 2.00, 4.38, 7.97 and 9.38, always below 10, and matches d - (d-2)^2\, \mathbb{E}\|X\|^{-2} to within simulation error.

   The positive-part version gives 1.26 and 3.98 at the first two points and is practically identical further away, where the factor is almost never negative (probability 4 \times 10^{-4} at \|\theta\| = 5). It dominates \delta^{\mathrm{JS}}, a classical result (Lehmann and Casella, 1998).

10. ★★ Mean field in both directions. Let p = \mathcal{N}(0, \Sigma) in two dimensions with unit variances and correlation \rho. Find the product of Gaussians minimizing \mathrm{KL}(q \,\|\, p), then the product of densities minimizing \mathrm{KL}(p \,\|\, q). Evaluate for \rho = 0.9.

    Solution: the first is given by Proposition 6.4: centred, variances 1 - \rho^2.

    For the second, with p_1, p_2 the marginals of p and H the entropy, write \mathrm{KL}(p \,\|\, q_1 q_2) = -H(p) + H(p_1) + H(p_2) + \mathrm{KL}(p_1 \,\|\, q_1) + \mathrm{KL}(p_2 \,\|\, q_2) . It is minimized at q_j = p_j, with unit variances; the minimum is the mutual information -\frac12 \log(1 - \rho^2).

    For \rho = 0.9: standard deviations 0.436 against 1, and a minimum of 0.830 nats.

11. ★★ Gibbs is Metropolis–Hastings with acceptance 1. Show that the Gibbs update of block j is a Metropolis–Hastings step whose proposal is the full conditional, and that its acceptance probability is always 1.

    Solution: the proposal changes only \theta_j, with density \pi(\theta_j' \mid \theta_{-j}). The Metropolis–Hastings ratio is \frac{\pi(\theta_j', \theta_{-j})\, \pi(\theta_j \mid \theta_{-j})}{\pi(\theta_j, \theta_{-j})\, \pi(\theta_j' \mid \theta_{-j})} = \frac{\pi(\theta_{-j})\, \pi(\theta_j' \mid \theta_{-j})\, \pi(\theta_j \mid \theta_{-j})}{\pi(\theta_{-j})\, \pi(\theta_j \mid \theta_{-j})\, \pi(\theta_j' \mid \theta_{-j})} = 1 . Theorem 5.4 then gives a second proof of Proposition 4.3.

12. ★★ Moments of the Dirichlet process. Let G \sim \mathrm{DP}(\alpha, G_0) with G_0 without atoms.

    (a) Show that two draws from G coincide with probability 1/(1+\alpha).

    (b) For disjoint A and B, compute \mathrm{Cov}(G(A), G(B)).

    (c) How many sticks leave an expected mass below 10^{-2} when \alpha = 5?

    Solution: (a) by the Pólya urn, the second draw equals the first with probability 1/(\alpha + 1).

    (b) With C the complement of A \cup B, the triple (G(A), G(B), G(C)) is Dirichlet, and \mathbb{E}[W_1 W_2] = a_1 a_2 / (a_0(a_0 + 1)) with a_0 = \alpha. Hence \mathrm{Cov}\big(G(A), G(B)\big) = -\frac{G_0(A)\, G_0(B)}{1 + \alpha} .

    (c) (5/6)^K \leq 10^{-2} requires K \geq 25.3, so 26 sticks.

13. ★★ Laplace approximations of a Beta posterior. For the posterior Beta(8, 4) of Example 1.5, compute the Laplace approximation in \theta and in the logit \varphi = \log(\theta/(1-\theta)), and compare the resulting 95 % intervals with the exact one.

    Solution: in \theta, the log density 7\log\theta + 3\log(1-\theta) has its mode at 0.7 and curvature 7/0.49 + 3/0.09 = 47.6, so the variance is 0.021 and the interval (0.416, 0.984).

    In \varphi, the Jacobian turns the density into \theta^8(1-\theta)^4, with mode \theta = 2/3, \varphi = \log 2, and curvature 12 \times \frac23 \times \frac13 = 2.67.

    The interval \log 2 \pm 1.96/\sqrt{2.67} maps to (0.376, 0.869) on the scale of \theta.

    The exact interval is (0.390, 0.891): the logit approximation is much closer, because the posterior is more Gaussian on an unbounded scale. Laplace depends on the parameterization.

14. ★★ A Gibbs sampler for regression with unknown noise (Python). In the evidence code of ch. 3, take degree 3, keep \alpha = 1, and put a Gamma(1, 0.01) prior (shape, rate) on the noise precision \beta. Derive the two full conditionals, implement the sampler, and report the posterior of the noise standard deviation.

    Hints and expected result: w \mid \beta, y is given by Proposition 3.1, and \beta \mid w, y \sim Gamma(1 + n/2,\; 0.01 + \|y - \Phi w\|^2/2).

    With 5000 iterations, the noise standard deviation has posterior mean about 0.32 and 90 % interval about (0.23, 0.44): larger than the 0.2 used to simulate. A cubic cannot represent \sin(\pi x) exactly (training RMSE 0.26 in ch. 3), and the prior shrinks the large coefficients. The posterior means of the coefficients are about (-0.07, 2.24, 0.04, -2.19).

15. ★★★ Overconfidence far from the data. For a binary classifier whose logit is w^\top h(x), with h the last hidden layer of a ReLU network, show that along almost every ray t \mapsto t x the MAP confidence tends to 1. Then show that a Gaussian posterior \mathcal{N}(m, S) on w, with S positive definite, combined with Proposition 3.6, keeps the confidence below 1.

    Hints and expected result: h is piecewise affine with finitely many regions, so for large t one has h(tx) = t\, u + c on a single region. The logit is t\, w^\top u + w^\top c, which tends to \pm\infty unless w^\top u = 0 (Hein, Andriushchenko and Bitterwolf, 2019).

    Under the Gaussian posterior, \mu_a \approx t\, m^\top u and \sigma_a^2 \approx t^2\, u^\top S u, so the moderated probability tends to \sigma\Big(\frac{m^\top u}{\sqrt{\pi\, u^\top S u / 8}}\Big) < 1 . This is the mechanism of Kristiadi, Hein and Hennig (2020).

16. ★★★ Cold posteriors on a toy problem (Python, open-ended). On the two-moons data of scikit-learn (make_moons, 200 training points, noise 0.2), sample a two-layer network with SGLD at temperatures T \in \{1, 0.3, 0.1\}. Compare test accuracy and log score. Then add “augmented” jittered copies of each training point, counted as independent observations, and compare again.

    Hints and expected result: tempering multiplies the gradient by 1/T and keeps the noise. On this small, well-specified problem, expect T = 1 to be competitive in log score.

    With augmentation counted as extra data, the likelihood is over-counted and the ranking of temperatures moves. Relate what you see to the explanations of Wenzel et al. (2020) and Izmailov et al. (2021), and report several seeds with error bars.

17. ★★★ Is the number of clusters consistent? (Python, open-ended). Run the DP mixture code of ch. 7 on samples of size 150, 600 and 1200 from the same three-component mixture. Track the posterior probability of exactly 3 clusters, and of 3 clusters holding at least 5 % of the points. Then implement the restaurant process of a mixture of finite mixtures (Miller and Harrison, 2018) and compare.

    Hints and expected result: the probability of exactly 3 clusters does not increase with n, because small clusters keep appearing, while the count of large clusters concentrates at 3.

    The mixture of finite mixtures replaces the weights n_c and \alpha by n_c + \gamma and a ratio of coefficients that depends on the prior on K.

    Its posterior on K concentrates at 3 as n grows. Profile the code before running large n.

Projects and assessment

The official assessment, per master-mva.com, is a group study of a research paper chosen from a list. Each group writes a small report and gives a short oral presentation in front of the class, in a student seminar during the last lecture. The description strongly encourages a dash of creativity: find a weak point, shortcoming or limitation of the paper and push in that direction, by extending a proof, implementing another feature or running different experiments. The Université Paris-Saclay course detail classifies the assessment as “Evaluation Continue non Intégrale”.

All the subjects below are paper studies in that spirit, with a weak point to push. The solo subjects are sized for one student; the group subjects match the official format. The written exam is not part of the official assessment: it is offered as a self-test, and none of these proposals is official material.

Solo projects

1. [Solo · ≈ 20 h] Stochastic-gradient Langevin dynamics: a paper study. Study Welling and Teh (2011) and reimplement SGLD (Algorithm 8.4).
   • Experiments: the two-parameter Gaussian mixture example of the paper (synthetic data), then Bayesian logistic regression on the breast-cancer data set bundled with scikit-learn.
   • Weak point to push: the bias of constant step sizes. Measure the error on posterior means and variances against a long MALA chain (ch. 5), as a function of step size and batch size.
   • Success criterion: a plot of the error against the step size showing the bias, and a comparison of effective sample size per second with MALA.
   • Deliverables: code, a 4–6 page report, a 10-minute talk. Chapters 5, 8.
2. [Solo · ≈ 20 h] Dropout, ensembles and Laplace: a paper study. Compare Monte Carlo dropout (Gal and Ghahramani, 2016) with deep ensembles (Lakshminarayanan, Pritzel and Blundell, 2017) and a last-layer Laplace approximation (Proposition 8.3).
   • Data: the one-dimensional problem with a gap of ch. 8, and the diabetes regression data set bundled with scikit-learn.
   • Weak point to push: a posterior should contract as data accumulate. Measure how each method’s predictive variance changes when the training set grows from 50 to 400 points.
   • Success criterion: a table of test log score, RMSE and 90 % interval coverage, and a plot of predictive variance against n for each method.
   • Deliverables: code, report, talk. Chapters 3, 6, 8.
3. [Solo · ≈ 15 h] Objective priors at work: interval estimation for a proportion. Reproduce the coverage figures of Brown, Cai and DasGupta (2001) for the Wald, Wilson, Clopper–Pearson and Jeffreys intervals.
   • Weak point to push: extend the study to a Poisson rate with the Jeffreys prior \lambda^{-1/2}, or to the difference of two proportions, where no interval is exact.
   • Success criterion: coverage curves over a fine grid of parameters, with an explanation of the oscillations and a recommendation backed by numbers.
   • Deliverables: code, report, talk. Chapters 1, 2.

Group projects

1. [Group of 2–3 · ≈ 30 h each] Dirichlet process mixtures: sampling against variational inference. Study Blei and Jordan (2006) together with Neal (2000).
   • Work split: one member implements the collapsed Gibbs sampler (Algorithm 7.9); one the truncated stick-breaking variational algorithm; one the mixture of finite mixtures of Miller and Harrison (2018).
   • Data: synthetic Gaussian mixtures, and the digits data set of scikit-learn projected on its first principal components.
   • Weak point to push: the number of clusters (ch. 9). Show how each method’s posterior on the number of clusters behaves as n grows.
   • Success criterion: held-out log predictive density and running time for each method, and a figure reproducing the inconsistency proved by Miller and Harrison (2013).
   • Deliverables: code, a 6-page report, a 15-minute seminar talk. Chapters 4, 6, 7, 9.
2. [Group of 2–3 · ≈ 30 h each] What do Bayesian neural network posteriors look like? A scaled-down replication. Study Izmailov, Vikram, Hoffman and Wilson (2021) and Wenzel et al. (2020).
   • Setting: a small multilayer perceptron on the scikit-learn digits data set, small enough for full-batch HMC to serve as reference.
   • Work split: one member runs the HMC reference (ch. 5); one SGLD at several temperatures (Algorithm 8.4, Definition 8.8); one mean-field VI and deep ensembles, and the shared evaluation.
   • Success criterion: a table of accuracy, log score and expected calibration error, the agreement of each method’s predictive distribution with HMC, and a cold-posterior curve with and without data augmentation.
   • Deliverables: code, report, seminar talk. Chapters 5, 6, 8.
3. [Group of 2–3 · ≈ 30 h each] Markov chain Monte Carlo for tall data. Study Bardenet, Doucet and Holmes (2017).
   • Setting: Bayesian logistic regression on 10^5 synthetic points.
   • Work split: one member implements full Metropolis–Hastings and naive subsampling; one the confidence sampler with control variates of the paper; one SGLD and the evaluation.
   • Weak point to push: the speed-up when the posterior is far from Gaussian, for instance with heavy-tailed covariates or a misspecified model.
   • Success criterion: the number of likelihood evaluations per iteration, and the error on posterior means and variances against full Metropolis–Hastings.
   • Deliverables: code, report, seminar talk. Chapters 5, 8.

Written exam (proposal)

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

1. Course question (4 points). Define the frequentist risk, the integrated risk and a Bayes rule. Prove that an action minimizing the posterior expected loss for each observation gives a Bayes rule. State two conditions under which a Bayes rule is admissible.

   Expected: Definitions 1.1–1.2; Theorem 1.3 with the Fubini exchange and the pointwise inequality; Proposition 1.7, uniqueness or continuity of risks with a full-support prior.

2. Computation (6 points). Daily counts (3, 1, 4, 0, 2) are Poisson with rate \lambda, and the prior is Gamma(1, 1).

   (a) Give the posterior. (b) Give the Bayes estimators under the squared loss and under the loss (\lambda - a)^2/\lambda. (c) Compute the predictive probability of a day with no event. (d) Redo (a) with the Jeffreys prior.

   Expected: (a) Gamma(11, 6). (b) The posterior mean 11/6 = 1.83; for the weighted loss, minimizing the posterior expectation \mathbb{E}[\lambda] - 2a + a^2\, \mathbb{E}[1/\lambda] gives a = \frac{1}{\mathbb{E}[1/\lambda \mid x]} = \frac{11 - 1}{6} = 1.67 . (c) (6/7)^{11} = 0.184. (d) The prior \lambda^{-1/2} gives Gamma(10.5, 5), with mean 2.1.

3. Short proof (5 points). Prove the ELBO identity and the coordinate-ascent update of mean-field variational inference. Deduce that for a Gaussian target the mean-field variances are the inverse diagonal entries of the precision matrix, and that they underestimate the marginal variances.

   Expected: Proposition 6.1, Theorem 6.2 with the KL to the tilted factor, Example 6.3, and Proposition 6.4 through the Schur complement.

4. Problem (5 points). Let G \sim \mathrm{DP}(\alpha, G_0) with G_0 without atoms.

   (a) Compute the mean and variance of G(A). (b) Derive the predictive rule of X_{n+1} given X_{1:n} from the posterior. (c) Show that the expected number of distinct values among n draws grows like \alpha\log n. (d) Explain in two sentences what changes with a Pitman–Yor process.

   Expected: (a) Proposition 7.3 through the Beta marginal. (b) Theorem 7.4 and Proposition 7.5. (c) The sum of \alpha/(\alpha + i - 1), compared with an integral. (d) A discount \sigma makes new clusters more likely as K grows, so the number of clusters grows like n^\sigma (Proposition 9.5).

Grading rubric (proposal)

References

Official references (master-mva.com)

• G. Parmigiani and L. Inoue. Decision Theory: Principles and Approaches. Wiley, 2009.
• C. P. Robert. The Bayesian Choice: From Decision-Theoretic Foundations to Computational Implementation, 2nd edition. Springer, 2007.
• K. P. Murphy. Machine Learning: A Probabilistic Perspective. MIT Press, 2012.
• S. Ghosal and A. W. van der Vaart. Fundamentals of Nonparametric Bayesian Inference. Cambridge University Press, 2017.

Textbooks

• J. O. Berger. Statistical Decision Theory and Bayesian Analysis, 2nd edition. Springer, 1985.
• J. M. Bernardo and A. F. M. Smith. Bayesian Theory. Wiley, 1994.
• C. M. Bishop. Pattern Recognition and Machine Learning. Springer, 2006.
• N. Chopin and O. Papaspiliopoulos. An Introduction to Sequential Monte Carlo. Springer, 2020.
• A. Gelman, J. B. Carlin, H. S. Stern, D. B. Dunson, A. Vehtari and D. B. Rubin. Bayesian Data Analysis, 3rd edition. CRC Press, 2013.
• N. L. Hjort, C. Holmes, P. Müller and S. G. Walker (eds.). Bayesian Nonparametrics. Cambridge University Press, 2010.
• E. L. Lehmann and G. Casella. Theory of Point Estimation, 2nd edition. Springer, 1998.
• S. Meyn and R. L. Tweedie. Markov Chains and Stochastic Stability, 2nd edition. Cambridge University Press, 2009.
• K. P. Murphy. Probabilistic Machine Learning: An Introduction. MIT Press, 2022.
• K. P. Murphy. Probabilistic Machine Learning: Advanced Topics. MIT Press, 2023.
• R. M. Neal. Bayesian Learning for Neural Networks. Lecture Notes in Statistics 118, Springer, 1996.
• J. Pitman. Combinatorial Stochastic Processes. Lecture Notes in Mathematics 1875, Springer, 2006.
• C. E. Rasmussen and C. K. I. Williams. Gaussian Processes for Machine Learning. MIT Press, 2006.
• C. P. Robert and G. Casella. Monte Carlo Statistical Methods, 2nd edition. Springer, 2004.
• L. J. Savage. The Foundations of Statistics. Wiley, 1954.
• A. W. van der Vaart. Asymptotic Statistics. Cambridge University Press, 1998.
• J. von Neumann and O. Morgenstern. Theory of Games and Economic Behavior. Princeton University Press, 1944.
• A. Wald. Statistical Decision Functions. Wiley, 1950.

Decision theory and the interpretations of probability (ch. 1–2)

• J. M. Bernardo. Reference posterior distributions for Bayesian inference. Journal of the Royal Statistical Society B, 41(2):113–147, 1979.
• J. O. Berger, J. M. Bernardo and D. Sun. The formal definition of reference priors. Annals of Statistics, 37(2):905–938, 2009.
• A. Birnbaum. On the foundations of statistical inference. Journal of the American Statistical Association, 57(298):269–306, 1962.
• L. D. Brown, T. T. Cai and A. DasGupta. Interval estimation for a binomial proportion. Statistical Science, 16(2):101–133, 2001.
• B. de Finetti. La prévision : ses lois logiques, ses sources subjectives. Annales de l’Institut Henri Poincaré, 7(1):1–68, 1937.
• B. Efron and C. Morris. Stein’s estimation rule and its competitors — an empirical Bayes approach. Journal of the American Statistical Association, 68(341):117–130, 1973.
• E. Hewitt and L. J. Savage. Symmetric measures on Cartesian products. Transactions of the American Mathematical Society, 80(2):470–501, 1955.
• W. James and C. Stein. Estimation with quadratic loss. Proceedings of the Fourth Berkeley Symposium on Mathematical Statistics and Probability, 1:361–379, 1961.
• H. Jeffreys. An invariant form for the prior probability in estimation problems. Proceedings of the Royal Society A, 186(1007):453–461, 1946.
• B. J. K. Kleijn and A. W. van der Vaart. The Bernstein–von-Mises theorem under misspecification. Electronic Journal of Statistics, 6:354–381, 2012.
• D. V. Lindley and L. D. Phillips. Inference for a Bernoulli process (a Bayesian view). The American Statistician, 30(3):112–119, 1976.

Bayesian supervised and unsupervised learning (ch. 3–4)

• C. M. Bishop. Bayesian PCA. Advances in Neural Information Processing Systems 11, 1999.
• D. M. Blei, A. Y. Ng and M. I. Jordan. Latent Dirichlet allocation. Journal of Machine Learning Research, 3:993–1022, 2003.
• G. E. P. Box. Sampling and Bayes’ inference in scientific modelling and robustness. Journal of the Royal Statistical Society A, 143(4):383–430, 1980.
• G. Celeux, M. Hurn and C. P. Robert. Computational and inferential difficulties with mixture posterior distributions. Journal of the American Statistical Association, 95(451):957–970, 2000.
• A. Gelman, X.-L. Meng and H. Stern. Posterior predictive assessment of model fitness via realized discrepancies. Statistica Sinica, 6(4):733–807, 1996.
• T. L. Griffiths and M. Steyvers. Finding scientific topics. Proceedings of the National Academy of Sciences, 101(suppl. 1):5228–5235, 2004.
• R. E. Kass and A. E. Raftery. Bayes factors. Journal of the American Statistical Association, 90(430):773–795, 1995.
• D. V. Lindley. A statistical paradox. Biometrika, 44(1–2):187–192, 1957.
• D. J. C. MacKay. Bayesian interpolation. Neural Computation, 4(3):415–447, 1992.
• D. J. C. MacKay. The evidence framework applied to classification networks. Neural Computation, 4(5):720–736, 1992.
• J. Rousseau and K. Mengersen. Asymptotic behaviour of the posterior distribution in overfitted mixture models. Journal of the Royal Statistical Society B, 73(5):689–710, 2011.
• G. Schwarz. Estimating the dimension of a model. Annals of Statistics, 6(2):461–464, 1978.
• M. Stephens. Dealing with label switching in mixture models. Journal of the Royal Statistical Society B, 62(4):795–809, 2000.
• L. Tierney and J. B. Kadane. Accurate approximations for posterior moments and marginal densities. Journal of the American Statistical Association, 81(393):82–86, 1986.
• M. E. Tipping and C. M. Bishop. Probabilistic principal component analysis. Journal of the Royal Statistical Society B, 61(3):611–622, 1999.

Monte Carlo methods (ch. 5)

• R. Bardenet, A. Doucet and C. Holmes. On Markov chain Monte Carlo methods for tall data. Journal of Machine Learning Research, 18(47):1–43, 2017.
• A. Beskos, N. Pillai, G. Roberts, J.-M. Sanz-Serna and A. Stuart. Optimal tuning of the hybrid Monte Carlo algorithm. Bernoulli, 19(5A):1501–1534, 2013.
• 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.
• P. Del Moral, A. Doucet and A. Jasra. Sequential Monte Carlo samplers. Journal of the Royal Statistical Society B, 68(3):411–436, 2006.
• S. Duane, A. D. Kennedy, B. J. Pendleton and D. Roweth. Hybrid Monte Carlo. Physics Letters B, 195(2):216–222, 1987.
• A. Durmus and É. Moulines. Nonasymptotic convergence analysis for the unadjusted Langevin algorithm. Annals of Applied Probability, 27(3):1551–1587, 2017.
• A. Gelman and D. B. Rubin. Inference from iterative simulation using multiple sequences. Statistical Science, 7(4):457–472, 1992.
• S. Geman and 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.
• W. K. Hastings. Monte Carlo sampling methods using Markov chains and their applications. Biometrika, 57(1):97–109, 1970.
• M. D. Hoffman and A. Gelman. The No-U-Turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo. Journal of Machine Learning Research, 15:1593–1623, 2014.
• N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller and E. Teller. Equation of state calculations by fast computing machines. Journal of Chemical Physics, 21(6):1087–1092, 1953.
• R. M. Neal. Annealed importance sampling. Statistics and Computing, 11(2):125–139, 2001.
• R. M. Neal. MCMC using Hamiltonian dynamics. In Handbook of Markov Chain Monte Carlo, Chapman & Hall/CRC, 2011.
• G. O. Roberts, A. Gelman and W. R. Gilks. Weak convergence and optimal scaling of random walk Metropolis algorithms. Annals of Applied Probability, 7(1):110–120, 1997.
• G. O. Roberts and J. S. Rosenthal. Optimal scaling of discrete approximations to Langevin diffusions. Journal of the Royal Statistical Society B, 60(1):255–268, 1998.
• A. Vehtari, A. Gelman, D. Simpson, B. Carpenter and 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.

Variational methods (ch. 6 and 9)

• P. Alquier and J. Ridgway. Concentration of tempered posteriors and of their variational approximations. Annals of Statistics, 48(3):1475–1497, 2020.
• D. M. Blei, A. Kucukelbir and J. D. McAuliffe. Variational inference: a review for statisticians. Journal of the American Statistical Association, 112(518):859–877, 2017.
• M. D. Hoffman, D. M. Blei, C. Wang and J. Paisley. Stochastic variational inference. Journal of Machine Learning Research, 14:1303–1347, 2013.
• M. I. Jordan, Z. Ghahramani, T. S. Jaakkola and L. K. Saul. An introduction to variational methods for graphical models. Machine Learning, 37(2):183–233, 1999.
• D. P. Kingma and M. Welling. Auto-encoding variational Bayes. International Conference on Learning Representations, 2014.
• M. Kuss and C. E. Rasmussen. Assessing approximate inference for binary Gaussian process classification. Journal of Machine Learning Research, 6:1679–1704, 2005.
• T. P. Minka. Expectation propagation for approximate Bayesian inference. Proceedings of Uncertainty in Artificial Intelligence, 2001.
• R. Ranganath, S. Gerrish and D. M. Blei. Black box variational inference. Proceedings of AISTATS, 2014.
• D. J. Rezende, S. Mohamed and D. Wierstra. Stochastic backpropagation and approximate inference in deep generative models. Proceedings of ICML, 2014.
• Y. Wang and D. M. Blei. Frequentist consistency of variational Bayes. Journal of the American Statistical Association, 114(527):1147–1161, 2019.
• Y. Yao, A. Vehtari, D. Simpson and A. Gelman. Yes, but did it work?: Evaluating variational inference. Proceedings of ICML, 2018.

Bayesian nonparametrics (ch. 7 and 9)

• C. E. Antoniak. Mixtures of Dirichlet processes with applications to Bayesian nonparametric problems. Annals of Statistics, 2(6):1152–1174, 1974.
• D. Blackwell and J. B. MacQueen. Ferguson distributions via Pólya urn schemes. Annals of Statistics, 1(2):353–355, 1973.
• D. M. Blei and M. I. Jordan. Variational inference for Dirichlet process mixtures. Bayesian Analysis, 1(1):121–143, 2006.
• P. Diaconis and D. Freedman. On the consistency of Bayes estimates. Annals of Statistics, 14(1):1–26, 1986.
• J. L. Doob. Application of the theory of martingales. In Le Calcul des Probabilités et ses Applications, Colloques Internationaux du CNRS 13, 23–27, 1949.
• T. S. Ferguson. A Bayesian analysis of some nonparametric problems. Annals of Statistics, 1(2):209–230, 1973.
• S. Ghosal, J. K. Ghosh and A. W. van der Vaart. Convergence rates of posterior distributions. Annals of Statistics, 28(2):500–531, 2000.
• S. Ghosal and A. W. van der Vaart. Entropies and rates of convergence for maximum likelihood and Bayes estimation for mixtures of normal densities. Annals of Statistics, 29(5):1233–1263, 2001.
• A. Y. Lo. On a class of Bayesian nonparametric estimates: I. Density estimates. Annals of Statistics, 12(1):351–357, 1984.
• J. W. Miller and M. T. Harrison. A simple example of Dirichlet process mixture inconsistency for the number of components. Advances in Neural Information Processing Systems 26, 2013.
• J. W. Miller and M. T. Harrison. Inconsistency of Pitman–Yor process mixtures for the number of components. Journal of Machine Learning Research, 15:3333–3370, 2014.
• J. W. Miller and M. T. Harrison. Mixture models with a prior on the number of components. Journal of the American Statistical Association, 113(521):340–356, 2018.
• R. M. Neal. Markov chain sampling methods for Dirichlet process mixture models. Journal of Computational and Graphical Statistics, 9(2):249–265, 2000.
• J. Pitman and M. Yor. The two-parameter Poisson–Dirichlet distribution derived from a stable subordinator. Annals of Probability, 25(2):855–900, 1997.
• L. Schwartz. On Bayes procedures. Zeitschrift für Wahrscheinlichkeitstheorie und verwandte Gebiete, 4(1):10–26, 1965.
• J. Sethuraman. A constructive definition of Dirichlet priors. Statistica Sinica, 4(2):639–650, 1994.
• Y. W. Teh. A hierarchical Bayesian language model based on Pitman–Yor processes. Proceedings of COLING/ACL, 2006.
• Y. W. Teh, M. I. Jordan, M. J. Beal and D. M. Blei. Hierarchical Dirichlet processes. Journal of the American Statistical Association, 101(476):1566–1581, 2006.
• A. W. van der Vaart and J. H. van Zanten. Rates of contraction of posterior distributions based on Gaussian process priors. Annals of Statistics, 36(3):1435–1463, 2008.

Bayesian deep learning (ch. 8)

• L. Aitchison. A statistical theory of cold posteriors in deep neural networks. International Conference on Learning Representations, 2021.
• C. Blundell, J. Cornebise, K. Kavukcuoglu and D. Wierstra. Weight uncertainty in neural networks. Proceedings of ICML, 2015.
• Y. Cho and L. K. Saul. Kernel methods for deep learning. Advances in Neural Information Processing Systems 22, 2009.
• E. Daxberger, A. Kristiadi, A. Immer, R. Eschenhagen, M. Bauer and P. Hennig. Laplace redux — effortless Bayesian deep learning. Advances in Neural Information Processing Systems 34, 2021.
• V. Fortuin, A. Garriga-Alonso, S. W. Ober, F. Wenzel, G. Rätsch, R. E. Turner, M. van der Wilk and L. Aitchison. Bayesian neural network priors revisited. International Conference on Learning Representations, 2022.
• Y. Gal and Z. Ghahramani. Dropout as a Bayesian approximation: representing model uncertainty in deep learning. Proceedings of ICML, 2016.
• T. Gneiting and A. E. Raftery. Strictly proper scoring rules, prediction, and estimation. Journal of the American Statistical Association, 102(477):359–378, 2007.
• C. Guo, G. Pleiss, Y. Sun and K. Q. Weinberger. On calibration of modern neural networks. Proceedings of ICML, 2017.
• M. Hein, M. Andriushchenko and J. Bitterwolf. Why ReLU networks yield high-confidence predictions far away from the training data and how to mitigate the problem. Proceedings of CVPR, 2019.
• P. Izmailov, S. Vikram, M. D. Hoffman and A. G. Wilson. What are Bayesian neural network posteriors really like? Proceedings of ICML, 2021.
• A. Kristiadi, M. Hein and P. Hennig. Being Bayesian, even just a bit, fixes overconfidence in ReLU networks. Proceedings of ICML, 2020.
• B. Lakshminarayanan, A. Pritzel and C. Blundell. Simple and scalable predictive uncertainty estimation using deep ensembles. Advances in Neural Information Processing Systems 30, 2017.
• J. Lee, Y. Bahri, R. Novak, S. S. Schoenholz, J. Pennington and J. Sohl-Dickstein. Deep neural networks as Gaussian processes. International Conference on Learning Representations, 2018.
• W. J. Maddox, T. Garipov, P. Izmailov, D. Vetrov and A. G. Wilson. A simple baseline for Bayesian uncertainty in deep learning. Advances in Neural Information Processing Systems 32, 2019.
• A. G. de G. Matthews, M. Rowland, J. Hron, R. E. Turner and Z. Ghahramani. Gaussian process behaviour in wide deep neural networks. International Conference on Learning Representations, 2018.
• Y. Ovadia et al. Can you trust your model’s uncertainty? Evaluating predictive uncertainty under dataset shift. Advances in Neural Information Processing Systems 32, 2019.
• Y. W. Teh, A. H. Thiery and S. J. Vollmer. Consistency and fluctuations for stochastic gradient Langevin dynamics. Journal of Machine Learning Research, 17(7):1–33, 2016.
• M. Vladimirova, J. Verbeek, P. Mesejo and J. Arbel. Understanding priors in Bayesian neural networks at the unit level. Proceedings of ICML, 2019.
• M. Welling and Y. W. Teh. Bayesian learning via stochastic gradient Langevin dynamics. Proceedings of ICML, 2011.
• F. Wenzel, K. Roth, B. S. Veeling, J. Świątkowski, L. Tran, S. Mandt, J. Snoek, T. Salimans, R. Jenatton and S. Nowozin. How good is the Bayes posterior in deep neural networks really? Proceedings of ICML, 2020.
• A. G. Wilson and P. Izmailov. Bayesian deep learning and a probabilistic perspective of generalization. Advances in Neural Information Processing Systems 33, 2020.

Research directions (ch. 10)

• P. Alquier. User-friendly introduction to PAC-Bayes bounds. Foundations and Trends in Machine Learning, 17(2):174–303, 2024.
• A. N. Angelopoulos and S. Bates. Conformal prediction: a gentle introduction. Foundations and Trends in Machine Learning, 16(4):494–591, 2023.
• R. Bardenet and A. Hardy. Monte Carlo with determinantal point processes. Annals of Applied Probability, 30(1):368–417, 2020.
• P. G. Bissiri, C. C. Holmes and S. G. Walker. A general framework for updating belief distributions. Journal of the Royal Statistical Society B, 78(5):1103–1130, 2016.
• G. Cardoso, Y. Janati El Idrissi, S. Le Corff and E. Moulines. Monte Carlo guided denoising diffusion models for Bayesian linear inverse problems. International Conference on Learning Representations, 2024.
• H. Chung, J. Kim, M. T. McCann, M. L. Klasky and J. C. Ye. Diffusion posterior sampling for general noisy inverse problems. International Conference on Learning Representations, 2023.
• K. Cranmer, J. Brehmer and G. Louppe. The frontier of simulation-based inference. Proceedings of the National Academy of Sciences, 117(48):30055–30062, 2020.
• V. Fortuin. Priors in Bayesian deep learning: a review. International Statistical Review, 90(3):563–591, 2022.
• P. Grünwald and T. van Ommen. Inconsistency of Bayesian inference for misspecified linear models, and a proposal for repairing it. Bayesian Analysis, 12(4):1069–1103, 2017.
• N. Hollmann, S. Müller, K. Eggensperger and F. Hutter. TabPFN: a transformer that solves small tabular classification problems in a second. International Conference on Learning Representations, 2023.
• S. Müller, N. Hollmann, S. Pineda Arango, J. Grabocka and F. Hutter. Transformers can do Bayesian inference. International Conference on Learning Representations, 2022.
• T. Papamarkou et al. Position: Bayesian deep learning is needed in the age of large-scale AI. Proceedings of ICML, 2024.
• S. M. Xie, A. Raghunathan, P. Liang and T. Ma. An explanation of in-context learning as implicit Bayesian inference. International Conference on Learning Representations, 2022.