28 - Bayesian Statistics: Priors, Posteriors, Conjugates and Why Regularisation Is a Belief

Machine Learning Engineer course · Module 02 · Mathematics and statistics · Lesson 13

Objective. To apply Bayes' theorem not to an event but to a parameter, that is, to treat what you do not know as a distribution that the data modify. By the end you can update a belief with a single formula, derive in full the two conjugates you always meet (Beta-Binomial and Normal-Normal), recognise that the MAP estimate under a Gaussian prior is ridge regression and under a Laplace prior is the lasso, and tell a credible interval from a confidence interval without cheating on the words.

A myth to dispel. «Bayesian statistics is subjective, frequentist statistics is objective.» No: both start from declared assumptions. The difference is that the Bayesian approach forces you to write your initial belief into a formula, where it can be seen and argued with; the frequentist approach hides it in the choice of model, of test family and of significance level. A written belief is more objective than an implicit one.

Introduction

In lesson 18, on probability, we used Bayes' theorem on events: the diagnostic test that comes back positive, the door the host opens. There the unknown was «which event occurred», and the answer was a number. In lesson 20, on statistics, we worked in an entirely frequentist way: estimators, standard error, confidence intervals, with the parameter treated as an unknown but fixed constant. This lesson makes the third move, the one that connects the two: it applies Bayes to parameters. The unknown is no longer an event, it is θ, and the answer is not a number but a distribution.

This is not a stylistic exercise. Two things already written in this course change appearance when you look at them from here. The first is the +1 correction from lesson 32 on Naive Bayes: it is not a trick against zero counts, it is a declared prior. The second is the L2 penalty from lesson 44 on overfitting and regularisation, presented there as a term that «discourages large coefficients»: it is exactly the logarithm of a Gaussian prior, and its weight λ has a closed form that says how much you believe the prior relative to the data.

The end result of this lesson. You will be able to take a problem with little data, write a prior you can defend, derive the posterior by hand when the conjugate allows it and by approximation when it does not, and say precisely what you are entitled to conclude.

Roadmap of the lesson.

1. From likelihood to posterior: the four pieces of the formula, and which of the four is the problem.
2. The Beta-Binomial and Normal-Normal conjugates, derived in full, and the Laplace correction that follows.
3. MAP versus MLE: ridge and lasso unmasked as priors.
4. Credible interval versus confidence interval, measured by simulation, and the posterior predictive.
5. An A/B test updated day by day, against the same data read by a frequentist.
6. When Bayes is the wrong tool, with the criteria for noticing beforehand.

The outline follows the thread that holds the lesson together: the formula with its four pieces, two conjugates derived by hand, the discovery that ridge and lasso are priors written in another notation, then the measurements by simulation and finally the limits. If you have already seen Bayes on events you can skip to the third item, but the step from conjugates to regularisation depends on the pseudo-counts of the second.

From likelihood to posterior

The analogy: the doctor who already knows the ward

Two doctors see the same patient, with the same symptom, and measure the same laboratory value. The first works in a tropical medicine ward, where the rare disease in question turns up every week. The second works in a suburban surgery, where it has never been seen. The first considers the diagnosis likely, the second considers it unlikely and orders a second test. Both are right, and not because medicine is a matter of opinion: because they start from a different prevalence, and prevalence is information, not prejudice.

In inference terms: the laboratory value is the likelihood, identical for the two; the prevalence of the ward is the prior, different; the conclusion is the posterior, different as a result. If you replace «ward» with «application domain» and «prevalence» with «what you know about the parameter before looking at this data», you have Bayesian statistics. And you also have its price: two people with different priors, on the same data, reach different conclusions, and the only way to argue about it is for the prior to be written down.

The formula, piece by piece

Let θ be the parameter (a probability, a mean, a vector of coefficients) and D the observed data. Bayes' theorem, applied to θ, says:

p(θ | D) = p(D | θ) · p(θ) / p(D)

and, since p(D) does not depend on θ, in practice it is almost always used in the proportional form:

p(θ | D) ∝ p(D | θ) · p(θ)

Read it as a sentence: the belief afterwards is the belief beforehand, corrected by how well each value of the parameter explains the data. The four pieces — prior p(θ), likelihood p(D | θ), posterior p(θ | D) and evidence p(D) — are defined in the glossary at the end. What needs fixing here is the asymmetry almost everyone stumbles over: prior and posterior are densities over θ and integrate to 1, the likelihood does not. As a function of θ, p(D | θ) is not a probability distribution and there is no reason for its area to be worth anything.

It is the evidence that makes the general case hard. That integral sign is harmless with one parameter and devastating with a thousand: if θ has 500 components, ∫ p(D | θ) p(θ) dθ is a 500-dimensional integral that no formula closes and no grid covers. All the machinery you will meet later — grids, Monte Carlo, Markov chains, variational inference — exists to get around that integral. The conjugates are the few lucky cases in which it can be solved by hand.

The Beta-Binomial conjugate, derived in full

The smallest case that contains everything: a biased coin, a conversion rate, a failure probability. The parameter θ ∈ [0, 1] is a probability of success, the data are k successes out of n independent trials. We need a prior on θ, that is, a density on the interval [0, 1]: the Beta family does exactly this.

Beta(θ; a, b) = θ^{a−1} (1 − θ)^{b−1} / B(a, b), with B(a, b) the normalising constant.

The binomial likelihood, as a function of θ, is

p(D | θ) = C(n, k) · θ^k (1 − θ)^{n−k}

where C(n, k) is the binomial coefficient, which does not depend on θ and will therefore end up in the constant. Let us multiply prior and likelihood, keeping only the factors that depend on θ:

p(θ | D) ∝ θ^k (1 − θ)^{n−k} · θ^{a−1} (1 − θ)^{b−1}

Powers with the same base add in the exponents:

p(θ | D) ∝ θ^{a + k − 1} (1 − θ)^{b + n − k − 1}

Now the decisive move, and it is worth saying exactly what it is: recognising the shape. That expression is θ raised to «something minus one» times (1 − θ) raised to «something else minus one»: it is the kernel of a Beta. Since there is only one density on [0, 1] proportional to that kernel, the normalising constant is forced and there is no need to compute it. Therefore

p(θ | D) = Beta(a + k, b + n − k) ∎

This is the operational definition of a conjugate: a family of priors such that, multiplied by that likelihood, it returns a prior of the same family. The evidence integral was not computed: it was avoided, because we knew in advance which density we would obtain.

The reading that matters: a and b are pseudo-counts

Look at the update formula: a receives the successes, b the failures. So a and b read as fake data you have already seen: a Beta(2, 2) prior declares «before starting, I behave as if I had already observed 2 successes and 2 failures». This makes the prior arguable in concrete terms: not «how confident am I», but «how many observations is my conviction worth». That is the right question to ask in a meeting, and it is answered with a number.

import numpy as np
from scipy import stats

a, b = 2.0, 2.0
n, k = 10, 7
a_post, b_post = a + k, b + n - k
post = stats.beta(a_post, b_post)
print("prior      : Beta(%.0f, %.0f)   mean = %.4f" % (a, b, a / (a + b)))
print("data       : %d successes out of %d   frequency = %.4f" % (k, n, k / n))
print("posterior  : Beta(%.0f, %.0f)   mean = %.4f" % (a_post, b_post, post.mean()))
print("posterior mode (a-1)/(a+b-2) =", round((a_post - 1) / (a_post + b_post - 2), 4))
print("maximum likelihood estimate k/n =", round(k / n, 4))
print("95%% credible interval            : [%.4f, %.4f]" % post.interval(0.95))
print("prior pseudo-counts: %.0f successes and %.0f failures already seen" % (a, b))

Output obtained (Python 3.14.4, SciPy 1.16.3):

prior      : Beta(2, 2)   mean = 0.5000
data       : 7 successes out of 10   frequency = 0.7000
posterior  : Beta(9, 5)   mean = 0.6429
posterior mode (a-1)/(a+b-2) = 0.6667
maximum likelihood estimate k/n = 0.7
95% credible interval            : [0.3857, 0.8614]
prior pseudo-counts: 2 successes and 2 failures already seen

What the code does. It applies the formula just derived — Beta(2, 2) plus 7 successes and 3 failures gives Beta(9, 5) — and prints three summaries of the same posterior: mean 9/14, mode 8/12 and the equal-tailed 95% credible interval.

What it tells us. The posterior mean 0.6429 lies between the prior mean (0.5) and the observed frequency (0.7), closer to the second because the data are worth 10 observations against the prior's 4. Mean and mode do not coincide (0.6429 against 0.6667), because the Beta is skewed when a ≠ b: choosing one or the other is a decision, not a detail. And the interval [0.386, 0.861] is very wide, the most useful message in the output: with ten trials you know almost nothing, and the posterior has the decency to say so.

Where the Naive Bayes +1 correction comes from

In lesson 32 Naive Bayes estimated the probability of each word in each class as count / total, and broke on a word never seen: a zero factor zeroed the product and with it the class. The remedy was to add 1 to all counts. Seen from here, that correction has a name and a derivation: the case with K outcomes instead of 2 has the Dirichlet as its conjugate, which stands to the multinomial as the Beta stands to the binomial, and its posterior is Dirichlet(α₁ + c₁, …, α_K + c_K). With α = (1, …, 1) — the uniform prior on the simplex — the posterior mean of category j is

(c_j + 1) / (N + K)

which is, character for character, the Laplace correction. It is not a plug against a division by zero: it is the mean of a posterior under the prior «every outcome is possible and I have seen one fake instance of each».

conteggi = np.array([12, 5, 0, 3])          # counts of 4 words in one class
K = len(conteggi)
N = conteggi.sum()
mle = conteggi / N
media_post = (conteggi + 1) / (N + K)       # mean of Dirichlet(1+c1, ..., 1+cK)
print("observed counts           :", conteggi)
print("MLE estimate (c/N)        :", np.round(mle, 4))
print("posterior mean (c+1)/(N+K):", np.round(media_post, 4))
print("sum of the means          :", round(float(media_post.sum()), 10))
print("the never-seen word goes from", mle[2], "to", round(float(media_post[2]), 4))

Output obtained:

observed counts           : [12  5  0  3]
MLE estimate (c/N)        : [0.6  0.25 0.   0.15]
posterior mean (c+1)/(N+K): [0.5417 0.25   0.0417 0.1667]
sum of the means          : 1.0
the never-seen word goes from 0.0 to 0.0417

What the code does. It takes four counts, one of which is zero, and compares the maximum likelihood estimate c/N with the posterior mean (c+1)/(N+K). The sum line checks that the posterior means are still a distribution: exactly 1.0.

What it tells us. The word never seen goes from 0.0 to 0.0417: it is no longer impossible, merely improbable. The price is visible in the first column, which drops from 0.60 to 0.5417: the probability given to the absent outcome is subtracted from the others, because the total stays 1. Anyone using the correction without knowing this is making that trade blind.

The rule to take away. Every time you «add a small number for stability» you are choosing a prior. The only difference between you and a Bayesian is whether you wrote it down.

How the posterior tightens

With a = b = 2 the pseudo-counts are 4: they are 29% of the total at n = 10, 3.8% at n = 100, 0.4% at n = 1000. The standard deviation of the posterior, meanwhile, falls as 1/√n.

The Normal-Normal conjugate, derived in full

The second conjugate worth knowing concerns a mean. We observe n values x₁, …, x_n from a normal with unknown mean μ and known variance σ², and we put on μ a normal prior N(μ₀, τ²). The known-variance assumption is a real and declared simplification: it is what gets the formula out in one page. The case with unknown variance as well exists, has the Normal-Inverse-Gamma as its conjugate, and is in Gelman ch. 3.

Let us write the logarithm of the product prior × likelihood, throwing away everything that does not depend on μ:

log p(μ | D) = −(1/2σ²) · Σᵢ (xᵢ − μ)² − (1/2τ²) · (μ − μ₀)² + const.

Let us expand the two squares and collect by powers of μ. From the first term: Σᵢ (xᵢ − μ)² = Σᵢ xᵢ² − 2μ Σᵢ xᵢ + n μ². From the second: (μ − μ₀)² = μ² − 2μμ₀ + μ₀². Putting them together and discarding the constants:

log p(μ | D) = −(1/2) · [ μ² (n/σ² + 1/τ²) − 2μ (n x̄/σ² + μ₀/τ²) ] + const.

where x̄ is the sample mean. Let us introduce two names that make everything readable: the precision is the inverse of the variance. Call λ₀ = 1/τ² the precision of the prior and λ_D = n/σ² the precision of the data, and set λ_n = λ₀ + λ_D. Then

log p(μ | D) = −(λ_n/2) · [ μ² − 2μ · (λ_D x̄ + λ₀ μ₀)/λ_n ] + const.

Now the completion of the square, worked out: to make the bracketed expression a perfect square (μ − m)² we need m = (λ_D x̄ + λ₀ μ₀)/λ_n, because (μ − m)² = μ² − 2μm + m² and the term m² is a constant we may add and subtract freely:

log p(μ | D) = −(λ_n/2) · (μ − m)² + const.

Exponentiating, p(μ | D) ∝ exp[ −(μ − m)²/(2/λ_n) ], which is the kernel of a normal with mean m and variance 1/λ_n. Therefore

p(μ | D) = N( (λ₀ μ₀ + λ_D x̄)/(λ₀ + λ_D) , 1/(λ₀ + λ_D) ) ∎

Read it as two sentences, because they are the two things to remember from the whole section. The posterior mean is the weighted average of prior mean and data mean, with weights equal to their respective precisions. And the posterior precision is the sum of the precisions: information adds up, and it adds up in the right coordinate, which is the inverse of the variance, not the variance.

The three conjugates met so far all read with the same key, and the right-hand column is the one you take into a meeting.

How to read the last column. The weight of the data grows linearly in n, that of the prior stays put: this is why the prior fades. And the number in the third row, σ²/τ², is no notational coincidence: it is exactly the λ of the ridge we are about to derive. The regularisation coefficient is the number of observations your belief is worth.

mu0, tau2 = 100.0, 25.0                     # prior: N(100, 25)
sigma2 = 16.0                               # known variance of the data
rng = np.random.default_rng(42)
dati = rng.normal(108.0, np.sqrt(sigma2), size=8)
n_oss = dati.size
prec_prior = 1.0 / tau2
prec_dati = n_oss / sigma2
mu_post = (prec_prior * mu0 + prec_dati * dati.mean()) / (prec_prior + prec_dati)
var_post = 1.0 / (prec_prior + prec_dati)
print("sample mean              :", round(float(dati.mean()), 4))
print("prior precision 1/t2     :", round(prec_prior, 6))
print("data precision n/s2      :", round(prec_dati, 6))
print("posterior mean           :", round(float(mu_post), 4))
print("posterior variance       :", round(float(var_post), 6))

# independent check: the same posterior computed by brute force on a grid
griglia = np.linspace(80, 130, 200001)
logpost = stats.norm.logpdf(griglia, mu0, np.sqrt(tau2))
logpost = logpost + stats.norm.logpdf(dati[:, None], griglia[None, :], np.sqrt(sigma2)).sum(axis=0)
dens = np.exp(logpost - logpost.max())      # subtract the maximum: avoids underflow
dens /= np.trapezoid(dens, griglia)
media_griglia = np.trapezoid(griglia * dens, griglia)
var_griglia = np.trapezoid((griglia - media_griglia) ** 2 * dens, griglia)
print("grid mean                :", round(float(media_griglia), 4))
print("grid variance            :", round(float(var_griglia), 6))
print("maximum difference       : %.2e" % max(abs(media_griglia - mu_post), abs(var_griglia - var_post)))

Output obtained (seed 42, NumPy 2.3.5):

sample mean              : 106.7571
prior precision 1/t2     : 0.04
data precision n/s2      : 0.5
posterior mean           : 106.2565
posterior variance       : 1.851852
grid mean                : 106.2565
grid variance            : 1.851852
maximum difference       : 1.42e-14

What the code does. The first half applies the closed form, the second checks it by brute force on a grid of 200,001 points. The detail worth stealing is dens = np.exp(logpost - logpost.max()): summing eight log-densities gets you to numbers like −10³, and exp(−1000) is zero in floating point; subtracting the maximum before exponentiating does not change the result, because it is a constant factor that normalisation removes.

What it tells us. The two routes agree to 1.42e-14, that is, to the last digit of double precision: the derivation is right. On the substance: the precision of the data is 0.5 against the prior's 0.04, so the data weigh twelve and a half times as much, and indeed the posterior mean 106.2565 is far closer to 106.7571 (the data) than to 100 (the prior). The posterior variance 1.85 is smaller than both τ² = 25 and σ²/n = 2: combining two pieces of information gives something more precise than either, and this is the operational content of the sum of precisions.

MAP versus MLE, and regularisation unmasked

A posterior is a distribution, but sometimes you need one number: to write a prediction into a system, to compare with a frequentist estimate, to initialise an optimisation. The most used summary is the point of maximum posterior density, the MAP estimate (maximum a posteriori). Since the logarithm is increasing and the evidence does not depend on θ:

θ_MAP = argmax_θ p(θ | D) = argmax_θ [ log p(D | θ) + log p(θ) ]

Compare it with the maximum likelihood estimate, θ_MLE = argmax_θ log p(D | θ). The difference is a single addend: the logarithm of the prior. And here is where the interesting thing of the whole lesson happens, because that addend, for the commonest priors, is an object you have already seen under another name.

Gaussian prior: the term that appears is ridge

Take a linear regression with Gaussian noise of known variance σ²: y = Xβ + ε, ε ~ N(0, σ²). The log-likelihood, up to constants, is −‖y − Xβ‖²/(2σ²). Put on each coefficient an independent Gaussian prior centred at zero, β_j ~ N(0, τ²): its logarithm, up to constants, is −‖β‖²/(2τ²). Adding:

log p(β | D) = −‖y − Xβ‖²/(2σ²) − ‖β‖²/(2τ²) + const.

Maximising is the same as minimising the opposite; multiplying by 2σ², which is a positive constant and does not move the minimum:

β_MAP = argmin_β [ ‖y − Xβ‖² + (σ²/τ²) · ‖β‖² ]

This is, character for character, the ridge regression of lesson 44, with

λ = σ²/τ²

Read the formula for λ, because it says everything lesson 44 could not say. The numerator is the noise in the data, the denominator is the width of your belief about the coefficients. You penalise heavily when the data are noisy (σ² large) or when you believe the coefficients are small (τ² small). And the case λ = 0 is not «no regularisation»: it is τ² → ∞, that is, the prior declaring every value of the coefficients equally likely, including 10¹². That is not the absence of a belief, it is a rather strange belief.

from scipy import optimize
from sklearn.linear_model import Ridge

rng = np.random.default_rng(0)
n_righe, d = 60, 5
X = rng.normal(size=(n_righe, d))
beta_vero = np.array([2.0, -1.5, 0.0, 0.8, -0.3])
sigma2_r, tau2_r = 1.0, 0.5
y = X @ beta_vero + rng.normal(0, np.sqrt(sigma2_r), size=n_righe)
lam = sigma2_r / tau2_r


def neg_log_post(beta):
    """Minus the logarithm of the posterior, up to constants."""
    res = y - X @ beta
    return res @ res / (2 * sigma2_r) + beta @ beta / (2 * tau2_r)


def grad(beta):
    """Analytic gradient: needed to reach a tight tolerance."""
    return -X.T @ (y - X @ beta) / sigma2_r + beta / tau2_r


sol = optimize.minimize(neg_log_post, np.zeros(d), jac=grad, method="BFGS",
                        options={"gtol": 1e-12})
ridge = Ridge(alpha=lam, fit_intercept=False).fit(X, y)
print("lambda = sigma2/tau2 =", lam)
print("MAP  (scipy.optimize):", np.round(sol.x, 6))
print("Ridge(alpha=lambda)  :", np.round(ridge.coef_, 6))
print("maximum difference   : %.2e" % np.abs(sol.x - ridge.coef_).max())
print("MLE (lambda = 0)     :", np.round(np.linalg.lstsq(X, y, rcond=None)[0], 6))

Output obtained (seed 0, scikit-learn 1.8.0):

lambda = sigma2/tau2 = 2.0
MAP  (scipy.optimize): [ 2.014519 -1.370516  0.105054  0.737658 -0.194408]
Ridge(alpha=lambda)  : [ 2.014519 -1.370516  0.105054  0.737658 -0.194408]
maximum difference   : 1.33e-15
MLE (lambda = 0)     : [ 2.147018 -1.448232  0.118071  0.794285 -0.202626]

What the code does. It builds a synthetic 60×5 regression with known true coefficients (data simulated by us, declared as such), computes the MAP by minimising −log p(β | D) with an analytic gradient, and compares it with Ridge(alpha=σ²/τ²). Two details that are not decorative: fit_intercept=False, because the prior is on the coefficients and not on an intercept, and gtol=1e-12, because with the default tolerance the difference would stay around 1e-6, enough to cast doubt on an identity that is in fact exact.

What it tells us. The two vectors agree to 1.33e-15: they are not «similar», they are the same thing computed in two ways. The last line adds the useful comparison: the MLE is larger in absolute value on every component (2.147 against 2.015). That is the shrinkage towards zero imposed by the prior, and on the third component — true coefficient 0 — it takes the estimate from 0.118 to 0.105: in the right direction, but without reaching zero. The reason is in the next subsection.

Laplace prior: the term that appears is the lasso

Let us change the prior and leave everything else identical. In place of the Gaussian we put on each coefficient a Laplace distribution centred at zero, p(β_j) ∝ exp(−|β_j| / b): same idea («the coefficients are near zero»), but with a sharp peak at the origin and heavier tails. Its logarithm, up to constants, is −‖β‖₁ / b, where ‖β‖₁ = Σ_j |β_j|. Repeating the earlier steps:

β_MAP = argmin_β [ ‖y − Xβ‖² + (2σ²/b) · ‖β‖₁ ]

which is the lasso, with penalty coefficient 2σ²/b. Before verifying it there is a conversion trap to defuse, because it makes an exact identity look wrong. sklearn.linear_model.Lasso does not minimise that quantity: it minimises (1/2n)·‖y − Xβ‖² + α·‖β‖₁, that is, the same thing divided by 2n. Equating the two ratios of penalty to data term gives 2nα = 2σ²/b, that is

α = σ² / (n · b)

Anyone skipping this step compares two differently regularised models, sees different coefficients and concludes that the correspondence does not hold. It does hold: it is the conversion that is wrong.

To see the difference that matters we need a design in which there is something to zero out: 80 rows, 12 columns of which 8 have a true coefficient of exactly zero.

from sklearn.linear_model import Lasso

rng = np.random.default_rng(3)
n_righe, d = 80, 12
X = rng.normal(size=(n_righe, d))
beta_vero = np.array([2.5, -1.8, 1.2, 0.9, 0, 0, 0, 0, 0, 0, 0, 0], dtype=float)
sigma2 = 1.0
y = X @ beta_vero + rng.normal(0, np.sqrt(sigma2), size=n_righe)

tau2 = 0.5                      # Gaussian prior  N(0, tau2)       -> ridge
scala_b = 0.05                  # Laplace prior   Laplace(0, b)    -> lasso
alpha_ridge = sigma2 / tau2                      # Ridge: ||y-Xb||^2 + alpha*||b||^2
alpha_lasso = sigma2 / (n_righe * scala_b)       # Lasso: (1/2n)||y-Xb||^2 + alpha*||b||_1

ridge = Ridge(alpha=alpha_ridge, fit_intercept=False).fit(X, y)
lasso = Lasso(alpha=alpha_lasso, fit_intercept=False, max_iter=200000).fit(X, y)
print("true coefficients :", np.round(beta_vero, 3))
print("ridge             :", np.round(ridge.coef_, 3))
print("lasso             :", np.round(lasso.coef_, 3))
print("exact zeros ridge :", int((ridge.coef_ == 0).sum()), "out of", d)
print("exact zeros lasso :", int((lasso.coef_ == 0).sum()), "out of", d)


# Is the lasso really the MAP? Evaluate the objective at 20000 random
# perturbations: if it is the minimum, none of them should do better.
def obiettivo_map(beta):
    res = y - X @ beta
    return res @ res / (2 * sigma2) + np.abs(beta).sum() / scala_b


perturbato = lasso.coef_ + rng.normal(0, 0.02, size=(20000, d))
val_perturbati = np.array([obiettivo_map(b) for b in perturbato])
print("ridge: smallest coefficient in absolute value = %.6f" % np.abs(ridge.coef_).min())
print("MAP objective at the lasso       : %.6f" % obiettivo_map(lasso.coef_))
print("MAP objective at the ridge       : %.6f" % obiettivo_map(ridge.coef_))
print("best of 20000 perturbations      : %.6f" % val_perturbati.min())
print("perturbations beating lasso     :",
      int((val_perturbati < obiettivo_map(lasso.coef_) - 1e-9).sum()))

# Optimality condition: |X_j^T r| / n <= alpha for every switched-off coefficient.
residuo = y - X @ lasso.coef_
correlazioni = np.abs(X.T @ residuo) / n_righe
spenti = lasso.coef_ == 0
print("coefficients off: max |X_j^T r|/n = %.6f  (threshold alpha = %.6f)"
      % (correlazioni[spenti].max(), alpha_lasso))
print("coefficients on : min |X_j^T r|/n = %.6f  (should be ~alpha)"
      % correlazioni[~spenti].min())

Output obtained (seed 3, scikit-learn 1.8.0):

true coefficients : [ 2.5 -1.8  1.2  0.9  0.   0.   0.   0.   0.   0.   0.   0. ]
ridge             : [ 2.006 -1.73   1.276  0.829  0.044 -0.148 -0.144  0.006  0.073  0.182
  0.016 -0.016]
lasso             : [ 1.867 -1.602  1.119  0.582  0.    -0.    -0.     0.     0.     0.008
  0.    -0.   ]
exact zeros ridge : 0 out of 12
exact zeros lasso : 7 out of 12
ridge: smallest coefficient in absolute value = 0.005879
MAP objective at the lasso       : 148.465847
MAP objective at the ridge       : 162.941647
best of 20000 perturbations      : 148.876569
perturbations beating lasso     : 0
coefficients off: max |X_j^T r|/n = 0.165989  (threshold alpha = 0.250000)
coefficients on : min |X_j^T r|/n = 0.249988  (should be ~alpha)

What the code does. It fits ridge and lasso to the same design with the two conversions just derived, counts the exact zeros (an == 0 comparison, not «smaller than a threshold»), and checks in two independent ways that the lasso is the MAP: it evaluates the objective at 20,000 vectors obtained by perturbing the solution, and it checks the optimality condition, by which a coefficient can sit at zero only if its correlation with the residual |Xⱼᵀr|/n does not exceed α.

What it tells us. The lasso zeroes exactly 7 coefficients out of 12, ridge none: its smallest is 0.005879, and it will stay so for any finite λ. The MAP check is clear-cut: objective 148.47 against ridge's 162.94, and zero perturbations out of 20,000 improve on it. The last pair of lines is the geometric explanation in two numbers: the switched-off coefficients have correlation with the residual of at most 0.166, below the 0.25 threshold; the switched-on ones have it at 0.249988, resting against the threshold. The lasso switches off everything that does not clear α, and keeps the rest just above it.

The lasso is not free. Of the 8 genuinely useless columns it zeroes 7: the ninth stays at 0.008, switched on by chance. And it shrinks the true coefficients more than ridge does — 0.582 against 0.829 where the truth is 0.9: that bias on the survivors is the price of selection, and in prediction it can cost more than sparsity gains. Why one zeroes and the other does not: the derivative of β² at zero is zero, so near the origin the L2 penalty stops pushing; the derivative of |β| at zero does not exist — there is a kink — and the push is ±1 right to the end, until the data resist with force at least α. The kink in the Laplace prior is variable selection.

Credible interval versus confidence interval

The two definitions, side by side, with a difference that is not a shade of language.

95% credible interval. An interval [l, u] such that ∫_l^u p(θ | D) dθ = 0.95: given this posterior, 95% of my belief about the parameter lies in there. It is a statement about θ, and it makes sense to speak of the probability that θ is inside.

95% confidence interval. A procedure which, applied to data generated repeatedly from the model, produces intervals containing the true value in 95% of cases: I know nothing about this particular interval; I know the method is wrong one time in twenty. It is a statement about the procedure, and θ is a constant: it is not legitimate to say «there is a 95% probability that θ is in here».

The sentence that separates them, worth memorising: the credible one is a probability about the parameter, the confidence one is a property of the procedure. Calling a credible interval a «confidence interval» is not terminological pedantry: it promises a coverage guarantee the former does not offer, as we are about to measure.

First fact: on plentiful data they are indistinguishable

With a uniform prior and a sample of reasonable size, the two procedures give almost the same numbers. Let us measure it: 10,000 binomial experiments with n = 40 and p = 0.30, and for each we build the 95% credible interval from Beta(1 + k, 1 + n − k) and the Wilson confidence interval, which is the correct way to build an interval for a proportion (the naive normal approximation p̂ ± 1.96·√(p̂(1−p̂)/n) is notoriously poor on extreme proportions).

def wilson(k, n, z=1.959963984540054):
    """Wilson confidence interval for a proportion."""
    centro = (k + z * z / 2) / (n + z * z)
    mezzo = z / (n + z * z) * np.sqrt(k * (n - k) / n + z * z / 4)
    return centro - mezzo, centro + mezzo


rng = np.random.default_rng(11)
prove, n_oss, p_vero = 10000, 40, 0.30
ks = rng.binomial(n_oss, p_vero, size=prove)       # 10,000 experiments in one go
cred = stats.beta(1 + ks, 1 + n_oss - ks).interval(0.95)
conf = wilson(ks, n_oss)
cop_cred = float(((cred[0] <= p_vero) & (p_vero <= cred[1])).mean())
cop_conf = float(((conf[0] <= p_vero) & (p_vero <= conf[1])).mean())
print("simulated experiments         :", prove)
print("credible interval coverage  : %.4f" % cop_cred)
print("Wilson interval coverage      : %.4f" % cop_conf)
print("mean credible width           : %.4f" % float((cred[1] - cred[0]).mean()))
print("mean Wilson width             : %.4f" % float((conf[1] - conf[0]).mean()))

Output obtained (seed 11, 10,000 experiments):

simulated experiments         : 10000
credible interval coverage  : 0.9478
Wilson interval coverage      : 0.9478
mean credible width           : 0.2706
mean Wilson width             : 0.2697

What the code does. It simulates the 10,000 experiments with a single vectorised call, builds both intervals in parallel — stats.beta accepts arrays of parameters — and computes the fraction containing 0.30.

What it tells us. The two coverages are identical to four decimals, 0.9478, and the mean widths differ by 0.0009. The result should be stated for what it is: with n = 40 and a uniform prior, choosing one procedure over the other changes nothing measurable. Anyone presenting Bayesian statistics as intrinsically more prudent or more powerful has no argument here. Note also that 0.9478 is not 0.95: on a discrete variable exact coverage is unattainable for either method, and it oscillates around the nominal value as n and p change.

Second fact: on little data the prior shows through

The real difference appears where the information in the data is scarce. Let us repeat with n = 5, and add an informative and wrong prior: Beta(8, 2), whose mean is 0.8 while the truth is 0.30. The comparison is three-way: wrong prior, flat prior, Wilson.

rng = np.random.default_rng(12)
n_pic = 5
ks = rng.binomial(n_pic, p_vero, size=prove)
cred_inf = stats.beta(8 + ks, 2 + n_pic - ks).interval(0.95)    # wrong prior
cred_pia = stats.beta(1 + ks, 1 + n_pic - ks).interval(0.95)    # flat prior
conf = wilson(ks, n_pic)
for nome, iv in (("credible, Beta(8,2)    ", cred_inf),
                 ("credible, Beta(1,1)    ", cred_pia),
                 ("Wilson                 ", conf)):
    cop = float(((iv[0] <= p_vero) & (p_vero <= iv[1])).mean())
    amp = float((iv[1] - iv[0]).mean())
    print("%s: coverage %.4f | mean width %.4f" % (nome, cop, amp))

Output obtained (seed 12, 10,000 experiments, n = 5):

credible, Beta(8,2)    : coverage 0.1649 | mean width 0.4584
credible, Beta(1,1)    : coverage 0.9694 | mean width 0.6011
Wilson                 : coverage 0.9694 | mean width 0.5905

What the code does. Identical to the previous one, with n cut to 5 and three intervals instead of two. The only substantive difference is the pair of prior parameters on the first line.

What it tells us. Here the fracture shows. The flat prior and Wilson have the same coverage, 0.9694, while the informative prior drops to 0.1649: the credible interval built with Beta(8, 2) contains the true value in fewer than one experiment in six. And look at the second column: it is also the narrowest of the three, 0.4584 against 0.59-0.60. An interval that is narrower and much more wrong is the worst possible combination, because it presents itself as more informative.

It must be said precisely what this is evidence of. It is not a flaw in the method: the credible interval faithfully reports 95% of a belief that includes 10 pseudo-observations centred on 0.8, and with n = 5 those pseudo-counts are worth twice the data. The 95% frequentist coverage is a guarantee the Bayesian method offers only on average over the prior: if p really were drawn from Beta(8, 2), it would return to nominal. The collapse to 0.1649 measures how much a prior you cannot defend costs.

The rule to take away. With plenty of data the choice of statistical school is a matter of taste and tooling. With little data it is the most important choice you make, and the prior must be negotiated before looking at the results, not after.

Posterior predictive

So far we have talked about the parameter. But almost always the real question is not «what is θ», it is «what will happen next time». The correct answer is not obtained by fixing θ at its best value and using that: it is obtained by averaging the predictions over the whole posterior. It is called the posterior predictive:

p(ỹ | D) = ∫ p(ỹ | θ) · p(θ | D) dθ

In words: for each plausible value of the parameter ask what it predicts, and weight its prediction by how plausible that value is. The comparison is with the plug-in prediction, p(ỹ | θ_MAP), which is what anyone does who takes a point estimate and pushes it into the model. The two coincide only if the posterior is concentrated at a point; otherwise the plug-in underestimates the uncertainty, because it throws away the variance of the parameter and keeps only that of the noise.

In the Beta-Binomial case the integral can be solved and gives a distribution with a name: the Beta-Binomial, available in SciPy as stats.betabinom. Let us test it where the difference shows, that is, with very little data: 4 successes out of 5, and the question «out of 10 future trials, how many successes?».

n_vis, k_vis, m = 5, 4, 10
ap, bp = 1 + k_vis, 1 + n_vis - k_vis          # posterior with a flat prior
pred = stats.betabinom(m, ap, bp)              # posterior predictive, closed form
p_map = (ap - 1) / (ap + bp - 2)               # mode of the posterior
plug = stats.binom(m, p_map)                   # plug-in prediction with the MAP
print("posterior: Beta(%d, %d)   MAP = %.4f" % (ap, bp, p_map))
print("predictive Beta-Binomial : mean %.4f  std %.4f  95%% interval [%d, %d]"
      % (pred.mean(), pred.std(), pred.ppf(0.025), pred.ppf(0.975)))
print("plug-in Binomial(MAP)    : mean %.4f  std %.4f  95%% interval [%d, %d]"
      % (plug.mean(), plug.std(), plug.ppf(0.025), plug.ppf(0.975)))

Output obtained:

posterior: Beta(5, 2)   MAP = 0.8000
predictive Beta-Binomial : mean 7.1429  std 2.0825  95% interval [3, 10]
plug-in Binomial(MAP)    : mean 8.0000  std 1.2649  95% interval [5, 10]

What the code does. It builds the posterior Beta(5, 2) from 4 successes out of 5 with a flat prior, then two predictive distributions for 10 future trials: the correct one, which integrates over the whole posterior, and the plug-in, which freezes θ at the mode 0.8. For both it prints mean, standard deviation and 95% interval.

What it tells us. The correct predictive standard deviation is 2.08, the plug-in one 1.26: the latter declares an uncertainty 39% lower than it is entitled to. And the 95% interval goes from [3, 10] to [5, 10]: the plug-in rules out the outcomes 3 and 4, perfectly possible if θ is 0.55 rather than 0.8, and with five observations 0.55 is far from ruled out. The mean shifts too, from 8.0 to 7.14, because the posterior Beta(5, 2) has mean 5/7 ≈ 0.714, lower than the mode.

Trade-off. The correct predictive costs an integral and returns wider intervals, which in a presentation look like a worse result: that is why almost nobody uses it, and it is an incentive failure — the narrow interval is not more precise, it is merely more self-assured. With large n the difference vanishes and the plug-in becomes legitimate; the threshold is checked by comparing the two standard deviations, as above.

The «uninformative» prior does not exist

At this point the temptation is obvious: if the prior is the problem, let us use one that says nothing. The natural candidate is the uniform prior, Beta(1, 1) on probabilities, «flat» on the interval. The objection, which goes back to Jeffreys, is that flat is not a property of the parameter but of the coordinates you write it in. If you are indifferent among all values of θ, you cannot simultaneously be indifferent among all values of log θ or of the logit φ = log(θ/(1−θ)): the change-of-variable rule for densities, p_φ(φ) = p_θ(θ) · |dθ/dφ|, forbids it.

The prior that is invariant exists and is called the Jeffreys prior: it is proportional to the square root of the Fisher information, p(θ) ∝ √I(θ), and for the binomial it gives Beta(1/2, 1/2). Let us check both claims numerically.

theta = np.linspace(1e-6, 1 - 1e-6, 200001)
phi = np.log(theta / (1 - theta))                      # logit
piatto_theta = np.ones_like(theta)
indotta_su_phi = piatto_theta * theta * (1 - theta)    # multiply by |dtheta/dphi|
jeff_theta = 1 / np.sqrt(theta * (1 - theta))          # Jeffreys in theta: Beta(1/2, 1/2)
jeff_su_phi = jeff_theta * theta * (1 - theta)         # transformed into phi
jeff_diretta = np.sqrt(theta * (1 - theta))            # Jeffreys computed in phi
rapporto = jeff_su_phi / jeff_diretta
print("flat prior on theta, induced density on phi at the points phi = -2, 0, 2:")
for punto in (-2.0, 0.0, 2.0):
    i = int(np.abs(phi - punto).argmin())
    print("  phi = %5.2f -> density %.6f" % (phi[i], indotta_su_phi[i]))
print("ratio transformed Jeffreys / direct Jeffreys: min %.10f  max %.10f"
      % (rapporto.min(), rapporto.max()))

Output obtained:

flat prior on theta, induced density on phi at the points phi = -2, 0, 2:
  phi = -2.00 -> density 0.104992
  phi = -0.00 -> density 0.250000
  phi =  2.00 -> density 0.104992
ratio transformed Jeffreys / direct Jeffreys: min 1.0000000000  max 1.0000000000

What the code does. It transforms the uniform prior into logit coordinates by multiplying by |dθ/dφ| = θ(1−θ) and evaluates its density at three points. It then applies the same transformation to Jeffreys and compares it with Jeffreys computed directly in φ, where the Fisher information is θ(1−θ).

What it tells us. The density induced by the uniform prior on the logit is 0.25 at the centre and 0.105 at φ = ±2: it is not flat, it is a bell favouring central values of θ. Anyone declaring «I use an uninformative prior on the probability» is, in logit coordinates — the ones any logistic regression works in — declaring an informative belief centred in the middle. Jeffreys, by contrast, passes the test: ratio 1.0000000000 at both the minimum and the maximum over 200,001 points, that is, exactly invariant.

Jeffreys is not a universal solution, though: in several dimensions it can be improper and lead to improper posteriors. The rule to take away: do not ask «is this prior informative?», ask «informative with respect to which parametrisation?».

The experiment: an A/B test updated day by day

Let us put everything to work on a realistic case. Two versions of a page, true conversion rates 0.040 and 0.046 — that is, B really is better, by 15% relative, which is a considerable improvement for a conversion. The data are simulated by us with default_rng(7), so we know the truth and can judge the conclusions. We observe the stream at four checkpoints — after 1, 10, 100 and 1000 observations per arm — and at each one we compute two things on the same data: the Bayesian quantity P(θ_B > θ_A | D), estimated by sampling from the two posteriors, and the p-value of the frequentist two-proportion test, which is the tool from lesson 6 on controlled experiments.

rng = np.random.default_rng(7)
p_a, p_b = 0.040, 0.046                    # TRUE rates: we know them because we simulate
N = 1000
osservazioni_a = rng.random(N) < p_a
osservazioni_b = rng.random(N) < p_b
campionatore = np.random.default_rng(99)

print(" n | k_A | k_B | P(B > A) | frequentist p-value")
for t in (1, 10, 100, 1000):
    ka, kb = int(osservazioni_a[:t].sum()), int(osservazioni_b[:t].sum())
    sa = campionatore.beta(1 + ka, 1 + t - ka, size=200000)      # posterior of A
    sb = campionatore.beta(1 + kb, 1 + t - kb, size=200000)      # posterior of B
    prob = float((sb > sa).mean())                               # P(B > A) by sampling
    pooled = (ka + kb) / (2 * t)
    se = np.sqrt(pooled * (1 - pooled) * 2 / t)
    if se == 0:
        pval = float("nan")                                      # no conversions: test undefined
    else:
        z = (kb / t - ka / t) / se
        pval = float(2 * stats.norm.sf(abs(z)))
    print("%4d | %3d | %3d |  %.4f  | %s" % (t, ka, kb, prob,
                                             "undefined" if np.isnan(pval) else "%.4f" % pval))

Output obtained (seed 7 for the data, seed 99 for the sampling):

 n | k_A | k_B | P(B > A) | frequentist p-value
   1 |   0 |   0 |  0.5004  | undefined
  10 |   1 |   0 |  0.2385  | 0.3049
 100 |   7 |   5 |  0.2841  | 0.5515
1000 |  44 |  37 |  0.2148  | 0.4272

What the code does. It generates the two streams once and reads progressively longer prefixes of them, like an experiment accumulating traffic. P(B > A) has no convenient closed form: it is estimated by sampling 200,000 values from each posterior and counting how often B exceeds A. It is the first use of sampling as an estimator rather than as a check, and it is the thread the next lesson picks up.

What it tells us, and it is not what we expected. B has the higher true rate, and yet P(B > A) falls: from 0.50 to 0.21. At 1000 observations per arm we have seen 44 conversions for A and 37 for B, that is 4.4% against 3.7%: the worse variant is winning. This is not a bug in the code and it is the most useful thing that could have happened in this simulation, because it shows the reality of an A/B test on low conversion rates. With p ≈ 0.04, the standard deviation of the number of conversions over 1000 trials is √(1000·0.04·0.96) ≈ 6.2: a true difference of 6 expected conversions is of the same order as the noise. A thousand observations per arm are not enough, and they were not even close.

The lesson about method is in the comparison of the two right-hand columns. The frequentist at 1000 observations reports p = 0.4272: «not significant», which in everyday language becomes «no difference» and is a false conclusion — absence of evidence is not evidence of absence. The Bayesian reports P(B > A) = 0.2148, which reads directly: I give B roughly a one in five chance of being better, declaring the direction and the weakness of the evidence at once. Note also the n = 1 row: the frequentist test is undefined, while the posterior returns 0.5004, that is, «I know nothing about this». Having a sensible answer when the data are zero is a property, not a flourish.

Trade-off, and the honesty it requires. Let no one conclude that the Bayesian method «finds» the difference: it does not find it, and on this sample it points the wrong way. No method can extract from 1000 observations information that is not there. The difference is that one of the two reports its uncertainty in an interpretable form and the other produces a binary threshold that invites a logical error. Anyone wanting to see the case in which B emerges must change the seed or increase N: that is an exercise in the mini-project, and the point of that exercise is to discover how many observations are really needed.

The grid posterior: the conjugate is a shortcut, not magic

A legitimate suspicion remains: that conjugates work because the prior was chosen to make the sums come out. Let us dismantle it by computing the same posterior in the most stupid and most general way possible — evaluating prior × likelihood on a dense grid and normalising by numerical integration — and comparing it with the closed form. We use arm B of the previous experiment, with all 1000 observations.

ka, kb = int(osservazioni_a.sum()), int(osservazioni_b.sum())
g = np.linspace(0, 1, 20001)
nonnorm = g ** kb * (1 - g) ** (N - kb)        # flat prior x likelihood
dens = nonnorm / np.trapezoid(nonnorm, g)      # normalising is the evidence, computed by hand
chiusa = stats.beta(1 + kb, 1 + N - kb)
print("maximum density difference: %.3e" % np.abs(dens - chiusa.pdf(g)).max())
print("grid mean         : %.6f" % float(np.trapezoid(g * dens, g)))
print("closed-form mean  : %.6f" % chiusa.mean())
print("mean WITHOUT normalising: %.3e  (a meaningless number)"
      % float(np.trapezoid(g * nonnorm, g)))

Output obtained:

maximum density difference: 3.809e-12
grid mean         : 0.037924
closed-form mean  : 0.037924
mean WITHOUT normalising: 1.023e-72  (a meaningless number)

What the code does. It builds the posterior by brute force on 20,001 points and compares it with Beta(1 + k, 1 + n − k). The last line deliberately computes the mean forgetting to divide by the integral.

What it tells us. The two densities agree to 3.8e-12 and the means to all six printed digits: the conjugate adds nothing, it merely avoids doing this computation. You can always fall back on the grid when the prior you actually want is not conjugate — the price is that it costs M^d points in d dimensions, and already at d = 5 with 100 points per axis that is ten billion evaluations. The last line is the commonest trap: without normalising you get 1.023e-72, a number that is not a mean, not a probability and not anything, but which a program prints without complaint.

When one number is not enough: the bimodal posterior

The advantage of the grid is that it shows the shape, and the shape is sometimes the information. Suppose a bimodal belief, perfectly sensible in practice: «either this channel converts at around 20%, or at around 80%, depending on which segment ended up in it». You write it as a mixture of two Betas. Then 5 successes out of 10 arrive, which is exactly the point in the middle and does not help to choose.

g = np.linspace(1e-6, 1 - 1e-6, 200001)
prior_mix = 0.5 * stats.beta(20, 80).pdf(g) + 0.5 * stats.beta(80, 20).pdf(g)
vero = g ** 5 * (1 - g) ** 5                   # likelihood: 5 successes out of 10
post_mix = prior_mix * vero
post_mix /= np.trapezoid(post_mix, g)
moda = g[int(post_mix.argmax())]
media = float(np.trapezoid(g * post_mix, g))
cdf = np.concatenate([[0.0], np.cumsum((post_mix[1:] + post_mix[:-1]) / 2 * np.diff(g))])
lo, hi = np.interp([0.025, 0.975], cdf, g)
print("mode  : %.4f" % moda)
print("mean  : %.4f" % media)
print("density at the mean: %.4f  density at the mode: %.4f"
      % (float(np.interp(media, g, post_mix)), float(post_mix.max())))
print("95%% credible interval: [%.4f, %.4f]" % (lo, hi))

Output obtained:

mode  : 0.7778
mean  : 0.5000
density at the mean: 0.0000  density at the mode: 5.0138
95% credible interval: [0.1648, 0.8352]

What the code does. It builds the posterior as the product of a mixture of two Betas and the likelihood of 5 successes out of 10, normalises it, and extracts mode, mean and credible interval (the last by inverting the distribution function computed from trapezoidal cumulative sums).

What it tells us. Here all three summaries lie. The mode is 0.7778 and suggests «the channel converts well»: it hides a second, almost equally tall peak around 0.22. The mean is 0.5000, and the next line is the most instructive in the output: the density at the mean is 0.0000 against 5.0138 at the mode — the mean falls in the valley between the two peaks, that is, at the value the posterior considers least plausible of all. And the interval [0.1648, 0.8352] is correct but covers almost everything, because it has to contain both humps and the valley between them.

The rule to take away. Before summarising a posterior with a number, look at it. If it is unimodal and symmetric, mean, median and mode coincide and the choice is irrelevant; in every other case the choice of summary is a substantive decision, and if the posterior is bimodal the honest answer is not a number but «there are two scenarios, here they are».

When Bayes is the wrong tool

Every lesson in this course has the section in which the tool just built is set aside. Here there are three cases, and they are common.

When you cannot defend the prior to a colleague. The test is operational: try saying out loud how many observations your prior is worth and where that number comes from. If the answer is «it seemed reasonable», you have put into the result information that does not exist, and we have measured it: coverage 0.1649 against 0.9694. Use a flat prior, or collect more data: an invented prior is worse than no prior, because it is invisible in the final result.

When the posterior is not computable and sampling does not pay for itself. Outside the conjugates you need iterative methods that cost orders of magnitude more than a point estimate, must be diagnosed, and introduce one more way of being wrong: a chain that has not converged produces plausible and false numbers. If the model has a million parameters, is retrained every night and is judged on prediction error, use the MAP — which is regularisation — and you have 90% of the benefit at 2% of the cost.

When you need a decision, not a distribution. «P(B > A) = 0.21» becomes a decision only when you add a loss function: how much it costs to switch variant by mistake and how much it costs not to switch. In regulated settings the threshold is often imposed from outside, and presenting a posterior where a pre-registered procedure is expected is a process problem, not a statistical one. The mature answer is to report both readings, as in the A/B test table.

Typical mistakes (symptom, cause, remedy)

1. Calling a credible interval a «confidence interval». Symptom: a slide says «we are 95% confident the rate is between 0.386 and 0.861», computed from a posterior. Remedy: «credibility» when the number comes from a posterior, «confidence» when it comes from a procedure with proven coverage, and the prior written next to the first.
2. A prior that dominates at small n and goes unnoticed at large n. Symptom: the result has been stable for months, then on a segment with little data it becomes absurd. Cause: the same prior weighs 2% on 1000 observations and 67% on 5. Remedy: always print the ratio of pseudo-counts to sample size.
3. Believing that «uninformative» is a property of the prior. Symptom: Beta(1, 1) is used and the work happens on the logit. Cause: flatness does not survive a change of variable. Remedy: ask «flat in which coordinates?»; if invariance is needed, Jeffreys.
4. Forgetting the normalisation in the grid approximation. Symptom: a «mean» of 1.023e-72, or probabilities that do not sum to 1. Cause: prior × likelihood is not a density. Remedy: divide by np.trapezoid(dens, griglia), with an assert on the integral.
5. Summarising a bimodal posterior with the mode alone. Symptom: the point estimate jumps abruptly when a little data is added. Cause: the maximum hops from one peak to the other. Remedy: plot the posterior; if it has two humps, report both.
6. Predicting with the MAP and declaring the model's uncertainty. Symptom: prediction intervals are too narrow and are violated more often than expected. Cause: the plug-in throws away the parameter variance (here 1.26 against 2.08). Remedy: the posterior predictive, or stating that the interval is conditional on θ = θ_MAP.

The six mistakes fall into three groups: two of language and of reading the prior (confidence in place of credibility, a prior that dominates at small n), two of computation (grid normalisation, summarising a bimodal with the mode), two of model (flatness that does not survive a change of variable, prediction with the MAP). When reviewing a Bayesian analysis it helps to go through the list in this order.

Best practice

• Write the prior in pseudo-counts. «Beta(2, 2)» says nothing to a meeting; «I behave as if I had already seen 2 successes and 2 failures» can be argued with by anyone.
• Always report sensitivity to the prior. Redo the calculation with at least two alternative priors, one flat and one opposed to yours. If the conclusions change, the result is the prior, not the data.
• Check the closed form with an independent method. A grid in one dimension costs milliseconds and here found agreement to 1.42e-14.
• Look at the posterior before summarising it. Mode, mean and median coincide only in the symmetric case, which is not the norm.
• Fix the seeds and declare the versions (default_rng(42), Python 3.14.4, NumPy 2.3.5, SciPy 1.16.3, scikit-learn 1.8.0): the simulations are reproducible only this way.
• Declare when the data are simulated. Every number in this lesson comes from synthetic data with a known true rate: that is what makes it possible to judge the conclusions, and it is what real data lack.

The six rules are meant to make every Bayesian analysis open to challenge and reproducible: the first two make the prior and its influence visible, the third and the fourth check the computation and the shape of the posterior, the last two make the numbers repeatable. The one that costs least and is worth most is prior sensitivity: two more computations.

Minimum glossary

• Prior (prior): the distribution of the parameter before seeing the data, p(θ).
• Likelihood (likelihood): p(D | θ) read as a function of θ; it is not a density in θ and does not integrate to 1.
• Posterior (posterior): the distribution of the parameter after the data, p(θ | D).
• Evidence (evidence, marginal likelihood): p(D) = ∫ p(D | θ) p(θ) dθ, the normalising constant; it is the integral that makes the general case hard.
• Conjugate (conjugate prior): a family of priors which, multiplied by a given likelihood, returns a posterior of the same family.
• Pseudo-count (pseudo-count): the parameter of a conjugate prior read as a number of fake observations already seen (a and b in the Beta).
• MAP (maximum a posteriori): the value of θ maximising the posterior; it coincides with the MLE when the prior is constant.
• Credible interval (credible interval): an interval containing a declared share of the posterior probability. A statement about the parameter.
• Confidence interval (confidence interval): a procedure whose frequency of covering the true value is declared. A statement about the procedure.
• Posterior predictive (posterior predictive): p(ỹ | D) = ∫ p(ỹ | θ) p(θ | D) dθ, the prediction that averages over all the parameter uncertainty.

The ten terms are ordered like the lesson: first the three pieces of the formula and the evidence, then conjugate, pseudo-count and MAP, finally the two intervals and the predictive. After reading prior, likelihood and posterior you can use the glossary as a reminder: the only asymmetry to remember is that the likelihood is not a density in the parameter.

Mini-project

Build a script that checks the three central claims of the lesson, each with an independent check that can fail.

1. Write aggiorna_beta(a, b, osservazioni) and verify that updating in one block and one observation at a time gives the same result. Print mean, mode and credible interval;
2. write map_gaussiano(X, y, sigma2, tau2) in closed form, (XᵀX + λI)⁻¹Xᵀy with λ = σ²/τ², and verify that it coincides with Ridge(alpha=λ, fit_intercept=False). Then let τ² tend to a huge number: it must converge to plain least squares;
3. write copertura_credibile(...), which measures by simulation the 95% coverage with n = 5 under three priors: flat, wrong and spot-on. Document the case in which it collapses, and answer: a flaw in the interval or the declared price of a belief?

The three requests check the three central claims of the lesson with three independent tests: block updating against one observation at a time, the Gaussian MAP against ridge, the simulated coverage under three priors. Each check must be able to fail, otherwise it verifies nothing.

Deliverable: the script prints the three blocks of results, and for the third a comment line saying which prior you would use in production and why.

Challenge

• Easy. With a Beta(2, 2) prior and 7 successes out of 10, write the posterior by hand and compute its mean and mode. Explain why they do not coincide and in which case they would.
• Medium. Derive the Normal-Normal conjugate by completing the square, without looking back at the lesson. Then check the limiting case τ² → ∞: it must return N(x̄, σ²/n), that is, the frequentist result.
• Hard. Prove that the MAP under a Laplace prior p(β_j) ∝ exp(−|β_j|/b) is the lasso, making the penalty coefficient explicit as a function of σ² and b. Verify it on a design with many irrelevant columns and count how many coefficients are exactly zero. Mind the normalisation of Lasso, which minimises (1/2n)‖y − Xβ‖² + α‖β‖₁: between its α and your MAP coefficient there is a factor of 2n.
• Expert. Leave the conjugates behind: a logistic likelihood with a Gaussian prior, a posterior with no closed form. Approximate it by sampling (a Metropolis algorithm with Gaussian steps is enough), then compare the mode of the samples with the numerical MAP. If they do not agree within sampling error the chain has not converged, and you check that by starting two chains from distant points.
• Enterprise. Take a metric from your organisation with little data: the adoption rate of a new feature, the frequency of a rare fault, the conversion of a newly opened channel. Write a prior you can defend, declaring it in pseudo-counts and saying where the number comes from. Report the posterior with its credible interval and attach the sensitivity table: the same conclusion under at least three priors, including the flat one and one opposed to yours. If the conclusion changes, the result is your belief, not the measurement.

The five levels climb from a computation by hand to an analysis with your own data: the first two consolidate conjugates and the posterior mean, the third links the lasso to the Laplace prior, the fourth leaves the conjugates with a sampler, the last brings the prior into pseudo-counts on a real metric of your organisation. Start from the level where you get the first computation wrong.

The code exercises of this lesson are in a single file, esercizi.py (about 10 KB). It holds the skeleton of the three mini-project exercises, with the parts to complete and, at the end, the commented solutions with the expected values. It requires Python 3.12+, NumPy 2.x, SciPy and scikit-learn. Run it with python3 esercizi.py: it works even with the parts to complete left empty and tells you which ones are missing.

Recap

Mental map. p(θ | D) ∝ p(D | θ) p(θ): the hard part is not the product, it is the evidence p(D) = ∫ p(D|θ)p(θ)dθ. The conjugate is the family you recognise inside the product, and the pseudo-counts say what the prior is worth (table above). MAP = argmax [log p(D|θ) + log p(θ)]: Gaussian prior ⇒ ridge with λ = σ²/τ², Laplace prior ⇒ lasso. Credible = a probability about the parameter; confidence = a property of the procedure: identical at n = 40 (0.9478 both), divergent at n = 5 with a wrong prior (0.1649 against 0.9694). Predictive = integrate over the posterior, do not freeze the MAP.

Cheat sheet.

3 concepts: the posterior is a distribution, not an estimate; the conjugate avoids the integral by recognising the shape; regularisation is a prior written in another notation. 3 mistakes: calling credibility «confidence»; using a prior you cannot defend on little data; predicting with the MAP while declaring the model's uncertainty. 3 best practices: declare the prior in pseudo-counts; publish the sensitivity to the prior; check the closed form with a grid. 1 exercise: the mini-project of this lesson.

Next lesson: sampling, that is, how to estimate what you cannot compute.

Questions to take further

• Outside the conjugates, how do you obtain samples from a distribution whose unnormalised shape alone is known? (Lesson 29, module 02: sampling and Monte Carlo methods.)
• What guarantees that a chain of dependent samples represents the right distribution, and how do you notice that it has not got there yet? (Lesson 30, module 02: Markov chains and the stationary distribution.)
• If the MAP under a Gaussian prior is ridge, which prior corresponds to dropout or to early stopping? (An open question: some correspondences are proven, others are analogies; the criterion is whether anyone has written down the density.)

The three questions are hooks to the lessons that follow: the first two introduce sampling and Markov chains, the third is an open question about which regularisers have a proven prior and which only an analogy. They help you decide what to read next.

The same series also includes: 27 - Entropy, Cross-Entropy and KL Divergence: Why Loss Functions Are All the Same Thing (the previous lesson), 29 - Sampling and Monte Carlo Methods: Estimating What You Cannot Compute (the next lesson).