Bayesian reasoning

Most of machine learning, as taught in earlier courses, produces a single answer: one set of weights, one prediction, one cluster label. But data are always finite and noisy, so many answers are consistent with them. Probabilistic modelling keeps track of that uncertainty explicitly. It describes how the data could have been generated, treats the unknown quantities as random variables, and uses the rules of probability to work out what the observed data say about them. The result is not a single answer but a distribution over answers, which tells us how confident to be, how to combine evidence from different sources, and how to make decisions when being wrong has costs.

This course develops that toolkit. This lesson sets out Bayesian reasoning: Bayes' theorem, priors and posteriors, conjugate models, prediction, decisions and model comparison. Later lessons build the standard distributions and estimation methods, graphical models, latent variables, Gaussian processes, Monte Carlo sampling, Bayesian regression, variational inference, model checking, state-space models and uncertainty in machine learning; the section at the end of this lesson shows how they fit together. The course assumes the probability and statistics of STAT-252 and the machine learning foundations of AI-210.

Bayes' theorem

For a hypothesis or parameter \(\theta\) and observed data \(D\), the definition of conditional probability gives Bayes' theorem:

Each term has a name. The prior \(p(\theta)\) expresses what is believed about \(\theta\) before seeing the data. The likelihood \(p(D \mid \theta)\), viewed as a function of \(\theta\), measures how well each value explains the data. The posterior \(p(\theta \mid D)\) is the updated belief. The evidence or marginal likelihood \(p(D)\) is the normalizing constant: it does not depend on \(\theta\), so for many purposes we write simply posterior \(\propto\) likelihood \(\times\) prior.

Bayes' theorem is uncontroversial as mathematics; the Bayesian approach is distinctive in applying it to parameters, treating them as uncertain quantities with distributions rather than fixed unknown constants. The frequentist approach of classical statistics instead reasons about the distribution of estimators over hypothetical repeated data sets. Both are used throughout machine learning, and they often agree when data are plentiful; they differ most when data are scarce, when prior knowledge is strong, or when the uncertainty itself is what matters.

Priors, likelihoods and posteriors

A Bayesian model has two parts: a likelihood describing how data arise given parameters, and a prior over the parameters. Consider estimating the probability \(\theta\) that a user clicks an advertisement. If \(n\) users see it and \(k\) click, the likelihood is binomial: \(p(k \mid \theta) \propto \theta^k(1 - \theta)^{n-k}\). A natural prior is the beta distribution, \(\mathrm{Beta}(\theta \mid a, b) \propto \theta^{a-1}(1 - \theta)^{b-1}\), whose mean is \(a/(a + b)\). Multiplying gives

The posterior is again a beta distribution: the data simply add the number of successes to \(a\) and failures to \(b\). The prior parameters act like pseudo-counts, so \(\mathrm{Beta}(1, 1)\), the uniform distribution, is worth two observations, and \(\mathrm{Beta}(20, 380)\) encodes a belief that the click rate is around 5%, held with the strength of 400 observations.

from scipy import stats

a, b = 2, 38                         # prior: mean 0.05, worth 40 pseudo-observations
for n, k in [(0, 0), (20, 3), (200, 18), (2000, 170)]:
    post = stats.beta(a + k, b + n - k)
    lo, hi = post.ppf([0.025, 0.975])
    mle = f"{k / n:.3f}" if n else "  -  "
    print(f"n={n:5d} k={k:4d}  MLE={mle}  posterior mean={post.mean():.3f}  "
          f"95% interval=({lo:.3f}, {hi:.3f})")

n=    0 k=   0  MLE=  -    posterior mean=0.050  95% interval=(0.006, 0.135)
n=   20 k=   3  MLE=0.150  posterior mean=0.083  95% interval=(0.028, 0.165)
n=  200 k=  18  MLE=0.090  posterior mean=0.083  95% interval=(0.052, 0.121)
n= 2000 k= 170  MLE=0.085  posterior mean=0.084  95% interval=(0.073, 0.097)

With no data the posterior is the prior. After 20 users and 3 clicks, the maximum likelihood estimate (MLE) of 0.15 is pulled strongly toward the prior mean, and the interval is wide. As data accumulate, the posterior concentrates around the observed rate and the prior matters less and less: with 2000 observations the prior's 40 pseudo-counts are nearly irrelevant.

Conjugate priors

A prior is conjugate to a likelihood when the posterior belongs to the same family as the prior, as with the beta and binomial. Conjugacy gives closed-form updates, which is why these pairs appear constantly in practice and as building blocks of larger models:

For the normal mean with known noise variance \(\sigma^2\), the posterior after \(n\) observations with sample mean \(\bar x\) is normal with

Precisions (inverse variances) add, and the posterior mean is a weighted average of the prior mean and the sample mean, each weighted by its precision. This "precision-weighted averaging" is the essence of how Bayesian updating combines information, and it reappears in Kalman filters, linear regression and Gaussian processes.

From posterior to answers

Point estimates

A full posterior is the complete answer, but decisions often need a single number. The posterior mean minimizes expected squared error, the posterior median minimizes expected absolute error, and the posterior mode is the maximum a posteriori (MAP) estimate:

MAP estimation is maximum likelihood plus a penalty from the prior, which connects Bayesian reasoning to regularization. A zero-mean Gaussian prior on regression weights gives \(\log p(\mathbf w) = -\|\mathbf w\|^2/(2\tau^2) + \text{const}\), so MAP is ridge regression with \(\lambda = \sigma^2/\tau^2\); a Laplace prior gives the lasso. Weight decay in neural networks is a Gaussian prior in disguise. MAP, however, discards the uncertainty and depends on the parameterization, so it is a convenience rather than the Bayesian answer.

Credible intervals

A 95% credible interval contains 95% of the posterior probability, so it supports the statement "given the model and data, \(\theta\) lies in this interval with probability 0.95". This is the interpretation people usually want, and which a frequentist confidence interval does not literally provide (a confidence interval procedure covers the true value in 95% of repeated experiments). With weak priors and plenty of data the two intervals are often numerically similar.

Posterior predictive distribution

Predictions should average over parameter uncertainty rather than plug in a single estimate:

For the beta-binomial model, the probability that the next user clicks is the posterior mean \((a + k)/(a + b + n)\). With a uniform prior this is Laplace's rule of succession, \((k + 1)/(n + 2)\), which never predicts probability exactly 0 or 1 from finite data, unlike the MLE. Posterior predictive distributions are wider than plug-in predictions, correctly reflecting that the parameters are not known exactly.

Bayesian decisions: an A/B test

Posteriors answer the questions decision makers actually ask. Suppose two versions of a web page are tested: A has 120 conversions from 2400 visitors and B has 145 from 2380. Rather than a p-value, a Bayesian analysis gives the probability that B is better, and the expected gain or loss from choosing it, by sampling from the two posteriors.

import numpy as np

rng = np.random.default_rng(0)
nA, kA = 2400, 120
nB, kB = 2380, 145
pA = rng.beta(1 + kA, 1 + nA - kA, size=200_000)
pB = rng.beta(1 + kB, 1 + nB - kB, size=200_000)

lift = pB / pA - 1
print(f"P(B better than A)      = {np.mean(pB > pA):.3f}")
print(f"median relative lift     = {np.median(lift):.1%}")
print(f"95% interval for lift    = ({np.quantile(lift, 0.025):.1%}, {np.quantile(lift, 0.975):.1%})")
print(f"expected loss choosing B = {np.mean(np.maximum(pA - pB, 0)):.5f}")
print(f"expected loss choosing A = {np.mean(np.maximum(pB - pA, 0)):.5f}")

P(B better than A)      = 0.950
median relative lift     = 21.7%
95% interval for lift    = (-3.7%, 54.0%)
expected loss choosing B = 0.00014
expected loss choosing A = 0.01104

The expected loss is the conversion rate we would give up, on average, by choosing a version if it is actually the worse one. A common decision rule is to stop the test when the expected loss of the leading option falls below a threshold of caring, for example 0.0005 (0.05 percentage points); here B's expected loss of 0.00014 is already below it, even though the interval for the lift still includes zero. Bayesian analyses still require care: peeking repeatedly at the results and stopping when they look good does not invalidate the posterior itself, but it does change the long-run error rates, and the prior and model assumptions must be stated.

Model comparison and the marginal likelihood

The evidence \(p(D) = \int p(D \mid \theta)\,p(\theta)\,d\theta\), ignored when computing a posterior, becomes central when comparing models. For two models \(M_1\) and \(M_2\), Bayes' theorem at the model level gives

The ratio of marginal likelihoods is the Bayes factor. The marginal likelihood is the probability the model assigned to the observed data before seeing them, averaged over its prior. A flexible model with many parameters spreads its probability over many possible data sets, so it assigns less to any particular one; it wins only if the data really need the flexibility. This automatic Occam's razor penalizes complexity without a separate penalty term. Its drawback is sensitivity to the prior: a vague prior spreads probability thinly and is penalized heavily, even when the posterior hardly depends on it.

import numpy as np
from scipy.special import betaln

# 100 coin flips with 61 heads. M1: fair coin. M2: unknown bias with a uniform prior.
n, k = 100, 61
log_ev_fair = n * np.log(0.5)
log_ev_free = betaln(1 + k, 1 + n - k) - betaln(1, 1)      # integral of theta^k (1-theta)^(n-k)
bf = np.exp(log_ev_free - log_ev_fair)
print(f"log evidence, fair coin   = {log_ev_fair:.3f}")
print(f"log evidence, free bias   = {log_ev_free:.3f}")
print(f"Bayes factor (free : fair) = {bf:.2f}")

# For comparison: the classical two-sided binomial test
from scipy.stats import binomtest
print(f"two-sided p-value          = {binomtest(k, n, 0.5).pvalue:.4f}")

log evidence, fair coin   = -69.315
log evidence, free bias   = -68.984
Bayes factor (free : fair) = 1.39
two-sided p-value          = 0.0352

With 61 heads in 100 flips, the classical test rejects a fair coin at the 5% level, but the Bayes factor is only about 1.4, barely favouring the biased coin. The two approaches answer different questions: the p-value measures how surprising the data are under the fair model alone, while the Bayes factor compares how well the two specific models predicted the data. A uniform prior on the bias wastes probability on values like 0.05 or 0.95, which the data clearly rule out. Computing the marginal likelihood is usually hard, since it is an integral over all parameters; Lessons 7 and 9 present methods that approximate it.

Hierarchical models and shrinkage

Real problems often involve many related parameters: conversion rates for hundreds of products, test scores for many schools, failure rates for many machines of the same type. Estimating each separately (no pooling) is noisy for groups with little data; assuming they are all equal (complete pooling) ignores real differences. A hierarchical model treats the group parameters as drawn from a common population distribution whose parameters are also learned:

The posterior for each \(\theta_j\) is pulled toward the population mean, strongly for small groups and weakly for large ones. This partial pooling, or shrinkage, reduces total estimation error and is one of the most practically valuable ideas in statistics: it explains why a product with 2 sales from 3 visitors should not be ranked above one with 900 from 10,000.

How this course is organized

The twelve lessons build from foundations to applications:

• Lesson 1 (this lesson): Bayesian reasoning, conjugate models, decisions and model comparison.
• Lesson 2: standard distributions, maximum likelihood and MAP estimation, the exponential family, the multivariate Gaussian and KL divergence.
• Lesson 3: graphical models, conditional independence, exact inference and hidden Markov models.
• Lesson 4: learning graphical models, loopy belief propagation, conditional random fields and causal inference.
• Lesson 5: latent variable models and the EM algorithm.
• Lesson 6: Gaussian processes for regression and classification.
• Lesson 7: Monte Carlo sampling, MCMC and its diagnostics.
• Lesson 8: Bayesian linear regression, generalized linear models and hierarchical regression.
• Lesson 9: variational inference and variational autoencoders.
• Lesson 10: probabilistic programming, the Bayesian workflow, predictive checks and model comparison.
• Lesson 11: state-space models, Kalman filters and particle filters.
• Lesson 12: uncertainty in machine learning, calibration, conformal prediction, Bayesian optimization and a capstone project.