28 - Statistica bayesiana: prior, posterior, coniugati e perché la regolarizzazione è una credenza

Corso Machine Learning Engineer · Modulo 02 · Matematica e statistica · Lezione 13

Obiettivo. Applicare il teorema di Bayes non a un evento ma a un parametro, cioè trattare ciò che non conosci come una distribuzione che i dati modificano. Alla fine sai aggiornare una credenza con una formula sola, derivare per intero i due coniugati che si incontrano sempre (Beta-Binomiale e Normale-Normale), riconoscere che la stima MAP con un prior gaussiano è la regressione ridge e con un prior di Laplace è il lasso, e distinguere un intervallo di credibilità da un intervallo di confidenza senza barare sulle parole.

Falso mito da sfatare. «La statistica bayesiana è soggettiva, quella frequentista è obiettiva.» No: entrambe partono da assunzioni dichiarate. La differenza è che la bayesiana ti obbliga a scrivere la tua credenza iniziale in una formula, dove si vede e si discute; la frequentista la nasconde nella scelta del modello, della famiglia di test e del livello di significatività. Una credenza scritta è più obiettiva di una credenza implicita.

Introduzione

Nella lezione 18, dedicata alla probabilità, abbiamo usato il teorema di Bayes su eventi: il test diagnostico che risulta positivo, la porta che il conduttore apre. Lì l'incognita era «quale evento si è verificato», e la risposta era un numero. Nella lezione 20 sulla statistica abbiamo lavorato in modo interamente frequentista: stimatori, errore standard, intervalli di confidenza, con il parametro trattato come una costante sconosciuta ma fissa. Questa lezione fa la terza mossa, quella che collega le due: applica Bayes ai parametri. L'incognita non è più un evento, è θ, e la risposta non è un numero ma una distribuzione.

Non è un esercizio di stile. Due cose già scritte in questo corso cambiano aspetto quando le guardi da qui. La prima è la correzione +1 della lezione 32 sul Naive Bayes: non è un trucco contro i conteggi a zero, è un prior dichiarato. La seconda è la penalizzazione L2 della lezione 44 su sovradattamento e regolarizzazione, presentata lì come un termine che «scoraggia i coefficienti grandi»: è esattamente il logaritmo di un prior gaussiano, e il suo peso λ ha una formula chiusa che dice quanto credi al prior rispetto ai dati.

Risultato finale di questa lezione. Saprai prendere un problema con pochi dati, scrivere un prior che sai difendere, derivare la posterior a mano quando il coniugato lo permette e per approssimazione quando non lo permette, e dire con precisione che cosa hai il diritto di concludere.

Roadmap della lezione.

1. Dalla verosimiglianza alla posterior: i quattro pezzi della formula e quale dei quattro è il problema.
2. I coniugati Beta-Binomiale e Normale-Normale, derivati per intero, e la correzione di Laplace che ne risulta.
3. MAP contro MLE: la ridge e il lasso smascherati come prior.
4. Intervallo di credibilità contro intervallo di confidenza, misurati per simulazione, e la posterior predittiva.
5. Un A/B test aggiornato giorno per giorno, contro lo stesso dato letto da un frequentista.
6. Quando Bayes è lo strumento sbagliato, con i criteri per accorgersene prima.

La scaletta segue il filo che tiene insieme la lezione: la formula con i suoi quattro pezzi, due coniugati derivati a mano, la scoperta che ridge e lasso sono prior scritti in un'altra notazione, poi le misure per simulazione e infine i limiti. Chi ha già visto Bayes sugli eventi può saltare alla terza voce, ma il passaggio dai coniugati alla regolarizzazione dipende dai pseudo-conteggi della seconda.

Dalla verosimiglianza alla posterior

L'analogia: il medico che conosce già il reparto

Due medici vedono lo stesso paziente, con lo stesso sintomo, e misurano lo stesso valore in laboratorio. Il primo lavora in un reparto di medicina tropicale, dove la malattia rara in questione si presenta ogni settimana. Il secondo lavora in un ambulatorio di periferia, dove non si è mai vista. Il primo considera la diagnosi probabile, il secondo la considera improbabile e ordina un secondo esame. Hanno ragione entrambi, e non perché la medicina sia un'opinione: perché partono da una prevalenza diversa, e la prevalenza è un'informazione, non un pregiudizio.

In termini di inferenza: il dato di laboratorio è la verosimiglianza, identica per i due; la prevalenza del reparto è il prior, diverso; la conclusione è la posterior, diversa di conseguenza. Se sostituisci «reparto» con «dominio applicativo» e «prevalenza» con «quello che sai del parametro prima di guardare questi dati», hai la statistica bayesiana. E hai anche il suo prezzo: due persone con prior diversi, sugli stessi dati, concludono cose diverse, e l'unico modo di discuterne è che il prior sia scritto.

La formula, pezzo per pezzo

Sia θ il parametro (una probabilità, una media, un vettore di coefficienti) e D i dati osservati. Il teorema di Bayes, applicato a θ, dice:

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

e, poiché p(D) non dipende da θ, nella pratica si usa quasi sempre nella forma proporzionale:

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

Leggila come una frase: la credenza dopo è la credenza prima, corretta da quanto bene ogni valore del parametro spiega i dati. I quattro pezzi — prior p(θ), verosimiglianza p(D | θ), posterior p(θ | D) ed evidenza p(D) — sono definiti nel glossario in fondo. Qui serve fissare l'asimmetria che fa inciampare quasi tutti: prior e posterior sono densità su θ e integrano a 1, la verosimiglianza no. Come funzione di θ, p(D | θ) non è una distribuzione di probabilità e non c'è motivo per cui la sua area debba valere qualcosa.

È l'evidenza a rendere difficile il caso generale. Quel simbolo di integrale è innocuo con un parametro e devastante con mille: se θ ha 500 componenti, ∫ p(D | θ) p(θ) dθ è un integrale in 500 dimensioni che nessuna formula chiude e nessuna griglia copre. Tutta la macchineria che incontrerai dopo — griglie, Monte Carlo, catene di Markov, inferenza variazionale — esiste per aggirare quell'integrale. I coniugati sono i pochi casi fortunati in cui si risolve a mano.

Il coniugato Beta-Binomiale, derivato per intero

Il caso più piccolo che contenga tutto: una moneta truccata, un tasso di conversione, una probabilità di guasto. Il parametro θ ∈ [0, 1] è una probabilità di successo, i dati sono k successi su n prove indipendenti. Serve un prior su θ, cioè una densità sull'intervallo [0, 1]: la famiglia Beta fa esattamente questo.

Beta(θ; a, b) = θ^{a−1} (1 − θ)^{b−1} / B(a, b), con B(a, b) la costante di normalizzazione.

La verosimiglianza binomiale, come funzione di θ, è

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

dove C(n, k) è il coefficiente binomiale, che non dipende da θ e quindi finirà nella costante. Moltiplichiamo prior e verosimiglianza, tenendo solo i fattori che dipendono da θ:

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

Le potenze con la stessa base si sommano negli esponenti:

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

Ora la mossa decisiva, e vale la pena dire esattamente che cos'è: riconoscere la forma. Quell'espressione è θ elevato a «qualcosa meno uno» per (1 − θ) elevato a «qualcos'altro meno uno»: è il nucleo di una Beta. Siccome esiste una sola densità su [0, 1] proporzionale a quel nucleo, la costante di normalizzazione è obbligata e non serve calcolarla. Quindi

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

Questa è la definizione operativa di coniugato: una famiglia di prior tale che, moltiplicata per quella verosimiglianza, restituisce un prior della stessa famiglia. L'integrale dell'evidenza non è stato calcolato: è stato evitato, perché sapevamo in anticipo quale densità avremmo ottenuto.

La lettura che conta: a e b sono pseudo-conteggi

Guarda la formula di aggiornamento: a riceve i successi, b gli insuccessi. Allora a e b si leggono come dati finti che hai già visto: un prior Beta(2, 2) dichiara «prima di cominciare, mi comporto come se avessi già osservato 2 successi e 2 insuccessi». Questo rende il prior discutibile in termini concreti: non «quanto sono confidente», ma «a quante osservazioni equivale la mia convinzione». È la domanda giusta da fare in riunione, e si risponde con un numero.

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)   media = %.4f" % (a, b, a / (a + b)))
print("dati       : %d successi su %d   frequenza = %.4f" % (k, n, k / n))
print("posterior  : Beta(%.0f, %.0f)   media = %.4f" % (a_post, b_post, post.mean()))
print("moda posteriore (a-1)/(a+b-2) =", round((a_post - 1) / (a_post + b_post - 2), 4))
print("stima di massima verosimiglianza k/n =", round(k / n, 4))
print("intervallo di credibilita al 95%% : [%.4f, %.4f]" % post.interval(0.95))
print("pseudo-conteggi del prior: %.0f successi e %.0f insuccessi gia visti" % (a, b))

Output ottenuto (Python 3.14.4, SciPy 1.16.3):

prior      : Beta(2, 2)   media = 0.5000
dati       : 7 successi su 10   frequenza = 0.7000
posterior  : Beta(9, 5)   media = 0.6429
moda posteriore (a-1)/(a+b-2) = 0.6667
stima di massima verosimiglianza k/n = 0.7
intervallo di credibilita al 95% : [0.3857, 0.8614]
pseudo-conteggi del prior: 2 successi e 2 insuccessi gia visti

Cosa fa il codice. Applica la formula appena derivata — Beta(2, 2) più 7 successi e 3 insuccessi dà Beta(9, 5) — e stampa tre riassunti della stessa posterior: media 9/14, moda 8/12 e intervallo di credibilità al 95% a code uguali.

Che cosa ci dice. La media posteriore 0,6429 sta fra la media del prior (0,5) e la frequenza osservata (0,7), più vicina alla seconda perché i dati valgono 10 osservazioni contro le 4 del prior. Media e moda non coincidono (0,6429 contro 0,6667), perché la Beta è asimmetrica quando a ≠ b: scegliere l'una o l'altra è una decisione, non un dettaglio. E l'intervallo [0,386 ; 0,861] è larghissimo, il messaggio più utile dell'output: con dieci prove non sai quasi niente, e la posterior ha la decenza di dirlo.

Da dove viene la correzione +1 del Naive Bayes

Nella lezione 32 il Naive Bayes stimava la probabilità di ciascuna parola in ciascuna classe come conteggio / totale, e si rompeva su una parola mai vista: un fattore zero azzerava il prodotto e con esso la classe. Il rimedio era aggiungere 1 a tutti i conteggi. Vista da qui, quella correzione ha un nome e una derivazione: il caso con K esiti invece di 2 ha come coniugato la Dirichlet, che sta alla multinomiale come la Beta sta alla binomiale, e la sua posterior è Dirichlet(α₁ + c₁, …, α_K + c_K). Con α = (1, …, 1) — il prior uniforme sul simplesso — la media posteriore della categoria j è

(c_j + 1) / (N + K)

che è, carattere per carattere, la correzione di Laplace. Non è un tappo contro una divisione per zero: è la media di una posterior sotto il prior «ogni esito è possibile e ne ho visto uno finto di ciascuno».

conteggi = np.array([12, 5, 0, 3])          # conteggi di 4 parole in una classe
K = len(conteggi)
N = conteggi.sum()
mle = conteggi / N
media_post = (conteggi + 1) / (N + K)       # media di Dirichlet(1+c1, ..., 1+cK)
print("conteggi osservati        :", conteggi)
print("stima MLE (c/N)           :", np.round(mle, 4))
print("media posteriore (c+1)/(N+K):", np.round(media_post, 4))
print("somma delle medie         :", round(float(media_post.sum()), 10))
print("la parola mai vista passa da", mle[2], "a", round(float(media_post[2]), 4))

Output ottenuto:

conteggi osservati        : [12  5  0  3]
stima MLE (c/N)           : [0.6  0.25 0.   0.15]
media posteriore (c+1)/(N+K): [0.5417 0.25   0.0417 0.1667]
somma delle medie         : 1.0
la parola mai vista passa da 0.0 a 0.0417

Cosa fa il codice. Prende quattro conteggi, uno dei quali a zero, e mette a confronto la stima di massima verosimiglianza c/N con la media posteriore (c+1)/(N+K). La riga della somma verifica che le medie posteriori siano ancora una distribuzione: 1.0 esatto.

Che cosa ci dice. La parola mai vista passa da 0,0 a 0,0417: non è più impossibile, è solo improbabile. Il prezzo è visibile nella prima colonna, che scende da 0,60 a 0,5417: la probabilità regalata all'esito assente viene sottratta agli altri, perché il totale resta 1. Chi usa la correzione senza saperlo sta facendo questo scambio al buio.

La regola da portare via. Ogni volta che «aggiungi un piccolo numero per stabilità» stai scegliendo un prior. L'unica differenza fra te e un bayesiano è se lo hai scritto.

Come si stringe la posterior

Con a = b = 2 gli pseudo-conteggi sono 4: pesano il 29% del totale a n = 10, il 3,8% a n = 100, lo 0,4% a n = 1000. La deviazione standard della posterior, intanto, scende come 1/√n.

Il coniugato Normale-Normale, derivato per intero

Il secondo coniugato da conoscere riguarda una media. Osserviamo n valori x₁, …, x_n da una normale di media μ sconosciuta e varianza σ² nota, e mettiamo su μ un prior normale N(μ₀, τ²). L'ipotesi di varianza nota è una semplificazione vera e dichiarata: serve a far uscire la formula in una pagina. Il caso con varianza anch'essa sconosciuta esiste, ha come coniugato la Normale-Gamma-inversa, e si trova in Gelman cap. 3.

Scriviamo il logaritmo del prodotto prior × verosimiglianza, buttando via tutto ciò che non dipende da μ:

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

Espandiamo i due quadrati e raccogliamo per potenze di μ. Dal primo termine: Σᵢ (xᵢ − μ)² = Σᵢ xᵢ² − 2μ Σᵢ xᵢ + n μ². Dal secondo: (μ − μ₀)² = μ² − 2μμ₀ + μ₀². Mettendo insieme e scartando le costanti:

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

dove x̄ è la media campionaria. Introduciamo due nomi che rendono tutto leggibile: la precisione è l'inverso della varianza. Chiamiamo λ₀ = 1/τ² la precisione del prior e λ_D = n/σ² la precisione dei dati, e poniamo λ_n = λ₀ + λ_D. Allora

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

Ora il completamento del quadrato, svolto: per rendere l'espressione in parentesi un quadrato perfetto (μ − m)² serve m = (λ_D x̄ + λ₀ μ₀)/λ_n, perché (μ − m)² = μ² − 2μm + m² e il termine m² è una costante che possiamo aggiungere e sottrarre liberamente:

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

Esponenziando, p(μ | D) ∝ exp[ −(μ − m)²/(2/λ_n) ], che è il nucleo di una normale di media m e varianza 1/λ_n. Quindi

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

Leggila come due frasi, perché sono le due cose da ricordare di tutta la sezione. La media posteriore è la media pesata di media del prior e media dei dati, con pesi uguali alle rispettive precisioni. E la precisione posteriore è la somma delle precisioni: l'informazione si somma, e si somma nella coordinata giusta, che è l'inverso della varianza, non la varianza.

I tre coniugati incontrati finora si leggono tutti con la stessa chiave, e la colonna di destra è quella che si porta in riunione.

Come leggere l'ultima colonna. Il peso dei dati cresce linearmente in n, quello del prior resta fermo: è per questo che il prior svanisce. E il numero della terza riga, σ²/τ², non è una coincidenza di notazione: è esattamente il λ della ridge che deriveremo fra poco. Il coefficiente di regolarizzazione è il numero di osservazioni che vale la tua credenza.

mu0, tau2 = 100.0, 25.0                     # prior: N(100, 25)
sigma2 = 16.0                               # varianza nota dei dati
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("media campionaria        :", round(float(dati.mean()), 4))
print("precisione del prior 1/t2:", round(prec_prior, 6))
print("precisione dei dati n/s2 :", round(prec_dati, 6))
print("media posteriore         :", round(float(mu_post), 4))
print("varianza posteriore      :", round(float(var_post), 6))

# verifica indipendente: la stessa posterior calcolata per forza bruta su una griglia
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())      # sottraggo il massimo: evita l'underflow
dens /= np.trapezoid(dens, griglia)
media_griglia = np.trapezoid(griglia * dens, griglia)
var_griglia = np.trapezoid((griglia - media_griglia) ** 2 * dens, griglia)
print("media su griglia         :", round(float(media_griglia), 4))
print("varianza su griglia      :", round(float(var_griglia), 6))
print("differenza massima       : %.2e" % max(abs(media_griglia - mu_post), abs(var_griglia - var_post)))

Output ottenuto (seed 42, NumPy 2.3.5):

media campionaria        : 106.7571
precisione del prior 1/t2: 0.04
precisione dei dati n/s2 : 0.5
media posteriore         : 106.2565
varianza posteriore      : 1.851852
media su griglia         : 106.2565
varianza su griglia      : 1.851852
differenza massima       : 1.42e-14

Cosa fa il codice. La prima metà applica la formula chiusa, la seconda la verifica per forza bruta su una griglia di 200.001 punti. Il dettaglio da rubare è dens = np.exp(logpost - logpost.max()): sommando otto log-densità si arriva a numeri come −10³, e exp(−1000) è zero in virgola mobile; sottrarre il massimo prima di esponenziare non cambia il risultato, perché è un fattore costante che la normalizzazione elimina.

Che cosa ci dice. I due percorsi coincidono a 1,42e-14, cioè all'ultima cifra della doppia precisione: la derivazione è giusta. Nel merito: la precisione dei dati è 0,5 contro 0,04 del prior, quindi i dati pesano dodici volte e mezza, e infatti la media posteriore 106,2565 è molto più vicina a 106,7571 (i dati) che a 100 (il prior). La varianza posteriore 1,85 è più piccola sia di τ² = 25 sia di σ²/n = 2: combinare due informazioni dà qualcosa di più preciso di entrambe, e questo è il contenuto operativo della somma delle precisioni.

MAP contro MLE, e la regolarizzazione smascherata

Una posterior è una distribuzione, ma a volte serve un numero: per scrivere una predizione in un sistema, per confrontarsi con una stima frequentista, per inizializzare un'ottimizzazione. Il riassunto più usato è il punto di massima densità posteriore, la stima MAP (maximum a posteriori). Poiché il logaritmo è crescente e l'evidenza non dipende da θ:

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

Confrontala con la stima di massima verosimiglianza, θ_MLE = argmax_θ log p(D | θ). La differenza è un solo addendo: il logaritmo del prior. E qui accade la cosa interessante di tutta la lezione, perché quell'addendo, per i prior più comuni, è un oggetto che hai già visto con un altro nome.

Prior gaussiano: il termine che compare è la ridge

Prendi una regressione lineare con rumore gaussiano di varianza σ² nota: y = Xβ + ε, ε ~ N(0, σ²). La log-verosimiglianza, a meno di costanti, è −‖y − Xβ‖²/(2σ²). Metti su ogni coefficiente un prior gaussiano indipendente e centrato in zero, β_j ~ N(0, τ²): il suo logaritmo, a meno di costanti, è −‖β‖²/(2τ²). Sommando:

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

Massimizzare equivale a minimizzare l'opposto; moltiplicando per 2σ², che è una costante positiva e non cambia il punto di minimo:

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

Questa è, carattere per carattere, la regressione ridge della lezione 44, con

λ = σ²/τ²

Leggi la formula di λ, perché dice tutto quello che la lezione 44 non poteva dire. Il numeratore è il rumore dei dati, il denominatore è l'ampiezza della tua credenza sui coefficienti. Penalizzi molto quando i dati sono rumorosi (σ² grande) oppure quando credi che i coefficienti siano piccoli (τ² piccolo). E il caso λ = 0 non è «nessuna regolarizzazione»: è τ² → ∞, cioè il prior che dichiara ogni valore dei coefficienti equiprobabile, compreso 10¹². Non è l'assenza di una credenza, è una credenza particolarmente strana.

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):
    """Meno il logaritmo della posterior, a meno di costanti."""
    res = y - X @ beta
    return res @ res / (2 * sigma2_r) + beta @ beta / (2 * tau2_r)


def grad(beta):
    """Gradiente analitico: serve per arrivare a una tolleranza severa."""
    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("differenza massima   : %.2e" % np.abs(sol.x - ridge.coef_).max())
print("MLE (lambda = 0)     :", np.round(np.linalg.lstsq(X, y, rcond=None)[0], 6))

Output ottenuto (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]
differenza massima   : 1.33e-15
MLE (lambda = 0)     : [ 2.147018 -1.448232  0.118071  0.794285 -0.202626]

Cosa fa il codice. Costruisce una regressione sintetica 60×5 con coefficienti veri noti (dati simulati da noi, dichiarati tali), calcola il MAP minimizzando −log p(β | D) con gradiente analitico, e lo confronta con Ridge(alpha=σ²/τ²). Due dettagli non decorativi: fit_intercept=False, perché il prior è sui coefficienti e non su un'intercetta, e gtol=1e-12, perché con la tolleranza di default la differenza resterebbe intorno a 1e-6, abbastanza da far dubitare di un'identità che invece è esatta.

Che cosa ci dice. I due vettori coincidono a 1,33e-15: non sono «simili», sono la stessa cosa calcolata in due modi. L'ultima riga aggiunge il confronto utile: il MLE è più grande in modulo su tutte le componenti (2,147 contro 2,015). È la contrazione verso zero imposta dal prior, e sulla terza componente — coefficiente vero 0 — porta la stima da 0,118 a 0,105: nella direzione giusta, ma senza arrivare a zero. Il perché sta nella sottosezione seguente.

Prior di Laplace: il termine che compare è il lasso

Cambiamo prior lasciando tutto il resto identico. Al posto della gaussiana mettiamo su ogni coefficiente una distribuzione di Laplace centrata in zero, p(β_j) ∝ exp(−|β_j| / b): stessa idea («i coefficienti stanno vicino a zero»), ma con un picco appuntito nell'origine e code più pesanti. Il suo logaritmo, a meno di costanti, è −‖β‖₁ / b, dove ‖β‖₁ = Σ_j |β_j|. Ripetendo i passaggi di prima:

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

che è il lasso, con coefficiente di penalizzazione 2σ²/b. Prima di verificarlo c'è una trappola di conversione da disinnescare, perché fa sembrare sbagliata un'identità esatta. sklearn.linear_model.Lasso non minimizza quella quantità: minimizza (1/2n)·‖y − Xβ‖² + α·‖β‖₁, cioè la stessa cosa divisa per 2n. Uguagliando i due rapporti fra penalizzazione e termine dei dati si ottiene 2nα = 2σ²/b, cioè

α = σ² / (n · b)

Chi salta questo passaggio confronta due modelli regolarizzati in modo diverso, vede coefficienti diversi e conclude che la corrispondenza non vale. Vale: è la conversione a essere sbagliata.

Per vedere la differenza che conta serve un disegno in cui ci sia qualcosa da azzerare: 80 righe, 12 colonne di cui 8 con coefficiente vero esattamente nullo.

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                      # prior gaussiano  N(0, tau2)      -> ridge
scala_b = 0.05                  # prior di Laplace 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("coefficienti veri :", np.round(beta_vero, 3))
print("ridge             :", np.round(ridge.coef_, 3))
print("lasso             :", np.round(lasso.coef_, 3))
print("zeri esatti ridge :", int((ridge.coef_ == 0).sum()), "su", d)
print("zeri esatti lasso :", int((lasso.coef_ == 0).sum()), "su", d)


# Il lasso è davvero il MAP? Valuto l'obiettivo in 20000 perturbazioni casuali:
# se è il minimo, nessuna deve fare meglio.
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: coefficiente piu piccolo in modulo = %.6f" % np.abs(ridge.coef_).min())
print("obiettivo MAP nel lasso          : %.6f" % obiettivo_map(lasso.coef_))
print("obiettivo MAP nella ridge        : %.6f" % obiettivo_map(ridge.coef_))
print("miglior perturbazione su 20000   : %.6f" % val_perturbati.min())
print("perturbazioni che battono il lasso:",
      int((val_perturbati < obiettivo_map(lasso.coef_) - 1e-9).sum()))

# Condizione di ottimalità: |X_j^T r| / n <= alpha per ogni coefficiente spento.
residuo = y - X @ lasso.coef_
correlazioni = np.abs(X.T @ residuo) / n_righe
spenti = lasso.coef_ == 0
print("coefficienti spenti: max |X_j^T r|/n = %.6f  (soglia alpha = %.6f)"
      % (correlazioni[spenti].max(), alpha_lasso))
print("coefficienti accesi: min |X_j^T r|/n = %.6f  (deve valere ~alpha)"
      % correlazioni[~spenti].min())

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

coefficienti veri : [ 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.   ]
zeri esatti ridge : 0 su 12
zeri esatti lasso : 7 su 12
ridge: coefficiente piu piccolo in modulo = 0.005879
obiettivo MAP nel lasso          : 148.465847
obiettivo MAP nella ridge        : 162.941647
miglior perturbazione su 20000   : 148.876569
perturbazioni che battono il lasso: 0
coefficienti spenti: max |X_j^T r|/n = 0.165989  (soglia alpha = 0.250000)
coefficienti accesi: min |X_j^T r|/n = 0.249988  (deve valere ~alpha)

Cosa fa il codice. Adatta ridge e lasso allo stesso disegno con le due conversioni appena derivate, conta gli zeri esatti (confronto == 0, non «minore di una soglia»), e verifica in due modi indipendenti che il lasso sia il MAP: valuta l'obiettivo in 20.000 vettori ottenuti perturbando la soluzione, e controlla la condizione di ottimalità, per cui un coefficiente può stare a zero solo se la sua correlazione col residuo |Xⱼᵀr|/n non supera α.

Che cosa ci dice. Il lasso azzera esattamente 7 coefficienti su 12, la ridge nessuno: il suo più piccolo è 0,005879, e lo resterà per qualunque λ finito. La verifica MAP è netta: obiettivo 148,47 contro 162,94 della ridge, e zero perturbazioni su 20.000 lo migliorano. L'ultima coppia di righe è la spiegazione geometrica in due numeri: i coefficienti spenti hanno correlazione col residuo al massimo 0,166, sotto la soglia 0,25; quelli accesi ce l'hanno a 0,249988, appoggiata alla soglia. Il lasso spegne tutto ciò che non supera α, e tiene il resto appena sopra.

Il lasso non è gratis. Delle 8 colonne davvero inutili ne azzera 7: la nona resta a 0,008, accesa per caso. E contrae i coefficienti veri più della ridge — 0,582 contro 0,829 dove la verità è 0,9: quella distorsione sui sopravvissuti è il prezzo della selezione, e in previsione può costare più di quanto la sparsità faccia guadagnare. Perché uno azzera e l'altro no: la derivata di β² in zero è zero, quindi vicino all'origine la penalizzazione L2 smette di spingere; la derivata di |β| in zero non esiste — c'è uno spigolo — e la spinta vale ±1 fino all'ultimo, finché non è il dato a opporsi con forza almeno α. Lo spigolo del prior di Laplace è la selezione delle variabili.

Intervallo di credibilità contro intervallo di confidenza

Le due definizioni, affiancate, con la differenza che non è una sfumatura di linguaggio.

Intervallo di credibilità al 95%. Un intervallo [l, u] tale che ∫_l^u p(θ | D) dθ = 0,95: data questa posterior, il 95% della mia credenza sul parametro sta lì dentro. È un'affermazione su θ, e ha senso parlare della probabilità che θ ci stia.

Intervallo di confidenza al 95%. Una procedura che, applicata a dati generati ripetutamente dal modello, produce intervalli contenenti il valore vero nel 95% dei casi: non so nulla su questo intervallo in particolare; so che il metodo sbaglia una volta su venti. È un'affermazione sulla procedura, e θ è una costante: non è lecito dire «c'è il 95% di probabilità che θ stia qui dentro».

La frase che le separa, da mandare a memoria: il credibile è una probabilità sul parametro, il confidenziale è una proprietà della procedura. Chiamare «intervallo di confidenza» un intervallo di credibilità non è pedanteria terminologica: promette una garanzia di copertura che il secondo non offre, come stiamo per misurare.

Primo fatto: su dati abbondanti sono indistinguibili

Con un prior uniforme e un campione di dimensione ragionevole, le due procedure danno quasi gli stessi numeri. Misuriamolo: 10.000 esperimenti binomiali con n = 40 e p = 0,30, e per ciascuno costruiamo l'intervallo credibile al 95% da Beta(1 + k, 1 + n − k) e l'intervallo di confidenza di Wilson, che è il modo corretto di fare un intervallo per una proporzione (l'approssimazione normale ingenua p̂ ± 1,96·√(p̂(1−p̂)/n) ha copertura notoriamente scadente sulle proporzioni estreme).

def wilson(k, n, z=1.959963984540054):
    """Intervallo di confidenza di Wilson per una proporzione."""
    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 esperimenti in un colpo
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("esperimenti simulati          :", prove)
print("copertura intervallo credibile: %.4f" % cop_cred)
print("copertura intervallo Wilson   : %.4f" % cop_conf)
print("ampiezza media credibile      : %.4f" % float((cred[1] - cred[0]).mean()))
print("ampiezza media Wilson         : %.4f" % float((conf[1] - conf[0]).mean()))

Output ottenuto (seed 11, 10.000 esperimenti):

esperimenti simulati          : 10000
copertura intervallo credibile: 0.9478
copertura intervallo Wilson   : 0.9478
ampiezza media credibile      : 0.2706
ampiezza media Wilson         : 0.2697

Cosa fa il codice. Simula i 10.000 esperimenti con una sola chiamata vettorizzata, costruisce i due intervalli in parallelo — stats.beta accetta array di parametri — e calcola la frazione che contiene 0,30.

Che cosa ci dice. Le due coperture sono identiche a quattro decimali, 0,9478, e le ampiezze medie differiscono di 0,0009. Il risultato va detto per quello che è: con n = 40 e un prior uniforme, scegliere l'una o l'altra procedura non cambia nulla di misurabile. Chi presenta la statistica bayesiana come intrinsecamente più prudente o più potente qui non ha argomenti. Da notare anche che 0,9478 non è 0,95: su una variabile discreta la copertura esatta non è raggiungibile per nessuno dei due metodi, e oscilla intorno al valore nominale al variare di n e p.

Secondo fatto: su pochi dati il prior si vede tutto

La differenza vera compare dove l'informazione dei dati è scarsa. Ripetiamo con n = 5, e aggiungiamo un prior informativo e sbagliato: Beta(8, 2), che ha media 0,8 mentre la verità è 0,30. Il confronto si fa a tre: prior sbagliato, prior piatto, 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)    # prior sbagliato
cred_pia = stats.beta(1 + ks, 1 + n_pic - ks).interval(0.95)    # prior piatto
conf = wilson(ks, n_pic)
for nome, iv in (("credibile con Beta(8,2)", cred_inf),
                 ("credibile con 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: copertura %.4f | ampiezza media %.4f" % (nome, cop, amp))

Output ottenuto (seed 12, 10.000 esperimenti, n = 5):

credibile con Beta(8,2): copertura 0.1649 | ampiezza media 0.4584
credibile con Beta(1,1): copertura 0.9694 | ampiezza media 0.6011
Wilson                 : copertura 0.9694 | ampiezza media 0.5905

Cosa fa il codice. Identico al precedente, con n ridotto a 5 e tre intervalli invece di due. L'unica differenza sostanziale è la coppia di parametri del prior nella prima riga.

Che cosa ci dice. Qui si vede la frattura. Il prior piatto e Wilson hanno la stessa copertura, 0,9694, mentre il prior informativo scende a 0,1649: l'intervallo credibile costruito con Beta(8, 2) contiene il valore vero in meno di un esperimento su sei. E attenzione alla seconda colonna: è anche il più stretto dei tre, 0,4584 contro 0,59-0,60. Un intervallo più stretto e molto più sbagliato è la combinazione peggiore possibile, perché si presenta come più informativo.

Va detto con precisione di che cosa questo sia la prova. Non è un difetto del metodo: l'intervallo credibile riporta fedelmente il 95% di una credenza che include 10 pseudo-osservazioni centrate su 0,8, e con n = 5 quei pseudo-conteggi valgono il doppio dei dati. La copertura frequentista al 95% è una garanzia che il metodo bayesiano offre solo mediando sul prior: se p fosse davvero estratto da Beta(8, 2), tornerebbe nominale. Il crollo a 0,1649 misura quanto costa un prior che non sai difendere.

La regola da portare via. Con tanti dati la scelta della scuola statistica è una questione di gusto e di strumenti. Con pochi dati è la scelta più importante che fai, e il prior va negoziato prima di guardare i risultati, non dopo.

Posterior predittiva

Fin qui abbiamo parlato del parametro. Ma quasi sempre la domanda vera non è «quanto vale θ», è «che cosa succederà la prossima volta». La risposta corretta non si ottiene fissando θ al suo valore migliore e usando quello: si ottiene mediando le predizioni su tutta la posterior. Si chiama posterior predittiva:

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

In parole: per ogni valore plausibile del parametro chiedi che cosa prevede, e pesa la sua previsione con quanto quel valore è plausibile. Il confronto è con la predizione plug-in, p(ỹ | θ_MAP), che è ciò che fa chiunque prenda una stima puntuale e la infili nel modello. Le due coincidono solo se la posterior è concentrata in un punto; altrimenti la plug-in sottostima l'incertezza, perché butta via la varianza del parametro e tiene solo quella del rumore.

Nel caso Beta-Binomiale l'integrale si risolve e dà una distribuzione con un nome: la Beta-Binomiale, disponibile in SciPy come stats.betabinom. Mettiamola alla prova dove la differenza si vede, cioè con pochissimi dati: 4 successi su 5, e la domanda «su 10 prove future, quanti successi?».

n_vis, k_vis, m = 5, 4, 10
ap, bp = 1 + k_vis, 1 + n_vis - k_vis          # posterior con prior piatto
pred = stats.betabinom(m, ap, bp)              # posterior predittiva, forma chiusa
p_map = (ap - 1) / (ap + bp - 2)               # moda della posterior
plug = stats.binom(m, p_map)                   # predizione plug-in col MAP
print("posterior: Beta(%d, %d)   MAP = %.4f" % (ap, bp, p_map))
print("predittiva Beta-Binomiale: media %.4f  dev.std %.4f  intervallo 95%% [%d, %d]"
      % (pred.mean(), pred.std(), pred.ppf(0.025), pred.ppf(0.975)))
print("plug-in Binomiale(MAP)   : media %.4f  dev.std %.4f  intervallo 95%% [%d, %d]"
      % (plug.mean(), plug.std(), plug.ppf(0.025), plug.ppf(0.975)))

Output ottenuto:

posterior: Beta(5, 2)   MAP = 0.8000
predittiva Beta-Binomiale: media 7.1429  dev.std 2.0825  intervallo 95% [3, 10]
plug-in Binomiale(MAP)   : media 8.0000  dev.std 1.2649  intervallo 95% [5, 10]

Cosa fa il codice. Costruisce la posterior Beta(5, 2) da 4 successi su 5 con prior piatto, poi due distribuzioni predittive per 10 prove future: quella corretta, che integra su tutta la posterior, e quella plug-in, che congela θ alla moda 0,8. Di entrambe stampa media, deviazione standard e intervallo al 95%.

Che cosa ci dice. La deviazione standard predittiva corretta è 2,08, quella plug-in 1,26: la seconda dichiara un'incertezza inferiore del 39% rispetto a quella che ha diritto di dichiarare. E l'intervallo al 95% passa da [3, 10] a [5, 10]: la plug-in esclude gli esiti 3 e 4, perfettamente possibili se θ è 0,55 anziché 0,8, e con cinque osservazioni 0,55 è tutt'altro che escluso. Anche la media si sposta, da 8,0 a 7,14, perché la posterior Beta(5, 2) ha media 5/7 ≈ 0,714, inferiore alla moda.

Trade-off. La predittiva corretta costa un integrale e restituisce intervalli più larghi, che in una presentazione sembrano un risultato peggiore: è il motivo per cui quasi nessuno la usa, ed è un errore di incentivi — l'intervallo stretto non è più preciso, è solo più sicuro di sé. Con n grande la differenza svanisce e la plug-in diventa legittima; la soglia si controlla confrontando le due deviazioni standard, come qui sopra.

Il prior «non informativo» non esiste

A questo punto la tentazione è ovvia: se il prior è il problema, usiamone uno che non dica nulla. Il candidato naturale è il prior uniforme, Beta(1, 1) sulle probabilità, «piatto» sull'intervallo. L'obiezione, che risale a Jeffreys, è che piatto non è una proprietà del parametro ma delle coordinate in cui lo scrivi. Se sei indifferente fra tutti i valori di θ, non puoi essere contemporaneamente indifferente fra tutti i valori di log θ o del logit φ = log(θ/(1−θ)): la regola di cambio di variabile per le densità, p_φ(φ) = p_θ(θ) · |dθ/dφ|, lo vieta.

Il prior che è invariante esiste e si chiama prior di Jeffreys: è proporzionale alla radice dell'informazione di Fisher, p(θ) ∝ √I(θ), e per la binomiale dà Beta(1/2, 1/2). Verifichiamo numericamente entrambe le affermazioni.

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)    # moltiplico per |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)         # trasformato in phi
jeff_diretta = np.sqrt(theta * (1 - theta))            # Jeffreys calcolato in phi
rapporto = jeff_su_phi / jeff_diretta
print("prior piatto su theta, densita indotta su phi nei punti phi = -2, 0, 2:")
for punto in (-2.0, 0.0, 2.0):
    i = int(np.abs(phi - punto).argmin())
    print("  phi = %5.2f -> densita %.6f" % (phi[i], indotta_su_phi[i]))
print("rapporto Jeffreys trasformato / Jeffreys diretto: min %.10f  max %.10f"
      % (rapporto.min(), rapporto.max()))

Output ottenuto:

prior piatto su theta, densita indotta su phi nei punti phi = -2, 0, 2:
  phi = -2.00 -> densita 0.104992
  phi = -0.00 -> densita 0.250000
  phi =  2.00 -> densita 0.104992
rapporto Jeffreys trasformato / Jeffreys diretto: min 1.0000000000  max 1.0000000000

Cosa fa il codice. Trasforma il prior uniforme nelle coordinate del logit moltiplicandolo per |dθ/dφ| = θ(1−θ) e ne valuta la densità in tre punti. Poi fa la stessa trasformazione su Jeffreys e la confronta con Jeffreys calcolato direttamente in φ, dove l'informazione di Fisher vale θ(1−θ).

Che cosa ci dice. La densità indotta dal prior uniforme sul logit vale 0,25 al centro e 0,105 a φ = ±2: non è piatta, è una campana che privilegia i valori centrali di θ. Chi dichiara «uso un prior non informativo sulla probabilità» sta, nelle coordinate del logit — quelle in cui lavora qualunque regressione logistica — dichiarando una credenza informativa e centrata. Jeffreys invece passa il test: rapporto 1,0000000000 sia al minimo sia al massimo su 200.001 punti, cioè esattamente invariante.

Jeffreys però non è una soluzione universale: in più dimensioni può essere impropria e portare a posterior improprie. La regola da portare via: non chiedere «questo prior è informativo?», chiedi «informativo rispetto a quale parametrizzazione?».

L'esperimento: un A/B test aggiornato giorno per giorno

Mettiamo tutto al lavoro su un caso realistico. Due varianti di una pagina, tassi di conversione veri 0,040 e 0,046 — cioè B è davvero migliore, del 15% relativo, che è un miglioramento considerevole per una conversione. I dati sono simulati da noi con default_rng(7), quindi sappiamo la verità e possiamo giudicare le conclusioni. Osserviamo il flusso a quattro tappe — dopo 1, 10, 100 e 1000 osservazioni per braccio — e a ogni tappa calcoliamo due cose sullo stesso dato: la quantità bayesiana P(θ_B > θ_A | D), stimata campionando dalle due posterior, e il p-value del test frequentista di confronto fra due proporzioni, che è lo strumento della lezione 6 sugli esperimenti controllati.

rng = np.random.default_rng(7)
p_a, p_b = 0.040, 0.046                    # tassi VERI: li conosciamo perché simuliamo noi
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) | p-value frequentista")
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 di A
    sb = campionatore.beta(1 + kb, 1 + t - kb, size=200000)      # posterior di B
    prob = float((sb > sa).mean())                               # P(B > A) per campionamento
    pooled = (ka + kb) / (2 * t)
    se = np.sqrt(pooled * (1 - pooled) * 2 / t)
    if se == 0:
        pval = float("nan")                                      # nessuna conversione: test non definito
    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,
                                             "non definito" if np.isnan(pval) else "%.4f" % pval))

Output ottenuto (seed 7 per i dati, seed 99 per il campionamento):

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

Cosa fa il codice. Genera i due flussi una volta sola e ne legge prefissi via via più lunghi, come un esperimento che accumula traffico. P(B > A) non ha una forma chiusa comoda: si stima campionando 200.000 valori da ciascuna posterior e contando quante volte B supera A. È il primo uso del campionamento come stimatore e non come verifica, ed è il filo che la prossima lezione raccoglie.

Che cosa ci dice, e non è quello che ci aspettavamo. B ha il tasso vero più alto, eppure P(B > A) scende: da 0,50 a 0,21. A 1000 osservazioni per braccio abbiamo osservato 44 conversioni per A e 37 per B, cioè 4,4% contro 3,7%: la variante peggiore sta vincendo. Non è un errore del codice ed è la cosa più utile che potesse succedere in questa simulazione, perché mostra la realtà di un A/B test su conversioni basse. Con p ≈ 0,04, la deviazione standard del numero di conversioni su 1000 prove è √(1000·0,04·0,96) ≈ 6,2: una differenza vera di 6 conversioni attese è dello stesso ordine del rumore. Mille osservazioni per braccio non bastano, e non ci andavano nemmeno vicino.

La lezione sul metodo sta nel confronto fra le due colonne di destra. Il frequentista a 1000 osservazioni riporta p = 0,4272: «non significativo», che nel linguaggio corrente diventa «nessuna differenza» ed è una conclusione falsa — l'assenza di prova non è prova di assenza. Il bayesiano riporta P(B > A) = 0,2148, che si legge direttamente: do a B circa una probabilità su cinque di essere migliore, e dichiara insieme la direzione e la debolezza dell'evidenza. Si noti anche la riga n = 1: il test frequentista è non definito, mentre la posterior restituisce 0,5004, cioè «non ne so nulla». Avere una risposta sensata quando i dati sono zero è una proprietà, non un vezzo.

Trade-off, e l'onestà che serve. Non si concluda che il metodo bayesiano «trova» la differenza: non la trova, e in questo campione indica la direzione sbagliata. Nessun metodo può estrarre da 1000 osservazioni un'informazione che non c'è. La differenza è che uno dei due riporta la propria incertezza in una forma interpretabile e l'altro produce una soglia binaria che invita a un errore logico. Chi volesse vedere il caso in cui B emerge deve cambiare seme o aumentare N: è un esercizio del mini-progetto, e il punto di quell'esercizio è scoprire quante osservazioni servano davvero.

La posterior a griglia: il coniugato è una scorciatoia, non una magia

Resta un sospetto legittimo: che i coniugati funzionino perché si è scelto il prior che fa tornare i conti. Smontiamolo calcolando la stessa posterior nel modo più stupido e più generale possibile — valutando prior × verosimiglianza su una griglia fitta e normalizzando per integrazione numerica — e confrontandola con la forma chiusa. Usiamo il braccio B dell'esperimento precedente, con tutte e 1000 le osservazioni.

ka, kb = int(osservazioni_a.sum()), int(osservazioni_b.sum())
g = np.linspace(0, 1, 20001)
nonnorm = g ** kb * (1 - g) ** (N - kb)        # prior piatto x verosimiglianza
dens = nonnorm / np.trapezoid(nonnorm, g)      # la normalizzazione è l'evidenza, calcolata a mano
chiusa = stats.beta(1 + kb, 1 + N - kb)
print("massima differenza fra densita: %.3e" % np.abs(dens - chiusa.pdf(g)).max())
print("media su griglia  : %.6f" % float(np.trapezoid(g * dens, g)))
print("media forma chiusa: %.6f" % chiusa.mean())
print("media SENZA normalizzare: %.3e  (numero privo di senso)"
      % float(np.trapezoid(g * nonnorm, g)))

Output ottenuto:

massima differenza fra densita: 3.809e-12
media su griglia  : 0.037924
media forma chiusa: 0.037924
media SENZA normalizzare: 1.023e-72  (numero privo di senso)

Cosa fa il codice. Costruisce la posterior per forza bruta su 20.001 punti e la confronta con Beta(1 + k, 1 + n − k). L'ultima riga calcola apposta la media dimenticando di dividere per l'integrale.

Che cosa ci dice. Le due densità coincidono a 3,8e-12 e le medie a tutte e sei le cifre stampate: il coniugato non aggiunge niente, evita solo di fare questo calcolo. Si può sempre ricadere sulla griglia quando il prior che vuoi davvero non è coniugato — il prezzo è che costa M^d punti in d dimensioni, e già a d = 5 con 100 punti per asse sono dieci miliardi di valutazioni. L'ultima riga è la trappola più comune: senza normalizzare si ottiene 1,023e-72, un numero che non è una media, non è una probabilità e non è nulla, ma che un programma stampa senza protestare.

Quando un solo numero non basta: la posterior bimodale

Il vantaggio della griglia è che mostra la forma, e la forma a volte è l'informazione. Supponiamo una credenza bimodale, perfettamente sensata in pratica: «o questo canale converte intorno al 20%, o intorno all'80%, a seconda di quale segmento ci è finito dentro». La si scrive come miscela di due Beta. Poi arrivano 5 successi su 10, che è il punto esattamente in mezzo e non aiuta a scegliere.

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                   # verosimiglianza: 5 successi su 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("moda  : %.4f" % moda)
print("media : %.4f" % media)
print("densita nella media: %.4f  densita nella moda: %.4f"
      % (float(np.interp(media, g, post_mix)), float(post_mix.max())))
print("intervallo di credibilita al 95%%: [%.4f, %.4f]" % (lo, hi))

Output ottenuto:

moda  : 0.7778
media : 0.5000
densita nella media: 0.0000  densita nella moda: 5.0138
intervallo di credibilita al 95%: [0.1648, 0.8352]

Cosa fa il codice. Costruisce la posterior come prodotto di una miscela di due Beta e della verosimiglianza di 5 successi su 10, la normalizza, e ne estrae moda, media e intervallo di credibilità (quest'ultimo invertendo la funzione di ripartizione calcolata per somme cumulate trapezoidali).

Che cosa ci dice. Qui tutti e tre i riassunti mentono. La moda è 0,7778 e suggerisce «il canale converte molto»: nasconde un secondo picco quasi altrettanto alto intorno a 0,22. La media è 0,5000, e la riga successiva è la più istruttiva dell'output: la densità nel punto della media è 0,0000 contro 5,0138 nella moda — la media cade nella valle fra i due picchi, cioè nel valore che la posterior considera meno plausibile di tutti. E l'intervallo [0,1648 ; 0,8352] è corretto ma copre quasi tutto, perché deve contenere entrambe le gobbe e la valle in mezzo.

La regola da portare via. Prima di riassumere una posterior con un numero, guardala. Se è unimodale e simmetrica, media, mediana e moda coincidono e la scelta è irrilevante; in ogni altro caso la scelta del riassunto è una decisione sostanziale, e se la posterior è bimodale la risposta onesta non è un numero ma «ci sono due scenari, eccoli».

Quando Bayes è lo strumento sbagliato

Ogni lezione di questo corso ha la sezione in cui lo strumento appena costruito viene messo da parte. Qui i casi sono tre, e sono frequenti.

Quando non sai difendere il prior davanti a un collega. Il test è operativo: prova a dire ad alta voce a quante osservazioni equivale il tuo prior e da dove viene quel numero. Se la risposta è «mi sembrava ragionevole», hai messo nel risultato un'informazione che non esiste, e l'abbiamo misurata: copertura 0,1649 contro 0,9694. Usa un prior piatto, o raccogli più dati: un prior inventato è peggio di nessun prior, perché è invisibile nel risultato finale.

Quando la posterior non è calcolabile e il campionamento non si ripaga. Fuori dai coniugati servono metodi iterativi che costano ordini di grandezza più di una stima puntuale, vanno diagnosticati e introducono un modo di sbagliare in più: una catena non arrivata a convergenza produce numeri plausibili e falsi. Se il modello ha un milione di parametri, viene riaddestrato ogni notte ed è giudicato sull'errore di previsione, usa il MAP — che è la regolarizzazione — e hai il 90% del beneficio al 2% del costo.

Quando serve una decisione, non una distribuzione. «P(B > A) = 0,21» diventa una decisione solo quando ci aggiungi una funzione di perdita: quanto costa cambiare variante per sbaglio e quanto costa non cambiarla. In contesti regolati la soglia è spesso imposta dall'esterno, e presentare una posterior dove è attesa una procedura pre-registrata è un problema di processo, non di statistica. La risposta matura è riportare entrambe le letture, come nella tabella dell'A/B test.

Errori tipici (sintomo, causa, rimedio)

1. Chiamare «intervallo di confidenza» un intervallo di credibilità. Sintomo: una slide dice «siamo confidenti al 95% che il tasso sia fra 0,386 e 0,861», calcolata da una posterior. Rimedio: «credibilità» quando il numero viene da una posterior, «confidenza» quando viene da una procedura con copertura dimostrata, e il prior scritto accanto al primo.
2. Un prior che domina su n piccolo e passa inosservato su n grande. Sintomo: il risultato è stabile da mesi, poi su un segmento con pochi dati diventa assurdo. Causa: lo stesso prior pesa il 2% su 1000 osservazioni e il 67% su 5. Rimedio: stampa sempre il rapporto fra pseudo-conteggi e dimensione del campione.
3. Credere che «non informativo» sia una proprietà del prior. Sintomo: si usa Beta(1, 1) e si lavora sul logit. Causa: la piattezza non sopravvive al cambio di variabile. Rimedio: chiedersi «piatto in quali coordinate?»; se serve invarianza, Jeffreys.
4. Dimenticare la normalizzazione nell'approssimazione a griglia. Sintomo: una «media» di 1,023e-72, o probabilità che non sommano a 1. Causa: prior × verosimiglianza non è una densità. Rimedio: dividere per np.trapezoid(dens, griglia), con un assert sull'integrale.
5. Riassumere con la sola moda una posterior bimodale. Sintomo: la stima puntuale cambia bruscamente aggiungendo pochi dati. Causa: il massimo salta da un picco all'altro. Rimedio: disegnare la posterior; se ha due gobbe, riportarle entrambe.
6. Predire con il MAP e dichiarare l'incertezza del modello. Sintomo: gli intervalli di previsione sono troppo stretti e vengono violati più spesso del previsto. Causa: la plug-in butta via la varianza del parametro (qui 1,26 contro 2,08). Rimedio: posterior predittiva, o dichiarare che l'intervallo è condizionato a θ = θ_MAP.

I sei errori si dividono in tre gruppi: due di linguaggio e di lettura del prior (confidenza al posto di credibilità, prior che domina a n piccolo), due di calcolo (normalizzazione della griglia, riassumere una bimodale con la moda), due di modello (piattezza che non sopravvive al cambio di variabile, previsione col MAP). Nella revisione di un'analisi bayesiana conviene scorrere la lista in quest'ordine.

Best practice

• Scrivi il prior in pseudo-conteggi. «Beta(2, 2)» non dice nulla a una riunione; «mi comporto come se avessi già visto 2 successi e 2 insuccessi» è discutibile da chiunque.
• Riporta sempre la sensibilità al prior. Rifai il conto con almeno due prior alternativi, uno piatto e uno opposto al tuo. Se le conclusioni cambiano, il risultato è il prior, non i dati.
• Verifica la forma chiusa con un metodo indipendente. Una griglia in una dimensione costa millisecondi e qui ha trovato un accordo a 1,42e-14.
• Guarda la posterior prima di riassumerla. Moda, media e mediana coincidono solo nel caso simmetrico, che non è la norma.
• Fissa i seed e dichiara le versioni (default_rng(42), Python 3.14.4, NumPy 2.3.5, SciPy 1.16.3, scikit-learn 1.8.0): le simulazioni sono riproducibili solo così.
• Dichiara quando i dati sono simulati. Tutti i numeri di questa lezione vengono da dati sintetici con un tasso vero noto: è ciò che permette di giudicare le conclusioni, ed è ciò che manca nei dati reali.

Le sei regole sono pensate per rendere discutibile e riproducibile ogni analisi bayesiana: le prime due rendono visibile il prior e la sua influenza, la terza e la quarta controllano il calcolo e la forma della posterior, le ultime due rendono ripetibili i numeri. Quella che costa meno e vale di più è la sensibilità al prior: due conti in più.

Glossario minimo

• Prior (prior): la distribuzione del parametro prima di vedere i dati, p(θ).
• Verosimiglianza (likelihood): p(D | θ) letta come funzione di θ; non è una densità in θ e non integra a 1.
• Posterior (posterior): la distribuzione del parametro dopo i dati, p(θ | D).
• Evidenza (evidence, marginal likelihood): p(D) = ∫ p(D | θ) p(θ) dθ, la costante di normalizzazione; è l'integrale che rende difficile il caso generale.
• Coniugato (conjugate prior): famiglia di prior che, moltiplicata per una data verosimiglianza, restituisce una posterior della stessa famiglia.
• Pseudo-conteggio (pseudo-count): il parametro di un prior coniugato letto come numero di osservazioni finte già viste (a e b nella Beta).
• MAP (maximum a posteriori): il valore di θ che massimizza la posterior; coincide col MLE quando il prior è costante.
• Intervallo di credibilità (credible interval): intervallo che contiene una quota dichiarata della probabilità posteriore. Affermazione sul parametro.
• Intervallo di confidenza (confidence interval): procedura la cui frequenza di copertura del valore vero è dichiarata. Affermazione sulla procedura.
• Posterior predittiva (posterior predictive): p(ỹ | D) = ∫ p(ỹ | θ) p(θ | D) dθ, la previsione che media su tutta l'incertezza del parametro.

I dieci termini sono ordinati come la lezione: prima i tre pezzi della formula e l'evidenza, poi coniugato, pseudo-conteggio e MAP, infine i due intervalli e la predittiva. Dopo aver letto prior, verosimiglianza e posterior puoi usare il glossario come promemoria: l'unica asimmetria da ricordare è che la verosimiglianza non è una densità nel parametro.

Mini-progetto

Costruisci uno script che verifichi le tre affermazioni centrali della lezione, ciascuna con un controllo indipendente che possa fallire.

1. Scrivi aggiorna_beta(a, b, osservazioni) e verifica che aggiornare in blocco e un'osservazione alla volta dia lo stesso risultato. Stampa media, moda e intervallo di credibilità;
2. scrivi map_gaussiano(X, y, sigma2, tau2) in forma chiusa, (XᵀX + λI)⁻¹Xᵀy con λ = σ²/τ², e verifica che coincida con Ridge(alpha=λ, fit_intercept=False). Poi fai tendere τ² a un numero enorme: deve convergere ai minimi quadrati puri;
3. scrivi copertura_credibile(...) che misura per simulazione la copertura al 95% con n = 5 sotto tre prior: piatto, sbagliato e azzeccato. Documenta il caso in cui crolla, e rispondi: difetto dell'intervallo o prezzo dichiarato di una credenza?

Le tre richieste verificano le tre affermazioni centrali della lezione con tre controlli indipendenti: l'aggiornamento in blocco contro quello a una osservazione per volta, il MAP gaussiano contro la ridge, la copertura simulata sotto tre prior. Ciascun controllo deve poter fallire, altrimenti non verifica nulla.

Consegna: lo script stampa i tre blocchi di risultati, e per il terzo una riga di commento che dica quale prior useresti in produzione e perché.

Challenge

• Easy. Con prior Beta(2, 2) e 7 successi su 10, scrivi la posterior a mano e calcolane media e moda. Spiega perché non coincidono e in quale caso coinciderebbero.
• Medium. Deriva il coniugato Normale-Normale completando il quadrato, senza riguardare la lezione. Verifica poi il caso limite τ² → ∞: deve restituire N(x̄, σ²/n), cioè il risultato frequentista.
• Hard. Dimostra che il MAP con prior di Laplace p(β_j) ∝ exp(−|β_j|/b) è il lasso, esplicitando il coefficiente di penalizzazione in funzione di σ² e b. Verificalo su un disegno con molte colonne irrilevanti e conta quanti coefficienti sono esattamente zero. Attenzione alla normalizzazione di Lasso, che minimizza (1/2n)‖y − Xβ‖² + α‖β‖₁: fra il suo α e il tuo coefficiente MAP c'è un fattore 2n.
• Expert. Esci dai coniugati: verosimiglianza logistica e prior gaussiano, posterior senza forma chiusa. Approssimala campionando (un Metropolis a passi gaussiani basta), poi confronta la moda dei campioni con il MAP numerico. Se non coincidono entro l'errore di campionamento la catena non è arrivata a convergenza, e lo verifichi facendo partire due catene da punti lontani.
• Enterprise. Prendi una metrica della tua organizzazione con pochi dati: il tasso di adozione di una funzionalità nuova, la frequenza di un guasto raro, la conversione di un canale appena aperto. Scrivi un prior che sai difendere, dichiarandolo in pseudo-conteggi e dicendo da dove viene il numero. Riporta la posterior con il suo intervallo di credibilità e allega la tabella di sensibilità: la stessa conclusione sotto almeno tre prior, incluso quello piatto e uno opposto al tuo. Se la conclusione cambia, il risultato è la tua credenza, non la misura.

I cinque livelli salgono da un conto a mano a un'analisi con i tuoi dati: i primi due consolidano coniugati e media posteriore, il terzo collega il lasso al prior di Laplace, il quarto esce dai coniugati con un campionatore, l'ultimo porta il prior in pseudo-conteggi su una metrica reale della tua organizzazione. Parti dal livello in cui sbagli il primo conto.

Gli esercizi di codice di questa lezione sono in un solo file, esercizi.py (circa 10 KB). Lo scheletro dei tre esercizi del mini-progetto, con le parti da completare e, in fondo, le soluzioni commentate con i valori attesi. Richiede Python 3.12+, NumPy 2.x, SciPy e scikit-learn. Eseguilo con python3 esercizi.py: gira anche con le parti da completare lasciate vuote e ti dice quali mancano.

Recap

Mappa mentale. p(θ | D) ∝ p(D | θ) p(θ): la difficoltà non è il prodotto, è l'evidenza p(D) = ∫ p(D|θ)p(θ)dθ. Il coniugato è la famiglia che si riconosce nel prodotto, e gli pseudo-conteggi dicono quanto vale il prior (tabella sopra). MAP = argmax [log p(D|θ) + log p(θ)]: prior gaussiano ⇒ ridge con λ = σ²/τ², prior di Laplace ⇒ lasso. Credibile = probabilità sul parametro; confidenza = proprietà della procedura: identiche a n = 40 (0,9478 entrambe), divergenti a n = 5 con prior sbagliato (0,1649 contro 0,9694). Predittiva = integra sulla posterior, non congela il MAP.

Cheat sheet.

3 concetti: la posterior è una distribuzione, non una stima; il coniugato evita l'integrale riconoscendo la forma; la regolarizzazione è un prior scritto in un'altra notazione. 3 errori: chiamare «confidenza» la credibilità; usare un prior che non si sa difendere su pochi dati; predire col MAP dichiarando l'incertezza del modello. 3 best practice: dichiarare il prior in pseudo-conteggi; pubblicare la sensibilità al prior; verificare la forma chiusa con una griglia. 1 esercizio: il mini-progetto di questa lezione.

Prossima lezione: il campionamento, cioè come si stima quello che non si riesce a calcolare.

Domande per andare avanti

• Fuori dai coniugati, come si ottengono campioni da una distribuzione di cui si conosce solo la forma non normalizzata? (Lezione 29, modulo 02: campionamento e metodi Monte Carlo.)
• Che cosa garantisce che una catena di campioni dipendenti rappresenti la distribuzione giusta, e come ci si accorge che non ci è ancora arrivata? (Lezione 30, modulo 02: catene di Markov e distribuzione stazionaria.)
• Se il MAP con prior gaussiano è la ridge, qual è il prior che corrisponde al dropout o alla fermata anticipata? (Domanda aperta: alcune corrispondenze sono dimostrate, altre sono analogie; il criterio è se qualcuno ha scritto la densità.)

Le tre domande sono ganci verso le lezioni che seguono: le prime due introducono campionamento e catene di Markov, la terza è una domanda aperta su quali regolarizzatori abbiano un prior dimostrato e quali solo un'analogia. Servono a decidere che cosa leggere dopo.

Nella stessa serie trovi anche: 27 - Entropia, entropia incrociata e divergenza KL: perché le funzioni di costo sono tutte la stessa cosa (la lezione precedente), 29 - Campionamento e metodi Monte Carlo: stimare quello che non si riesce a calcolare (la lezione successiva).