11 - Stima Bayesiana: prior, posterior, MAP

Hai mostrato un nuovo pulsante a 20 visitatori e 8 hanno cliccato. La frazione osservata è 0,40, ma tu vuoi sapere altro: quanto puoi fidarti di quel numero, e che cosa cambia se in passato i pulsanti simili convertivano al 10%? La statistica frequentista risponde con uno stimatore e un intervallo di confidenza; quella bayesiana risponde con una distribuzione di probabilità sul valore vero, costruita combinando ciò che sapevi prima con ciò che hai misurato. In questo articolo, l'undicesimo della serie Probabilità e Statistica, calcoli quella distribuzione per intero su un caso in cui si scrive a mano, la confronti con la stima di massima verosimiglianza (MLE, il valore che rende più probabile il dato, visto nell'articolo sugli stimatori) e poi la fai calcolare a un campionatore quando la forma chiusa non c'è. Ogni risultato è verificato con un calcolo esatto e con una simulazione a seme fisso.

Dal teorema di Bayes al posterior

Chiama θ (theta) il parametro che non conosci, qui la probabilità vera che un visitatore clicchi, e D i dati. Il prior p(θ) descrive quanto è plausibile ciascun valore di θ prima di guardare i dati. La verosimiglianza p(D | θ) dice quanto sarebbero probabili quei dati se θ fosse proprio quel valore. Il teorema di Bayes li compone nel posterior: p(θ | D) = p(D | θ) · p(θ) / p(D). Il denominatore p(D), detto evidenza, è la somma (un integrale) del numeratore su tutti i θ possibili e serve solo a far sì che l'area valga 1: per ottenere la forma del posterior basta quindi moltiplicare verosimiglianza e prior.

Con k successi su n prove indipendenti la verosimiglianza è θ^k (1 − θ)^n−k. Se scegli come prior una distribuzione Beta(a, b), cioè una densità proporzionale a θ^a−1 (1 − θ)^b−1 sull'intervallo da 0 a 1, il prodotto ha la stessa forma con gli esponenti sommati: il posterior è Beta(a + k, b + n − k). Un prior con questa proprietà, per cui il posterior resta nella stessa famiglia, si dice coniugato. Leggi a e b come successi e fallimenti «virtuali» osservati prima del test: Beta(1, 1) è la densità uniforme, equivale a non aver visto nulla; Beta(10, 10) equivale a venti prove passate con metà successi.

Un posterior calcolato per intero: 8 successi su 20

Il programma qui sotto parte da Beta(1, 1), arriva a Beta(9, 13) e calcola le stime puntuali e gli intervalli. Prima però fa una verifica che non dipende dalla formula: costruisce il posterior a forza bruta, moltiplicando verosimiglianza e prior su una griglia di mille valori di θ e normalizzando. Se la formula coniugata fosse sbagliata, le due curve non coinciderebbero; invece differiscono per meno di 10^-14. L'HPD (highest posterior density) è l'intervallo più corto che contiene il 95% della probabilità; la funzione hpd lo trova provando tutti gli intervalli con quella massa e tenendo il più stretto. L'intervallo equal-tail (ETI) taglia invece il 2,5% da ciascun lato. Il Wilson è l'intervallo frequentista per una proporzione, calcolato con la formula di Wilson (1927).

import numpy as np
from scipy import stats

n, k = 20, 8                      # 8 successi su 20 prove
a0, b0 = 1, 1                     # prior uniforme Beta(1, 1)
a1, b1 = a0 + k, b0 + (n - k)     # posterior coniugato: Beta(9, 13)
post = stats.beta(a1, b1)

# Controllo: posterior calcolato a forza bruta su una griglia (verosimiglianza x prior, poi normalizzo)
theta = np.linspace(0.0005, 0.9995, 1000)
dx = theta[1] - theta[0]
grezzo = theta**k * (1 - theta)**(n - k) * stats.beta.pdf(theta, a0, b0)
griglia = grezzo / (grezzo.sum() * dx)
print(f"max |griglia - Beta(9,13)| = {np.abs(griglia - post.pdf(theta)).max():.1e}")

# Stime puntuali ed intervalli: forma chiusa
mle = k / n
mean = post.mean()
median = post.ppf(0.5)
mode = (a1 - 1) / (a1 + b1 - 2)           # MAP con prior uniforme = MLE
eti = (post.ppf(0.025), post.ppf(0.975))

def hpd(dist, massa=0.95, passi=20000):
    # tra tutti gli intervalli con probabilita' 'massa' prendo il piu' corto
    p_basso = np.linspace(0, 1 - massa, passi)
    lo = dist.ppf(p_basso)
    hi = dist.ppf(p_basso + massa)
    j = np.argmin(hi - lo)
    return lo[j], hi[j]

h = hpd(post)
z = stats.norm.ppf(0.975)
centro = (k + z**2 / 2) / (n + z**2)
semi = z * np.sqrt(k * (n - k) / n + z**2 / 4) / (n + z**2)
print(f"MLE {mle:.4f}  MAP {mode:.4f}  mediana {median:.4f}  media {mean:.4f}  sd {post.std():.4f}")
print(f"ETI 95%  [{eti[0]:.4f}, {eti[1]:.4f}]  ampiezza {eti[1]-eti[0]:.4f}")
print(f"HPD 95%  [{h[0]:.4f}, {h[1]:.4f}]  ampiezza {h[1]-h[0]:.4f}")
print(f"Wilson 95% [{centro-semi:.4f}, {centro+semi:.4f}]")

# Simulazione a seme fisso: deve confermare il calcolo esatto entro l'errore di Monte Carlo
rng = np.random.default_rng(2026)
N = 200_000
x = rng.beta(a1, b1, N)
print(f"simulazione ({N} estrazioni): media {x.mean():.4f} +- {x.std()/np.sqrt(N):.4f}, "
      f"q2.5% {np.quantile(x, .025):.4f}, q97.5% {np.quantile(x, .975):.4f}")

# Predittiva: probabilita' di almeno 5 successi in 10 prove nuove
pp = stats.betabinom(10, a1, b1)
plug = stats.binom(10, mle)
y = rng.binomial(10, x)                    # uno scenario per ogni estrazione del posterior
print(f"P(Y>=5) predittiva {1 - pp.cdf(4):.4f}  simulata {np.mean(y >= 5):.4f}  plug-in MLE {1 - plug.cdf(4):.4f}")
print(f"varianza di Y: predittiva {pp.var():.3f}  plug-in {plug.var():.3f}")

# Output:
# max |griglia - Beta(9,13)| = 9.1e-15
# MLE 0.4000  MAP 0.4000  mediana 0.4063  media 0.4091  sd 0.1025
# ETI 95%  [0.2182, 0.6156]  ampiezza 0.3974
# HPD 95%  [0.2129, 0.6097]  ampiezza 0.3968
# Wilson 95% [0.2188, 0.6134]
# simulazione (200000 estrazioni): media 0.4093 +- 0.0002, q2.5% 0.2180, q97.5% 0.6156
# P(Y>=5) predittiva 0.4033  simulata 0.4043  plug-in MLE 0.3669
# varianza di Y: predittiva 3.363  plug-in 2.400

Il MAP (maximum a posteriori, il punto più alto del posterior) vale 0,4000 e coincide con l'MLE: con prior uniforme il posterior è proporzionale alla verosimiglianza, e il massimo non si sposta. Media (0,4091) e mediana (0,4063) sono più alte perché la Beta(9, 13) ha la coda destra più lunga. L'intervallo equal-tail [0,2182; 0,6156] e il Wilson [0,2188; 0,6134] sono quasi identici: con 20 prove e un prior piatto i due paradigmi dicono quasi la stessa cosa. La simulazione con 200.000 estrazioni conferma il calcolo esatto: media 0,4093 contro 0,4091 (errore standard 0,0002) e quantili uguali alla terza cifra. L'ultima parte del programma calcola la predittiva: la probabilità di ottenere almeno 5 click in altri 10 visitatori. Vale 0,4033 integrando sull'incertezza di θ, contro 0,3669 che si ottiene «fidandosi» di 0,40 come se fosse il valore vero. La varianza del numero di click è 3,36 contro 2,40: il plug-in sottostima la dispersione perché ignora che θ non è noto.

Quale stima puntuale scegliere

Il posterior è la risposta completa, ma spesso serve un solo numero. Quale, dipende da che cosa ti costa sbagliare, cioè dalla funzione di perdita:

• perdita quadratica (l'errore pesa come il suo quadrato): la stima migliore è la media del posterior;
• perdita assoluta (l'errore pesa in proporzione): la mediana;
• perdita 0-1 (conta solo indovinare esattamente): il MAP, la moda.

Non esiste quindi una stima «bayesiana» unica: sono tre risposte a tre domande diverse, e coincidono solo quando il posterior è simmetrico e unimodale. Per la Beta(9, 13) differiscono di pochi millesimi, ma per un posterior molto asimmetrico, o con due picchi, la scelta cambia il numero riportato. Una proprietà da ricordare del MAP: dipende da come parametrizzi il problema. Se stimi θ oppure il suo logaritmo dell'odds ottieni due mode che non si corrispondono, mentre la mediana resta la stessa sotto ogni trasformazione monotona (la media no).

MAP e regolarizzazione: la stessa formula vista da due parti

In una regressione lineare y = Xβ + rumore, con rumore gaussiano di deviazione σ, cercare il MAP quando ogni coefficiente ha prior gaussiano N(0, τ²) equivale a minimizzare la somma dei quadrati più una penalità sui coefficienti, cioè la regressione ridge, con parametro λ = σ² / τ². Il prior spiega perché la penalità funziona: dice che i coefficienti grandi sono a priori poco plausibili. Il programma verifica l'uguaglianza numericamente. Calcola la soluzione chiusa della ridge e poi cerca il massimo del log-posterior con un ottimizzatore generico (BFGS) che non conosce quella formula. Se la teoria è giusta, i due risultati coincidono. Per la regolarizzazione dal lato del machine learning vedi Overfitting e regolarizzazione.

import numpy as np
from scipy import optimize

rng = np.random.default_rng(7)
n, p = 40, 3
X = rng.standard_normal((n, p))
beta_vero = np.array([2.0, -1.0, 0.0])
sigma, tau = 1.5, 1.0                       # rumore dei dati e deviazione del prior N(0, tau^2)
y = X @ beta_vero + sigma * rng.standard_normal(n)

lam = sigma**2 / tau**2                     # il legame ridge <-> prior
ridge = np.linalg.solve(X.T @ X + lam * np.eye(p), X.T @ y)

def meno_log_posterior(b):                  # a meno di costanti
    return 0.5 * np.sum((y - X @ b) ** 2) / sigma**2 + 0.5 * np.sum(b**2) / tau**2

map_num = optimize.minimize(meno_log_posterior, np.zeros(p), method="BFGS", options={"gtol": 1e-10}).x
ols = np.linalg.lstsq(X, y, rcond=None)[0]
print("OLS (MLE)  ", np.round(ols, 4))
print("ridge      ", np.round(ridge, 4))
print("MAP numerico", np.round(map_num, 4), " scarto max", f"{np.abs(ridge - map_num).max():.1e}")

# Output:
# OLS (MLE)   [ 2.0581 -0.9495  0.6028]
# ridge       [ 1.9052 -0.896   0.5908]
# MAP numerico [ 1.9052 -0.896   0.5908]  scarto max 1.2e-08

I due vettori coincidono fino alla quarta cifra mostrata e lo scarto massimo è 1,2·10^-8, cioè la tolleranza dell'ottimizzatore. Si vede anche l'effetto del prior: tutti e tre i coefficienti vengono tirati verso lo zero rispetto ai minimi quadrati (2,0581 diventa 1,9052, 0,6028 diventa 0,5908). La scelta di τ è la scelta dell'intensità della penalità: un prior stretto (τ piccolo) significa λ grande e coefficienti più compressi. Questa equivalenza vale per il MAP, non per il posterior intero: la ridge restituisce solo il punto più alto e perde l'incertezza attorno a esso.

Che cosa promette davvero un intervallo di credibilità

Un intervallo di credibilità al 95% è un insieme di valori di θ a cui il posterior assegna probabilità 0,95. La lettura «θ è lì dentro con probabilità 0,95» è legittima, ma è condizionata al prior scelto: la garanzia non riguarda tutti i valori possibili di θ allo stesso modo. Il test seguente lo mostra. Nel primo caso θ viene pescato dal prior Beta(2, 2) stesso e poi si simulano i dati: l'intervallo copre il valore vero nel 95,0% dei casi, come promesso. Nei tre casi successivi θ è fissato e il prior è sbagliato per quel θ: la copertura diventa 92,4% se θ = 0,05, 96,8% se θ = 0,20, 95,7% se θ = 0,50. L'errore di simulazione è circa 0,001, quindi gli scarti dal 95% sono reali.

import numpy as np
from scipy import stats

rng = np.random.default_rng(99)
n = 20

# --- 1) Che cosa promette un intervallo di credibilita' ---------------------
a0, b0 = 2, 2                               # prior Beta(2, 2)
def copertura(thetas, k):
    lo = stats.beta.ppf(0.025, a0 + k, b0 + n - k)
    hi = stats.beta.ppf(0.975, a0 + k, b0 + n - k)
    return np.mean((lo <= thetas) & (thetas <= hi))

R = 100_000
th = rng.beta(a0, b0, R)                    # theta pescato dal prior stesso
print(f"copertura se theta ~ prior   : {copertura(th, rng.binomial(n, th)):.4f}  (atteso 0.95, errore ~0.0007)")
for t_fisso in (0.05, 0.20, 0.50):
    th = np.full(R, t_fisso)
    print(f"copertura se theta = {t_fisso:.2f}   : {copertura(th, rng.binomial(n, th)):.4f}")

# --- 2) Quanto pesa il prior: n=10, k=7 --------------------------------------
print("prior        media   ETI 95%")
for a, b in [(1, 1), (0.5, 0.5), (10, 10), (1, 10)]:
    d = stats.beta(a + 7, b + 3)
    print(f"Beta({a},{b})".ljust(12), f"{d.mean():.3f}  [{d.ppf(.025):.3f}, {d.ppf(.975):.3f}]")

# Output:
# copertura se theta ~ prior   : 0.9502  (atteso 0.95, errore ~0.0007)
# copertura se theta = 0.05   : 0.9241
# copertura se theta = 0.20   : 0.9676
# copertura se theta = 0.50   : 0.9571
# prior        media   ETI 95%
# Beta(1,1)    0.667  [0.390, 0.891]
# Beta(0.5,0.5) 0.682  [0.394, 0.907]
# Beta(10,10)  0.567  [0.389, 0.736]
# Beta(1,10)   0.381  [0.191, 0.592]

Questo contrasto con l'intervallo di confidenza frequentista è il punto da non confondere. L'intervallo di confidenza garantisce una frequenza di copertura qualunque sia θ, ma non dà alcuna probabilità sul θ di questo esperimento. L'intervallo di credibilità dà quella probabilità, ma con la copertura frequentista che dipende da θ e dal prior. La seconda parte del programma mostra quanto pesa il prior quando i dati sono pochi: con 7 successi su 10, il prior Beta(1, 1) porta la media a 0,667 e il Beta(10, 10), che «ricorda» venti prove equilibrate, a 0,567; un prior Beta(1, 10), che sosteneva successi rari, la porta a 0,381 e l'intervallo [0,191; 0,592] arriva a escludere 0,7. Con migliaia di osservazioni la differenza tra questi prior svanisce, perché la verosimiglianza domina: è il contenuto del teorema di Bernstein-von Mises, che in forma intuitiva dice che il posterior diventa gaussiano e centrato sull'MLE.

Un test A/B svolto, e la trappola delle sbirciate

Il confronto classico: la variante A ha 84 click su 1200 visite (7,0%), la B ne ha 99 su 1100 (9,0%). Con prior uniforme i posterior sono Beta(85, 1117) e Beta(100, 1002). La domanda «qual è la probabilità che B sia migliore?» è l'integrale di P(θ_B > t) pesato dalla densità di θ_A in t. Il programma lo calcola con un'integrazione numerica accurata (quad) e lo confronta con una simulazione da 400.000 coppie di estrazioni. In coda risponde a un'altra domanda, che gli articoli sul test A/B di solito saltano: se A e B sono identiche e guardi i dati ogni 100 visite per braccio, fermandoti appena P(B > A) supera 0,95, quante volte dichiari una vittoria che non c'è?

import numpy as np
from scipy import stats, integrate

kA, nA, kB, nB = 84, 1200, 99, 1100
pa, pb = stats.beta(1 + kA, 1 + nA - kA), stats.beta(1 + kB, 1 + nB - kB)
# esatto: integro P(thetaB > t) * densita' di thetaA
esatto = integrate.quad(lambda t: pb.sf(t) * pa.pdf(t), 0, 1, epsabs=1e-12, limit=500)[0]
rng = np.random.default_rng(11)
M = 400_000
a, b = pa.rvs(M, random_state=rng), pb.rvs(M, random_state=rng)
d = b - a
print(f"P(B>A) esatta {esatto:.4f}  simulata {np.mean(b > a):.4f} +- {np.sqrt(esatto*(1-esatto)/M):.4f}")
for nome, dist in (("A", pa), ("B", pb)):
    print(f"{nome}: media posterior {dist.mean():.4f}  ETI 95% [{dist.ppf(.025):.4f}, {dist.ppf(.975):.4f}]")
print(f"differenza: media {d.mean():.4f}  ETI 95% [{np.quantile(d, .025):.4f}, {np.quantile(d, .975):.4f}]")

# Fermarsi appena P(B>A) > 0.95: quanto spesso si "vince" per caso? ---
def prob_b_migliore(kA, nA, kB, nB):        # approssimazione normale dei due posterior Beta
    def mv(k, n):
        a_, b_ = 1 + k, 1 + n - k
        return a_ / (a_ + b_), a_ * b_ / ((a_ + b_) ** 2 * (a_ + b_ + 1))
    (ma, va), (mb, vb) = mv(kA, nA), mv(kB, nB)
    return stats.norm.cdf((mb - ma) / np.sqrt(va + vb))

rng = np.random.default_rng(2027)
p_vero, passo, tappe, ripetizioni = 0.07, 100, 20, 4000   # A e B IDENTICHE
falsi_con_sbirciate = falsi_a_fine = 0
for _ in range(ripetizioni):
    ca = np.cumsum(rng.binomial(passo, p_vero, tappe))
    cb = np.cumsum(rng.binomial(passo, p_vero, tappe))
    n_cum = passo * np.arange(1, tappe + 1)
    pr = prob_b_migliore(ca, n_cum, cb, n_cum)
    falsi_con_sbirciate += np.any(pr > 0.95)
    falsi_a_fine += pr[-1] > 0.95
e1, e2 = falsi_con_sbirciate / ripetizioni, falsi_a_fine / ripetizioni
print(f"vittorie false, 20 sbirciate: {e1:.3f} +- {np.sqrt(e1*(1-e1)/ripetizioni):.3f}")
print(f"vittorie false, un solo sguardo a 2000 visite per braccio: {e2:.3f} +- {np.sqrt(e2*(1-e2)/ripetizioni):.3f}")

# Output:
# P(B>A) esatta 0.9615  simulata 0.9607 +- 0.0003
# A: media posterior 0.0707  ETI 95% [0.0569, 0.0859]
# B: media posterior 0.0907  ETI 95% [0.0745, 0.1084]
# differenza: media 0.0200  ETI 95% [-0.0023, 0.0425]
# vittorie false, 20 sbirciate: 0.197 +- 0.006
# vittorie false, un solo sguardo a 2000 visite per braccio: 0.049 +- 0.003

La probabilità che B batta A è 0,9615, e la simulazione dà 0,9607: lo scarto è 0,0008, circa 2,6 errori standard, quindi dentro i 3 che usiamo come soglia (con altri tre semi provati lo scarto è rimasto fra 0,0000 e 0,0007; l'integrale numerico è il valore di riferimento). La differenza tra i due tassi ha media 0,0200 e intervallo di credibilità al 95% [−0,0023; 0,0425]: comprende lo zero. Non c'è contraddizione con il 96%: P(B > A) = 0,96 significa che lo zero sta nel 4% di sinistra, mentre un intervallo bilaterale al 95% ne lascia fuori il 2,5% per lato. Chi legge «96%» e annuncia che il nuovo pulsante vince con certezza sta dicendo più di quanto i dati sostengano.

L'ultima riga è la correzione a un luogo comune. È vero che il posterior, calcolato su tutti i dati visti fino a un momento, resta una descrizione coerente dell'incertezza anche se hai deciso quando fermarti guardandolo. Ma una regola di arresto come «fermati quando P(B > A) > 0,95» produce vittorie false con una frequenza che dipende da quante volte guardi: nella simulazione sono il 19,7% (errore 0,6%) con venti sbirciate, contro il 4,9% (errore 0,3%) se guardi una sola volta a 2000 visite per braccio. Per una decisione che costa, fissa in anticipo la soglia e le dimensioni, oppure verifica col simulatore le frequenze di errore della tua regola. Per il disegno degli esperimenti vedi Test and learn: esperimenti controllati e A/B test.

Quando il posterior non si calcola: Metropolis

Fuori dai modelli coniugati l'integrale dell'evidenza è intrattabile, ma per campionare dal posterior non serve conoscerlo: basta poterlo valutare a meno di una costante. L'algoritmo di Metropolis (1953) costruisce una catena di Markov, cioè una sequenza in cui ogni valore dipende solo dal precedente, che ha il posterior come distribuzione di lungo periodo (le catene sono l'oggetto dell'articolo successivo). A ogni passo propone un nuovo valore sommando un rumore gaussiano di ampiezza passo, e lo accetta con probabilità min(1, rapporto dei posterior); se lo rifiuta, ripete il valore corrente. Qui il bersaglio è la Beta(9, 13) già nota, così puoi confrontare l'esito con la risposta esatta. Il listato mostra anche le due diagnosi che servono sempre: R̂ di Gelman e Rubin (1992), che confronta la varianza dentro e fra quattro catene partite da punti diversi e vale circa 1 se hanno raggiunto la stessa distribuzione, e la dimensione campionaria effettiva (ESS), il numero di estrazioni indipendenti equivalenti alle estrazioni correlate prodotte.

import numpy as np
from scipy import stats

# Metropolis a passo casuale sul posterior Beta(9, 13) -------------------
def log_post(t):                            # log di theta^8 (1-theta)^12, -inf fuori da (0,1)
    return np.where((t > 0) & (t < 1), 8*np.log(np.clip(t, 1e-300, 1)) + 12*np.log(np.clip(1 - t, 1e-300, 1)), -np.inf)

def metropolis(passo, start, T, seed):
    r = np.random.default_rng(seed)
    x = np.empty(T); x[0] = start; acc = 0
    lp = float(log_post(x[0]))
    for i in range(1, T):
        prop = x[i-1] + passo * r.standard_normal()
        lpp = float(log_post(prop))
        if np.log(r.random()) < lpp - lp:   # accetto con prob. min(1, rapporto)
            x[i], lp, acc = prop, lpp, acc + 1
        else:
            x[i] = x[i-1]
    return x, acc / (T - 1)

def ess(x):                                 # somma delle autocorrelazioni finche' restano positive
    x = x - x.mean(); T = len(x)
    f = np.fft.rfft(x, 2*T); ac = np.fft.irfft(f*np.conj(f))[:T]; ac /= ac[0]
    s = 0.0
    for lag in range(1, T - 1, 2):          # coppie di ritardi, regola di Geyer
        g = ac[lag] + ac[lag+1]
        if g < 0: break
        s += g
    return T / (1 + 2*s)

def rhat(catene):                           # Gelman-Rubin classico, catene di pari lunghezza
    m, T = catene.shape
    W = catene.var(axis=1, ddof=1).mean()
    B = T * catene.mean(axis=1).var(ddof=1)
    return np.sqrt(((T-1)/T * W + B/T) / W)

esatto = stats.beta(9, 13)
T, burn = 20_000, 2_000
print("passo   accettazione  media    sd      ESS   R-hat")
for passo in (0.005, 0.3, 5.0):
    cat = np.array([metropolis(passo, s, T, 100 + j)[0][burn:] for j, s in enumerate((0.05, 0.3, 0.7, 0.95))])
    acc = np.mean([metropolis(passo, 0.3, T, 100 + j)[1] for j in range(4)])
    print(f"{passo:<7} {acc:>10.3f}  {cat.mean():.4f}  {cat.std():.4f}  {sum(ess(c) for c in cat):>6.0f}  {rhat(cat):.3f}")
print(f"esatto           {esatto.mean():.4f}  {esatto.std():.4f}")

# Output:
# passo   accettazione  media    sd      ESS   R-hat
# 0.005        0.985  0.4129  0.1036      66  1.066
# 0.3          0.388  0.4076  0.1019   16934  1.000
# 5.0          0.027  0.4099  0.1019    1245  1.001
# esatto           0.4091  0.1025

Il passo 0,3 funziona: accettazione 38,8%, media 0,4076 e deviazione 0,1019 contro 0,4091 e 0,1025 esatti, R̂ = 1,000, ESS circa 16.900 su 72.000 campioni (quattro catene da 18.000 dopo lo scarto dei primi 2.000). Il passo 0,005 accetta quasi tutto (98,5%) ma si muove a piccoli scatti: ESS 66, R̂ 1,066 e media 0,4129 sbagliata, perché le catene non hanno ancora esplorato la distribuzione. Il passo 5,0 propone salti enormi che quasi sempre cadono dove il posterior è trascurabile e vengono rifiutati (accettazione 2,7%): la catena resta ferma e l'ESS scende a circa 1.250. Il lavoro sul passo è quindi un compromesso fra accettare spesso e muoversi lontano. Per modelli con molti parametri si usano versioni più efficienti, come l'Hamiltonian Monte Carlo implementato in Stan o PyMC, non eseguite qui perché le librerie non sono nella macchina di prova: il principio, e le diagnosi, sono quelli mostrati.

Esercizio svolto e prova da fare da solo

Un laboratorio vuole stimare la frazione di pezzi difettosi. Con prior Beta(2, 18) (media 0,10, equivalente a 20 pezzi passati) osserva 3 difetti in 25 pezzi. Il posterior è Beta(2 + 3, 18 + 22) = Beta(5, 40): media 5/45 = 0,111, MAP (5 − 1)/(5 + 40 − 2) = 4/43 = 0,0930. L'MLE sarebbe 3/25 = 0,12. Il posterior pesa 25/45 = 56% sui dati e 44% sul prior: media = 0,56 · 0,12 + 0,44 · 0,10 = 0,111, come previsto.

Che cosa sai fare adesso

Sai trasformare un prior e dei dati in un posterior coniugato, distinguere MAP, mediana e media sapendo a quale perdita corrisponde ciascuna, e leggere un intervallo di credibilità per ciò che promette e non promette. Sai collegare il MAP gaussiano alla ridge, calcolare P(B > A) in un test A/B senza prenderla per una certezza, e sai che fermarsi alla prima soglia superata gonfia le vittorie false. Quando la forma chiusa manca, sai scrivere un Metropolis e controllarne passo, R̂ ed ESS. Il prossimo articolo, sulle catene di Markov, spiega perché quel campionatore converge; il bootstrap mostra l'alternativa frequentista per costruire intervalli senza prior. Per la versione in ottica machine learning vedi Statistica bayesiana: prior, posterior, coniugati. L'articolo precedente è la regressione lineare.

Nella stessa serie trovi anche Analisi dati esplorativa (EDA) con Python.