44  Valutare i modelli bayesiani

Introduzione

I capitoli precedenti ci hanno fornito due strumenti concettuali fondamentali: l’entropia, che misura l’incertezza di una distribuzione, e la divergenza di Kullback–Leibler (\(D_{\text{KL}}\)), che quantifica la distanza informazionale tra due distribuzioni. È giunto il momento di fare il passo successivo: utilizzare questi strumenti per valutare e confrontare i modelli statistici in ambito bayesiano.

Tutto parte da una domanda cruciale: quanto sarà capace il nostro modello di prevedere dati mai visti? Un modello valido non è quello che si adatta perfettamente ai dati osservati, ma quello che sa generalizzare efficacemente a situazioni future. Questa distinzione tra adattamento (fitting) e generalizzazione è il cuore della valutazione predittiva.

Immaginiamo di sviluppare un test psicometrico per stimare l’ansia pre-esame. Non ci interessa solo che il modello descriva bene il campione usato per costruirlo; vogliamo essere sicuri che fornirà previsioni affidabili anche per nuovi studenti. Scegliere tra modelli statistici, in questo senso, è come scegliere tra test psicologici: cerchiamo lo strumento che offre la maggiore affidabilità predittiva.

Un approccio intuitivo per valutare le prestazioni predittive è calcolare l’errore quadratico medio (MSE)¹ sui dati a disposizione, presupponendo che sarà simile per dati futuri. Tuttavia, questa strategia trascura un elemento essenziale: l’incertezza delle stime. Le nostre previsioni dipendono da parametri di cui non conosciamo il valore esatto, ma solo la distribuzione a posteriori. È la distribuzione predittiva a posteriori che, propagando questa incertezza, fornisce una visione completa e probabilistica delle previsioni.

Sorge però un problema: valutare la predittiva sugli stessi dati usati per stimare il modello porta a un ottimismo eccessivo e al rischio di overfitting. Per una stima realistica della capacità di generalizzazione, serve un criterio che valuti la performance predittiva su dati non osservati.

È qui che entra in gioco l’Expected Log Predictive Density (ELPD). L’ELPD misura, attraverso il log-score, quanto bene un modello prevede dati futuri generati dalla distribuzione vera. In sostanza, rappresenta l’accuratezza predittiva attesa. Un ELPD più alto indica una migliore capacità predittiva. (Anticipiamo un collegamento teorico che formalizzeremo dopo: massimizzare l’ELPD equivale a minimizzare la divergenza di Kullback-Leibler tra il modello e il processo generativo dei dati sottostante.)

La difficoltà pratica risiede nel fatto che l’ELPD, definito come un’aspettativa rispetto alla distribuzione vera – e ignota – dei dati, non è calcolabile in modo diretto. È necessario pertanto ricorrere a delle approssimazioni numeriche. La più comune e teoricamente solida è la cross-validazione leave-one-out (LOO-CV). Questo metodo valuta, per ogni singola osservazione, la capacità predittiva del modello quando questa viene esclusa dal set di addestramento. Sommando i contributi log-verosimili di tutte le osservazioni così valutate, si ottiene una stima robusta dell’ELPD out-of-sample, notevolmente più resistente al fenomeno dell’overfitting rispetto a stime basate sui dati di addestramento.

Questo capitolo illustra come la cross-validazione leave-one-out (LOO) fornisca un criterio rigoroso per il confronto tra modelli bayesiani, consentendo di mitigare il rischio di sovradattamento e orientando la selezione verso i modelli dotati delle migliori capacità predittive generalizzabili. Attraverso un itinerario concettuale progressivo – dal log-score come misura di accuratezza, alla LPPD (log-posterior predictive density) come stima in-sample, fino alla stima dell’ELPD (Expected Log Predictive Density) mediante PSIS-LOO – mostreremo il fondamento di questi strumenti nella teoria dell’informazione e la loro applicazione pratica nel confronto tra modelli.

Panoramica del capitolo

Per valutare la qualità predittiva dei modelli, seguiremo un percorso strutturato che dalla teoria conduce alla pratica applicativa:

1. Fondamenti della valutazione predittiva

• Distribuzione predittiva a posteriori: propagare l’incertezza parametrica nelle previsioni.
• Log-score: valutare l’accuratezza predittiva per ogni singola osservazione.

2. Dalla teoria alla stima pratica

• Da LPPD a ELPD: distinguere la bontà di adattamento (in-sample) dalla capacità di generalizzazione (out-of-sample).
• Fondamento teorico: collegare l’ELPD alla minimizzazione della divergenza di Kullback-Leibler.

3. Implementazione e confronto

• Stima con PSIS-LOO: approssimare efficientemente l’ELPD, con diagnostiche per l’affidabilità.
• Confronto tra modelli: utilizzare ΔELPD e il suo errore standard per identificare il modello con le migliori prestazioni predittive, bilanciando prestazioni e parsimonia.

ConsiglioPrerequisiti

• Per comprendere appieno questo capitolo è utile leggere il capitolo 7 Ulysses’ Compass di Statistical Rethinking (McElreath (2020)).

AttenzionePreparazione del Notebook

here::here("code", "_common.R") |> 
  source()

library(gt)
library(conflicted)
library(brms)
library(loo)
conflicts_prefer(rstan::loo)

ConsiglioIn sintesi

ELPD in pillole

• Definizione: Misura della capacità predittiva out-of-sample di un modello.
• Calcolo teorico: Media del log-score atteso per nuove osservazioni. Valori più alti indicano modelli più predittivi.
• Fondamento teorico: Massimizzare l’ELPD equivale a minimizzare la divergenza di Kullback-Leibler dalla distribuzione vera dei dati.
• Applicazione pratica: Stimato efficientemente mediante cross-validazione leave-one-out (PSIS-LOO), senza necessità di nuovi campioni.

Dal problema alla metrica: predizione bayesiana e log-score {#sec-predizione-logscore}

ImportanteIntuizione e definizione dell’ELPD

L’obiettivo: prevedere il futuro

Immaginiamo che il nostro modello, addestrato sui dati \(y\), produca una distribuzione predittiva \(q(\tilde{y} \mid y)\). L’obiettivo non è solo adattarsi ai dati osservati, ma predire bene dati nuovi \(\tilde{y}\), generati dalla distribuzione vera e sconosciuta \(p(\tilde{y})\).

La domanda chiave diventa: quanto è efficace \(q\) nel prevedere \(p\)?

La metrica: il log-score

La risposta risiede nel log-score, uno scoring rule proprio che assegna il punteggio massimo quando \(q\) coincide con \(p\).

Definiamo quindi l’Expected Log Predictive Density (ELPD): \[ \mathrm{ELPD} = \mathbb{E}_{p(\tilde y)}\!\left[\log q(\tilde y \mid y)\right]. \] Questa quantità rappresenta il valore atteso del log-score che il modello assegna a una nuova osservazione \(\tilde{y}\) generata da \(p\).

Regola fondamentale: massimizzare l’ELPD equivale a ridurre al minimo la distanza informazionale dalla verità.

Giustificazione teorica

Il legame con la teoria dell’informazione è esplicito attraverso la divergenza di Kullback–Leibler (KL):

\[ D_{\mathrm{KL}}(p \parallel q) \;=\; \underbrace{\mathbb{E}_{p}[\log p(\tilde y)]}_{\text{entropia di } p} \;-\; \underbrace{\mathbb{E}_{p}[\log q(\tilde y \mid y)]}_{\text{ELPD}}. \]

Dato che l’entropia di \(p\) è una costante intrinseca ai dati, minimizzare la divergenza KL equivale direttamente a massimizzare l’ELPD. Scegliere il modello con l’ELPD più alto significa quindi selezionare la distribuzione predittiva più vicina alla verità in senso informazionale.

44.1 Dalla Teoria alla stima

NotaDalla teoria alla pratica

La sfida pratica e la soluzione (LOO)

Il calcolo diretto dell’ELPD è impossibile perché \(p(\tilde{y})\) è sconosciuta. La soluzione è stimarlo dai dati disponibili, evitando di usare la stessa osservazione sia per l’addestramento che per il test.

Il metodo di riferimento è la Leave-One-Out Cross-Validation (LOO):

\[ \widehat{\mathrm{ELPD}}_{\mathrm{LOO}} = \sum_{i=1}^{n} \log q(y_i \mid y_{-i}), \] dove \(q(y_i \mid y_{-i})\) è la densità predittiva per \(y_i\) quando il modello è addestrato su tutti i dati tranne \(y_i\).

Questo approccio fornisce una stima quasi non distorta dell’ELPD, concettualmente identica a quella teorica ma calcolabile senza conoscere la distribuzione vera.

ConsiglioEfficienza computazionale: PSIS-LOO

Ricalcolare il modello \(n\) volte è computazionalmente proibitivo. La soluzione è PSIS-LOO (Pareto Smoothed Importance Sampling LOO), che approssima efficientemente la stima LOO senza dover rifittare il modello per ogni osservazione.

In aggiunta, PSIS-LOO fornisce diagnostiche Pareto-\(k\) per valutare l’affidabilità dell’approssimazione e identificare osservazioni eccessivamente influenti.

44.1.0.1 Interpretazione e confronto dell’ELPD

ConsiglioInterpretazione e confronto

Come interpretare l’ELPD

• Scala: L’ELPD è una somma su tutte le \(n\) osservazioni. Valori più alti indicano un modello migliore.
• Confronto: Per due modelli \(A\) e \(B\), si calcola \(\Delta ELPD = ELPD_A - ELPD_B\). Il suo errore standard aiuta a valutare la credibilità del miglioramento.
• Normalizzazione: Per dataset di dimensioni diverse, utilizzare l’ELPD per osservazione (\(ELPD/n\)) per confronti equi.

Nota sulla devianza: Alcuni software (es. Stan) riportano la devianza (\(-2 \times ELPD\)). In questo caso, valori più bassi indicano un modello migliore. Verificare sempre la convenzione utilizzata.

NotaEsempio: Confronto tra due modelli

Supponiamo di confrontare due modelli, \(A\) e \(B\), sullo stesso dataset (\(n = 200\)).

Modello	ELPD totale	ELPD/osservazione	ΔELPD vs B	SE(ΔELPD)
A	-520	-2.60	+10	4.5
B	-530	-2.65	—	—

Interpretazione:

• Il modello A ha un ELPD più alto.
• ΔELPD = 10, con SE = 4.5. Poiché la differenza è circa 2 volte l’errore standard, il miglioramento di A rispetto a B è credibile.

NotaEsempio con la convenzione della Devianza (Stan)

Alcuni software, come Stan, riportano la devianza: \(\text{Devianza} = -2 \times ELPD\).

Modello	ELPD totale	Devianza	ΔELPD vs B	SE(ΔELPD)
A	-520	1040	+10	4.5
B	-530	1060	—	—

Interpretazione (convenzione Stan):

• Un valore di devianza più basso indica un modello migliore.
• Il modello A (devianza 1040) è preferibile al modello B (devianza 1060).
• La differenza di devianza (20 punti) conferma un ΔELPD = 10, con SE = 4.5.

44.2 Dal problema alla metrica: la predizione bayesiana

Prima di addentrarci nei dettagli tecnici, facciamo un passo indietro per chiarire perché la valutazione predittiva bayesiana è diversa da quella frequentista classica.

Nell’approccio classico, una volta stimati i parametri (ad esempio con la massima verosimiglianza), li trattiamo come valori “veri” e fissi. Se \(\hat{\theta}\) è la nostra stima puntuale, la previsione per un nuovo dato \(\tilde{y}\) è:

\[ p(\tilde{y} \mid \hat{\theta}). \]

In questo modo si ragiona come se \(\hat{\theta}\) fosse noto con certezza. Ma nella realtà la stima è solo una “fotografia” dei dati raccolti: se ripetessimo l’esperimento, otterremmo una stima diversa. Questo approccio, quindi, ignora l’incertezza sui parametri.

Nell’approccio bayesiano, invece, riconosciamo che i parametri non sono numeri fissi ma quantità sconosciute. Dopo aver osservato i dati, la nostra conoscenza su \(\theta\) è descritta da un’intera distribuzione — la distribuzione a posteriori \(p(\theta \mid y)\) — che riflette tutte le ipotesi plausibili sui parametri date le osservazioni.

La previsione per un nuovo dato \(\tilde{y}\) diventa quindi:

\[ p(\tilde{y} \mid y) = \int p(\tilde{y} \mid \theta)\, p(\theta \mid y)\, d\theta. \tag{44.1}\]

Questa è la distribuzione predittiva a posteriori: una media pesata delle previsioni possibili, dove il peso è dato da quanto ciascun valore di \(\theta\) è coerente con i dati osservati.

NotaAnalogia didattica

Immaginiamo di dover prevedere il tempo di domani.

• L’approccio “classico” è come chiedere a un solo meteorologo, prendere la sua previsione e considerarla come se fosse la verità.
• L’approccio “bayesiano” è come consultare molti meteorologi diversi, dare più fiducia a chi ha fatto buone previsioni in passato, e costruire una previsione che tenga conto di tutte le opinioni plausibili.

Il risultato non è una previsione unica, ma una distribuzione di possibili scenari, che ci dice sia cosa aspettarci più probabilmente sia quanto siamo incerti.

44.2.1 La distribuzione predittiva posteriore

Abbiamo già incontrato la distribuzione predittiva posteriore nel capitolo sul modello Beta–Binomiale: è lo strumento che, nell’approccio bayesiano, consente di prevedere nuovi dati incorporando sia la struttura del modello sia l’incertezza sui parametri.

La logica è elegante nella sua semplicità: dopo aver osservato i dati \(y\), non otteniamo un singolo “miglior” valore dei parametri, ma una distribuzione posteriore \(p(\theta \mid y)\) che riflette i valori plausibili di \(\theta\) e la nostra incertezza.

Esempio concreto:
Uno psicologo che stima il livello medio di ansia in una popolazione, invece di affermare “la media è 4.7”, dirà: “il valore più plausibile è 4.7, ma è ragionevole che sia tra 4.2 e 5.1”, riflettendo la variabilità della distribuzione a posteriori.

44.2.1.1 Comprendere l’integrale predittivo

L’Equazione 44.1 può sembrare astratta, ma si capisce bene se la leggiamo passo passo:

1. \(p(\tilde{y} \mid \theta)\): se conoscessimo il vero valore di \(\theta\), questa sarebbe la nostra previsione per \(\tilde{y}\).
   
2. \(p(\theta \mid y)\): descrive la nostra incertezza su quale sia il vero valore di \(\theta\).
   
3. L’integrale: combina le previsioni possibili, pesandole secondo la plausibilità di ciascun valore di \(\theta\).

NotaIntuizione dettagliata

Se conoscessimo \(\theta\), la predizione sarebbe:

\[ p(\tilde y \mid \theta). \] Ma \(\theta\) non lo conosciamo: abbiamo solo la distribuzione posteriore \(p(\theta \mid y)\).
La predittiva posteriore si ottiene quindi come media pesata:

\[ p(\tilde y \mid y) = \int p(\tilde y \mid \theta)\, p(\theta \mid y)\, d\theta. \]

Esempio binomiale.
Se i dati futuri sono generati da una Binomiale con \(m\) prove e parametro \(\theta\):

\[ p(\tilde y = x \mid \theta) = \binom{m}{x}\,\theta^x (1-\theta)^{m-x}, \]

allora la distribuzione predittiva diventa:

\[ p(\tilde y = x \mid y) = \int \binom{m}{x}\,\theta^x (1-\theta)^{m-x}\,p(\theta\mid y)\,d\theta. \]

• Con prior Beta e dati binomiali, questo integrale ha una soluzione chiusa: la distribuzione Beta–Binomiale.
  
• In generale, non c’è formula chiusa: servono approssimazioni numeriche (ottenute mediante il metodo basato su griglia o mediante il metodo MCMC).

44.2.1.2 Il metodo su griglia

Nella maggior parte dei casi l’integrale che definisce la distribuzione predittiva posteriore non può essere risolto in forma chiusa. Un modo semplice e intuitivo per affrontare il problema è approssimarlo con una somma, utilizzando una griglia di valori del parametro \(\theta\). L’idea di fondo è trasformare un calcolo continuo (l’integrale) in un calcolo discreto (una somma pesata).

Il procedimento segue quattro passaggi:

1. Definire una griglia di valori possibili per \(\theta\). Ad esempio, scegliamo \(J=1000\) punti equispaziati tra 0 e 1: \[ \theta_1, \theta_2, \dots, \theta_J. \]

2. Calcolare la distribuzione posteriore su ogni punto della griglia. Per ciascun \(\theta_j\), valutiamo la densità posteriore \(p(\theta_j \mid y)\) e normalizziamo in modo che la somma dei valori sia uguale a 1. I risultati diventano i pesi \(w_j\): \[ w_j = \frac{p(\theta_j \mid y)}{\sum_{\ell=1}^J p(\theta_\ell \mid y)}. \]

3. Combinare le distribuzioni predittive condizionate. Per ogni possibile esito futuro \(x\), calcoliamo la probabilità predittiva come media pesata delle predittive condizionate: \[ p(\tilde y = x \mid y) \;\approx\; \sum_{j=1}^J w_j \, p(\tilde y = x \mid \theta_j). \]

4. Ottenere la distribuzione predittiva. Ripetendo il calcolo per tutti i possibili valori di \(x\), otteniamo una distribuzione completa (pmf) per i dati futuri. Da questa distribuzione possiamo:

   • ricavare probabilità,
   • generare campioni simulati di \(\tilde y\),
   • confrontare le previsioni del modello con i dati osservati.

Da ricordare:

• La distribuzione predittiva posteriore non è la media di singoli valori, ma una vera e propria distribuzione di probabilità sui dati futuri.
• Il metodo su griglia è il modo più semplice per implementare questa logica: discretizza i possibili valori di \(\theta\), li pesa in base alla loro plausibilità e combina le predizioni corrispondenti.
• In problemi più complessi, la stessa idea viene realizzata in maniera più efficiente con l’MCMC: invece di una griglia fissa, utilizziamo campioni \(\theta^{(s)}\) estratti dalla distribuzione posteriore.

NotaEsempio numerico completo (R)

# ------------------------------------------------------------------
# Predittiva posteriore per Binomiale con metodo su griglia 
# ------------------------------------------------------------------

# 1) Dati osservati e prior
k <- 10   # successi osservati
n <- 50   # prove osservate totali
a <- 1    # prior Beta(a, b)
b <- 1
m <- 10   # numero di prove future che vogliamo prevedere (per costruire la pmf di tilde y)
J <- 1001 # numero di punti della griglia su theta (dispari per includere 0, 0.5, 1)

# 2) Griglia su theta (da 0 a 1)
theta <- seq(0, 1, length.out = J)

# 3) Densità posteriore NON normalizzata su ogni punto di griglia
#    Posteriore Beta(a + k, b + n - k) -> densità proporzionale a theta^(a+k-1) * (1-theta)^(b+n-k-1)
post_unnorm <- theta^(a + k - 1) * (1 - theta)^(b + n - k - 1)

# 4) Normalizziamo per ottenere pesi che sommano a 1 (w_j)
w <- post_unnorm / sum(post_unnorm)

# (controllo) I pesi devono sommare a 1
sum(w)  # ~ 1
#> [1] 1

# 5) Costruiamo la pmf predittiva p(tilde y = x | y) come media pesata delle Binomiali condizionate
x_vals  <- 0:m                # tutti i possibili esiti futuri in m prove
pred_pmf <- rep(NA_real_, m+1)

# Per ogni x, sommiamo w_j * P(tilde y = x | theta_j)
for (i in 1:(m+1)) {
  x <- x_vals[i]
  # vettore di probabilità condizionate per tutti i theta_j
  p_x_given_theta <- dbinom(x, size = m, prob = theta)
  # media pesata rispetto ai pesi posteriori w_j
  pred_pmf[i] <- sum(w * p_x_given_theta)
}

# (controllo) La pmf deve sommare a 1
sum(pred_pmf)  # ~ 1
#> [1] 1

# 6) Tabella finale con la pmf predittiva
pred_df <- data.frame(x = x_vals, p = pred_pmf)
head(pred_df)   # prime righe
#>   x      p
#> 1 0 0.1139
#> 2 1 0.2506
#> 3 2 0.2762
#> 4 3 0.1995
#> 5 4 0.1040
#> 6 5 0.0407

Confronto con la formula chiusa (Beta–Binomiale) per verifica didattica:

# 7) Parametri della posteriore coniugata
a_post <- a + k
b_post <- b + n - k

# 8) Funzione Beta-Binomiale (pmf): choose(m, x) * Beta(x+a_post, m-x+b_post) / Beta(a_post, b_post)
dbetabinom <- function(x, m, a, b) {
  choose(m, x) * beta(x + a, m - x + b) / beta(a, b)
}

# 9) pmf teorica Beta-Binomiale
bb_pmf <- rep(NA_real_, m+1)
for (i in 1:(m+1)) {
  bb_pmf[i] <- dbetabinom(x_vals[i], m, a_post, b_post)
}

# 10) Confronto numerico: scarto massimo tra griglia e formula chiusa
max_abs_err <- max(abs(pred_pmf - bb_pmf))
max_abs_err   # con J grande dovrebbe essere molto vicino a 0
#> [1] 1.92e-15

# Affianca le due colonne per ispezione
data.frame(x = x_vals, grid = pred_pmf, beta_binom = bb_pmf) |> head()
#>   x   grid beta_binom
#> 1 0 0.1139     0.1139
#> 2 1 0.2506     0.2506
#> 3 2 0.2762     0.2762
#> 4 3 0.1995     0.1995
#> 5 4 0.1040     0.1040
#> 6 5 0.0407     0.0407

Cosa osservare:

• w contiene i pesi posteriori sui diversi valori di \(\theta\); sommano a 1.
• Per ogni possibile esito futuro x, la pmf predittiva pred_pmf[i] è una media pesata delle probabilità binomiali condizionate a ciascun theta[j].
• Con prior Beta e dati binomiali, la soluzione chiusa (Beta–Binomiale) coincide con l’approssimazione su griglia (errore massimo vicino a 0 se J è grande).

44.2.1.3 Riepilogo notazionale

Per sintesi useremo talvolta la forma compatta \(p(\cdot \mid y)\) per indicare la predittiva posteriore. Quando vogliamo evidenziare la previsione marginale per una singola osservazione \(y_i\) scriviamo:

\[ p(y_i \mid y) = \int p(y_i \mid \theta)\, p(\theta \mid y)\, d\theta. \]

Idea chiave: la predittiva posteriore propaga l’incertezza sui parametri nelle previsioni. Questo rende la valutazione predittiva coerente con il principio bayesiano e utilizzabile nel confronto tra modelli.

NotaMappa concettuale

Quantità	Significato	Uso principale
\(p(y_i \mid \theta)\)	Verosimiglianza	Calcolo predittivo
\(p(\theta \mid y)\)	Distribuzione posteriore	Ponderazione
\(p(y_i \mid y)\)	Predizione media bayesiana	Log-score, LPPD
\(p(y_i \mid y_{-i})\)	Predizione LOO (leave-one-out)	ELPD
\(p(\tilde{y} \mid y)\)	Distribuzione predittiva completa	Divergenza KL, confronto modelli

44.2.2 Il problema della valutazione predittiva

Ora che sappiamo come costruire previsioni bayesiane, affrontiamo la domanda centrale: come possiamo valutare la qualità di queste previsioni?

44.2.2.1 Dilemma teorico: fitting vs generalizzazione

Quando costruiamo un modello statistico, non ci interessa solo spiegare i dati già osservati, ma vogliamo sapere quanto bene il modello riuscirà a predire dati futuri.

In termini formali, distinguiamo tra due distribuzioni:

• la vera distribuzione dei dati futuri \(p(\tilde{y})\), che purtroppo non conosciamo mai;
• la distribuzione predittiva del modello \(p(\tilde{y} \mid y)\), cioè le previsioni del modello dopo aver visto i dati osservati \(y\).

L’obiettivo è quindi misurare quanto \(p(\tilde{y} \mid y)\) si avvicina a \(p(\tilde{y})\).

44.2.2.2 La divergenza KL come misura di distanza

La divergenza di Kullback–Leibler (KL) è lo strumento che la teoria dell’informazione ci mette a disposizione per quantificare questa distanza:

\[ D_{\text{KL}}(p \parallel q) \;=\; \mathbb{E}_{p} \left[ \log \frac{p(\tilde{y})}{q(\tilde{y} \mid y)} \right]. \tag{44.2}\]

Interpretazione intuitiva:

• se il modello \(q\) assegna probabilità elevate agli stessi eventi che sono probabili nella realtà \(p\), allora \(D_{\text{KL}}\) sarà piccolo;
• se invece il modello “sbaglia bersaglio”, attribuendo alta probabilità a eventi che in realtà sono rari, \(D_{\text{KL}}\) sarà grande.

NotaIdea chiave

La divergenza KL misura quanta informazione perdiamo quando usiamo le predizioni del modello \(q\) al posto della distribuzione vera \(p\).

44.2.2.3 L’ostacolo pratico e la motivazione dell’ELPD

C’è però un problema fondamentale: la distribuzione vera \(p(\tilde{y})\) non la conosciamo mai.

Abbiamo soltanto i dati osservati. Non possiamo quindi calcolare direttamente la distanza tra \(p\) e \(q\). Questo rende necessarie strategie di approssimazione, che vedremo nei prossimi capitoli: la validazione incrociata e criteri predittivi come l’Expected Log Predictive Density (ELPD).

44.2.2.4 Esempio intuitivo: la moneta truccata

Immagina una moneta truccata che produce testa il 70% delle volte. Vogliamo confrontare due modelli:

• Modello A: assume che la moneta sia equa (\(p=0.5\));
• Modello B: assume che la probabilità di testa sia \(p=0.65\).

Se sapessimo che la probabilità “vera” è \(p=0.7\), sarebbe chiaro che il Modello B è più vicino alla realtà. La divergenza di Kullback–Leibler servirebbe proprio a quantificare quanta informazione perdiamo scegliendo A invece di B.

ImportanteIl punto cruciale

Nella pratica non conosciamo mai la probabilità vera della moneta. Abbiamo solo i dati osservati (i lanci). Per decidere quale modello predice meglio, dobbiamo basarci su criteri che stimano indirettamente la capacità predittiva: ELPD e LOO sono tra i più efficaci.

44.3 Dal log-score alla LPPD: accuratezza in-sample

Abbiamo definito la distribuzione predittiva posteriore e visto il problema teorico della valutazione. Ora introduciamo il primo strumento pratico: il log-score.

44.3.0.1 La domanda di base: quanto il modello scommette bene?

Per ogni osservazione nei nostri dati ci chiediamo: quanto il nostro modello considerava plausibile proprio questo risultato? Il log-score fornisce la risposta.

44.3.0.2 Definizione formale del log-score

Il log-score per un’osservazione \(y_i\) è:

\[ \text{Log-score}(y_i) = \log p(y_i \mid y) , \tag{44.3}\] dove \(p(y_i \mid y)\) è la distribuzione predittiva posteriore: \[ p(y_i \mid y) = \int p(y_i \mid \theta)\, p(\theta \mid y)\, d\theta . \]

44.3.0.3 Interpretazione: probabilità e fiducia del modello

Il log-score misura quanto bene il modello avrebbe “puntato” sul risultato osservato:

• se \(p(y_i \mid y)\) è alto (l’osservazione era plausibile), il log-score è vicino a 0 → buono;
• se \(p(y_i \mid y)\) è basso (l’osservazione era quasi esclusa), il log-score diventa molto negativo → scarso.

NotaPerché usare il logaritmo?

Il logaritmo trasforma prodotti di probabilità in somme, rendendo i calcoli più stabili. Inoltre penalizza fortemente i casi in cui il modello assegna probabilità quasi nulla a ciò che accade davvero.

44.3.0.4 Dal singolo dato alla LPPD

Per valutare l’intero dataset, sommiamo i log-score di tutte le osservazioni:

\[ S = \sum_{i=1}^n \log p(y_i \mid y). \tag{44.4}\]

Più \(S\) è alto (meno negativo), più il modello “scommette bene” sui dati osservati (in-sample).

44.3.0.5 Parametri fissi vs incerti: due filosofie

44.3.0.5.1 Approccio classico: parametri fissi

Nell’impostazione frequentista tradizionale fissiamo i parametri a una stima puntuale, tipicamente \(\hat{\theta}\), e calcoliamo:

\[ \log p(y_i \mid \hat{\theta}). \] Questo ignora completamente l’incertezza sui parametri.

44.3.0.5.2 Approccio bayesiano: parametri incerti

Nell’approccio bayesiano riconosciamo che i parametri restano incerti anche dopo aver visto i dati. Perciò mescoliamo le previsioni su tutti i valori plausibili:

\[ \begin{align} \log p(y_i \mid y) &= \log \int p(y_i \mid \theta)\, p(\theta \mid y)\, d\theta . \end{align} \tag{44.5}\]

NotaDifferenza chiave

• Frequentista: “Quanto sono plausibili i dati se i parametri valgono esattamente \(\hat{\theta}\)?”
• Bayesiano: “Quanto sono plausibili i dati, in media, considerando tutti i valori di \(\theta\) compatibili con i dati osservati?”

La seconda prospettiva è più realistica, perché include l’incertezza sui parametri.

44.3.0.6 Implementazione pratica (MCMC)

L’integrale dell’Equazione 44.5 di solito non ha soluzione chiusa. Con i campioni MCMC \(\theta^{(1)},\dots,\theta^{(S)}\), lo stimiamo così:

1. Likelihood per campione \(p(y_i \mid \theta^{(s)})\) per \(s=1,\dots,S\).
2. Media aritmetica \[ p(y_i \mid y) \;\approx\; \frac{1}{S}\sum_{s=1}^S p(y_i \mid \theta^{(s)}). \tag{44.6}\]

Infine si prende il logaritmo per ottenere il log-score.

44.3.0.7 La LPPD come accuratezza in-sample

Ripetendo il procedimento per tutte le osservazioni e sommando otteniamo la Log Pointwise Predictive Density (LPPD):

\[ \text{LPPD} = \sum_{i=1}^n \log \left[ \frac{1}{S} \sum_{s=1}^S p(y_i \mid \theta^{(s)}) \right]. \tag{44.7}\]

• Log-score classico: valuta le predizioni con un solo \(\hat{\theta}\) → ignora l’incertezza.
• LPPD bayesiana: media le predizioni su tutti i parametri plausibili → incorpora l’incertezza in modo naturale.

44.3.0.8 Limiti della LPPD e rischio overfitting

La LPPD è calcolata sugli stessi dati usati per stimare il modello. Un modello molto flessibile può adattarsi perfettamente anche al rumore, ottenendo una LPPD alta in-sample ma poco utile per generalizzare.

Analogia: è come valutare uno studente facendogli ripetere gli stessi esercizi studiati a memoria: punteggio alto, ma poca garanzia di riuscire con esercizi nuovi.

Per questo serve una misura out-of-sample: l’ELPD stimato con la validazione incrociata (LOO), che vedremo subito.

NotaEsempio pratico: calcolo della LPPD

Supponiamo un singolo dato \(y_i=3\) successi su \(n=5\). Abbiamo tre possibili valori di \(\theta\) dalla posteriore, con pesi didattici:

• \(\theta^{(1)}=0.3\), \(w^{(1)}=0.2\)
• \(\theta^{(2)}=0.5\), \(w^{(2)}=0.5\)
• \(\theta^{(3)}=0.7\), \(w^{(3)}=0.3\)

Per ciascun \(\theta^{(s)}\) calcoliamo \(p(y_i \mid \theta^{(s)})\), poi facciamo la media pesata (eq. Equazione 44.6), e infine il log.

y_i  <- 3
n_i  <- 5
theta_vals        <- c(0.3, 0.5, 0.7)
posterior_weights <- c(0.2, 0.5, 0.3)

# 1. Likelihood per ogni campione
likelihoods <- dbinom(y_i, size = n_i, prob = theta_vals)

# 2. Media pesata delle likelihood
p_yi_given_y <- sum(posterior_weights * likelihoods)

# 3. Log-score (per un solo dato)
log_score_i <- log(p_yi_given_y)

cbind(theta_vals, likelihoods)
#>      theta_vals likelihoods
#> [1,]        0.3       0.132
#> [2,]        0.5       0.312
#> [3,]        0.7       0.309
p_yi_given_y
#> [1] 0.275
log_score_i
#> [1] -1.29

Con più osservazioni \({y_i}_{i=1}^n\), la LPPD è la somma dei log-score punto per punto (Equazione 44.7).

44.4 Dalla LPPD all’ELPD: valutazione out-of-sample

44.4.1 Il problema dell’overfitting spiegato

Immagina di voler valutare la capacità di uno studente di riconoscere emozioni nei volti.

1. Scenario A: lo testi sempre con le stesse fotografie viste durante l’allenamento.
2. Scenario B: lo testi con fotografie nuove di persone mai viste prima.

Nel primo caso, lo studente probabilmente ottiene un punteggio alto, ma non sappiamo se ha davvero imparato a riconoscere le emozioni o se sta semplicemente ricordando quelle immagini specifiche. Il secondo scenario è più onesto: misura la capacità di generalizzazione.

Questo è esattamente il passaggio da LPPD a ELPD:

• LPPD (in-sample) ≈ Scenario A: valuta quanto bene il modello “ricorda” i dati visti.
• ELPD (out-of-sample) ≈ Scenario B: valuta quanto bene il modello generalizza a casi nuovi.

44.4.2 Guardare oltre i dati osservati

Quando valutiamo un modello, non basta chiedersi quanto bene spiega i dati già visti. La domanda giusta è: quanto bene predirebbe dati nuovi, mai osservati?

La Expected Log Predictive Density (ELPD) risponde proprio a questa domanda. La logica è la stessa della LPPD (si sommano log-verosimiglianze punto per punto), ma con una differenza cruciale: per valutare \(y_i\) non usiamo \(y_i\) nell’addestramento.

La tecnica operativa è la Leave-One-Out (LOO):

\[ \widehat{\mathrm{ELPD}}_{\mathrm{LOO}} \;=\; \sum_{i=1}^n \log p(y_i \mid y_{-i}), \tag{44.8}\] dove \(y_{-i}\) indica il dataset a cui è stata tolta l’osservazione \(i\). In questo modo ogni punto è valutato fuori campione, proprio come nello Scenario B.

44.4.3 Un esempio concreto per chiarire la differenza

Supponi di costruire un modello che predice i punteggi di memoria a breve termine a partire dal livello di concentrazione.

• Con LPPD, valuti il modello sugli stessi studenti usati per stimarlo: “Quanto bene spiega questi dati noti?”
• Con ELPD (LOO), togli uno studente alla volta, stimi il modello sugli altri e predici il punteggio dello studente escluso: “Quanto bene predirebbe un nuovo studente?”

44.4.3.1 Procedura operativa

1. Esclusione: si rimuove la prima osservazione \(y_1\) dal dataset.
2. Stima del modello: si adatta il modello ai dati rimanenti \(y_{-1}\).
3. Valutazione predittiva: si calcola la densità predittiva \(p(y_1 \mid y_{-1})\) e il corrispondente logaritmo.
4. Iterazione: si ripete la procedura per ogni osservazione \(i = 1, \dots, n\).
5. Aggregazione: si sommano tutti i valori di log-verosimiglianza ottenuti, ricavando così \(\widehat{\mathrm{ELPD}}_{\mathrm{LOO}}\).

ConsiglioPunti chiave operativi

• LPPD misura l’accuratezza sui dati di stima e può sovrastimare le prestazioni per modelli flessibili (rischio di overfitting).
• ELPD (LOO) fornisce una stima più realistica della capacità di generalizzazione del modello.
• Nella pratica, riadattare il modello \(n\) volte risulta computazionalmente oneroso: si utilizza PSIS-LOO per stimare \(\widehat{\mathrm{ELPD}}_{\mathrm{LOO}}\) a partire da un singolo adattamento, con diagnostiche (Pareto-k) per valutarne l’affidabilità.
• In presenza di dati dipendenti o gerarchici, l’unità di “leave-one-out” deve corrispondere al cluster (ad esempio, il soggetto), non alla singola osservazione.

44.5 L’ELPD e la divergenza di Kullback–Leibler

Abbiamo visto che l’ELPD fornisce una misura empirica della capacità predittiva di un modello. Ma perché proprio l’ELPD è la scelta corretta? La giustificazione profonda poggia su un concetto centrale della teoria dell’informazione: la divergenza di Kullback–Leibler (KL).

44.5.1 Cosa misura la divergenza KL?

La divergenza KL, \(D_{\text{KL}}\), quantifica la distanza informazionale tra:

• la distribuzione vera dei dati futuri \(p(\tilde{y})\), che rappresenta la realtà (sconosciuta);
• la distribuzione predittiva del modello \(q(\tilde{y}\mid y)\), che rappresenta la nostra approssimazione.

Formalmente:

\[ D_{\text{KL}}(p \parallel q) \;=\; \mathbb{E}_p \left[ \log \frac{p(\tilde{y})}{q(\tilde{y} \mid y)} \right], \] dove l’aspettativa \(\mathbb{E}_p\) è presa rispetto alla distribuzione vera \(p\).

44.5.2 Scomposizione della divergenza KL

Possiamo riscrivere la formula come:

\[ D_{\text{KL}}(p \parallel q) \;=\; \underbrace{\mathbb{E}_{p}\big[\log p(\tilde{y})\big]}_{-H(p)} \;-\; \underbrace{\mathbb{E}_{p}\big[\log q(\tilde{y} \mid y)\big]}_{\text{accuratezza predittiva (ELPD teorica)}}. \tag{44.9}\]

Nota: \(H(p) = -\mathbb{E}_{p}[\log p(\tilde{y})]\), quindi \(\mathbb{E}_{p}[\log p(\tilde{y})] = -H(p)\).

• (1) Entropia di \(p\): rappresenta l’incertezza intrinseca della distribuzione vera. È una costante, indipendente dal modello.
• (2) Accuratezza predittiva: misura quanto il modello \(q\) assegna probabilità alte agli eventi che davvero accadono. Questo termine è esattamente ciò che stimiamo con l’ELPD.

44.5.3 Il collegamento fondamentale

Poiché l’entropia di \(p\) è costante, minimizzare la divergenza KL equivale a massimizzare l’ELPD.

Date due distribuzioni predittive \(q_A\) e \(q_B\), la differenza nelle rispettive divergenze KL si riduce a:

\[ D_{\text{KL}}(p \parallel q_A) - D_{\text{KL}}(p \parallel q_B) = \mathrm{ELPD}(q_B) - \mathrm{ELPD}(q_A). \]

44.5.4 Conclusione teorica fondamentale

\[ \text{Massimizzare l’ELPD} \;\;\equiv\;\; \text{Minimizzare la divergenza KL dalla verità}. \]

In altre parole, scegliere il modello con l’ELPD più alto significa scegliere quello la cui distribuzione predittiva è, in media, più vicina alla realtà.

Questo principio unifica teoria dell’informazione e pratica statistica: ELPD non è solo un criterio empirico, ma la traduzione operativa del principio di minimizzazione della distanza dalla verità.

Nota

Messaggio chiave. L’ELPD è la stima empirica di (meno) la divergenza KL. Più alto è l’ELPD → più bassa è la divergenza KL → migliore è la capacità predittiva.

44.5.4.1 Esempio numerico: ELPD e KL

Supponiamo di voler predire il numero di “teste” in \(n=10\) lanci di una moneta:

• Distribuzione vera: \(p(y) = \text{Binom}(n=10, p=0.6)\).
• Modello candidato: \(q(y) = \text{Binom}(n=10, q=0.5)\).

Calcoliamo l’ELPD per il modello candidato e confrontiamolo con quello del “modello vero”:

n <- 10
p <- 0.6   # probabilità vera
q <- 0.5   # probabilità del modello candidato

y_vals <- 0:n
p_y    <- dbinom(y_vals, size = n, prob = p)      # distribuzione vera
log_q  <- log(dbinom(y_vals, size = n, prob = q)) # log-predittiva modello q
log_p  <- log(dbinom(y_vals, size = n, prob = p)) # log-predittiva modello p

# ELPD
elpd_q <- sum(p_y * log_q)
elpd_p <- sum(p_y * log_p)

# Divergenza KL
kl_pq  <- sum(p_y * (log_p - log_q))

cat(sprintf("ELPD modello candidato (q=0.5): %.4f\n", elpd_q))
#> ELPD modello candidato (q=0.5): -2.0549
cat(sprintf("ELPD modello vero      (q=0.6): %.4f\n", elpd_p))
#> ELPD modello vero      (q=0.6): -1.8536
cat(sprintf("Differenza ELPD (vero - candidato): %.4f\n", elpd_p - elpd_q))
#> Differenza ELPD (vero - candidato): 0.2014
cat(sprintf("KL(p || q): %.4f\n", kl_pq))
#> KL(p || q): 0.2014

Risultati attesi:

• L’ELPD del modello candidato è più basso rispetto a quello del modello vero.
• La differenza tra i due ELPD coincide con la divergenza KL \(D_{\text{KL}}(p|q)\).

Questo mostra sia algebricamente sia numericamente che massimizzare ELPD = minimizzare KL.

44.6 Dalla teoria alla pratica: LOO-CV e PSIS-LOO

Come abbiamo stabilito, l’ELPD rappresenta la misura ideale per valutare la capacità predittiva di un modello, grazie alla sua diretta connessione con la divergenza di Kullback-Leibler. Tuttavia, la sua definizione teorica richiede il calcolo di un’aspettativa rispetto alla vera distribuzione dei dati futuri \(p(\tilde{y})\), che nella pratica risulta sempre ignota.

In che modo è possibile stimare questa quantità? La soluzione metodologicamente più solida è offerta dalla Leave-One-Out Cross-Validation (LOO-CV), una procedura che consente di approssimare l’ELPD teorico utilizzando esclusivamente i dati osservati.

44.6.0.1 Idea e procedura della LOO-CV

L’approccio Leave-One-Out si basa su un concetto intuitivo ma potente: utilizzare ciascuna osservazione del dataset come dato “nuovo”, valutando la capacità del modello - stimato sulle restanti osservazioni - di prevederla adeguatamente.

La procedura si articola nei seguenti passi:

1. Isolamento: si seleziona un’osservazione \(y_i\) dal dataset.
2. Esclusione: la si rimuove temporaneamente dal campione di stima.
3. Stima: si adatta il modello ai dati rimanenti \(y_{-i}\).
4. Valutazione: si calcola la densità predittiva per l’osservazione esclusa, \(q(y_i \mid y_{-i})\).
5. Iterazione: si ripete il processo per tutte le \(n\) osservazioni.
6. Aggregazione: si sommano i logaritmi delle densità predittive ottenute.

Formalmente:

\[ \mathrm{ELPD}_{\mathrm{LOO}} \;=\; \sum_{i=1}^n \log q(y_i \mid y_{-i}) . \tag{44.10}\] Il risultato è una stima out-of-sample genuina: ciascuna osservazione viene valutata da un modello stimato in sua assenza, garantendo così una valutazione imparziale della capacità predittiva.

44.6.1 Perché funziona

La LOO-CV sostituisce l’aspettativa teorica rispetto a \(p(\tilde{y})\) con una media empirica sulle osservazioni reali. Ogni \(y_i\) è trattata come nuovo dato generato da \(p\), e la media dei log-score fuori campione fornisce una stima della capacità predittiva attesa:

\[ \mathrm{ELPD}_{\mathrm{LOO}} \;\approx\; \mathbb{E}_{p}\!\left[\log q(\tilde{y}\mid y)\right]. \tag{44.11}\]

44.6.2 Confrontare i modelli con LOO-CV

Una volta ottenuta la stima dell’ELPD-LOO per ciascun modello, il confronto avviene calcolando la differenza:

\[ \Delta \mathrm{ELPD} \;=\; \mathrm{ELPD}_{\mathrm{LOO}}(M_1) - \mathrm{ELPD}_{\mathrm{LOO}}(M_2). \tag{44.12}\]

Un valore positivo di \(\Delta \mathrm{ELPD}\) indica una migliore capacità predittiva del modello \(M_1\) rispetto a \(M_2\). Per una valutazione completa, è fondamentale accompagnare questa differenza con il suo errore standard: secondo una regola empirica consolidata, un vantaggio pari ad almeno due errori standard può considerarsi credibile.

Coerenza metodologica: nei confronti tra modelli, è essenziale mantenere la stessa definizione puntuale della verosimiglianza (identica unità di osservazione o cluster) attraverso tutti i modelli analizzati.

44.6.3 Overfitting e vantaggi della validazione LOO

La valutazione di un modello sugli stessi dati utilizzati per la stima tende a sovrastimarne sistematicamente le prestazioni, un fenomeno noto come overfitting. La validazione Leave-One-Out (LOO) supera questa limitazione strutturale: ciascuna osservazione \(y_i\) viene valutata da un modello stimato escludendo proprio \(y_i\), producendo così una stima più realistica della capacità predittiva del modello su nuovi dati.

44.6.3.1 L’approssimazione PSIS-LOO e le diagnostiche

Il refitting ripetuto del modello per ciascuna osservazione risulta spesso proibitivo dal punto di vista computazionale. Il metodo Pareto-Smoothed Importance Sampling LOO (PSIS-LOO) supera questa limitazione stimando \(\mathrm{ELPD}_{\mathrm{LOO}}\) a partire da un’unica campionatura MCMC, offrendo al contempo strumenti diagnostici per valutarne l’affidabilità.

• Diagnostica Pareto-\(k\) (interpretazione orientativa):
  • \(k \leq 0.5\): stima affidabile
  • \(0.5 < k \leq 0.7\): attendibilità moderata, interpretare con cautela
  • \(k > 0.7\): possibile instabilità → considerare alternative come K-fold CV, LOO per cluster o refitting mirato
• Vantaggi principali: efficienza computazionale, stima puntuale dell’ELPD con errore standard, e controllo qualità attraverso il parametro \(k\).

In R, le metodologie descritte sono implementate nel pacchetto loo, che risulta già integrato negli ambienti brms e rstanarm attraverso le funzioni loo() e loo_compare(). Per gli utenti che lavorano direttamente con Stan, è sufficiente calcolare la log-verosimiglianza punto per punto all’interno del blocco generated quantities e successivamente passare l’output alla funzione loo::loo().

NotaValidazione per cluster

In presenza di dati dipendenti (come nei studi EMA, misure ripetute o trial multipli per soggetto), la procedura LOO deve essere applicata a livello di cluster. Ciò significa che per ogni validazione si esclude l’intero gruppo di osservazioni correlate (ad esempio, tutti i dati di un soggetto o di una sessione sperimentale), preservando così la struttura di dipendenza dei dati.

44.6.3.2 Sintesi dei criteri predittivi

• L’ELPD quantifica la capacità predittiva del modello su dati non osservati.
• In assenza della distribuzione vera \(p(\tilde{y})\), la procedura LOO-CV ne fornisce una stima empirica.
• Le differenze nell’ELPD-LOO tra modelli approssimano le relative divergenze di Kullback-Leibler.
• L’implementazione PSIS-LOO garantisce efficienza computazionale e include diagnostiche di affidabilità (parametro \(k\) di Pareto).
• Linea guida operativa: privilegiare i modelli con ELPD-LOO più elevato, riportando sistematicamente \(\Delta \mathrm{ELPD}\) con il corrispondente errore standard, e integrando la valutazione con considerazioni sulla parsimonia e l’interpretabilità.

NotaEsempio: confronto ELPD-LOO tra due modelli

Questo esempio mostra come passare dalla definizione teorica dell’ELPD alla stima pratica via LOO, nel caso elementare Beta–Bernoulliano.

Dati: \(y=\{1,1,1,0,1\}\) (cinque prove indipendenti: quattro successi, un insuccesso).

Modello A (Bayesiano): \(Y\sim\mathrm{Bernoulli}(\theta)\), prior \(\theta\sim \mathrm{Beta}(1,1)\). Per LOO, per ogni \(i\):

• escludi \(y_i\);
• posteriore \(\theta \mid y_{-i} \sim \mathrm{Beta}(1+s_{-i},\,1+(n-1)-s_{-i})\);
• predittiva per \(y_i\): \(q(y_i=1\mid y_{-i}) = \dfrac{\alpha}{\alpha+\beta}\).

Modello B (confronto): moneta equa fissa (\(q=0.5\)).

# Dati
y <- c(1, 1, 1, 0, 1)
n <- length(y)

# Log-predittiva LOO per Modello A (Beta(1,1) + Bernoulli)
loo_log_pred_beta <- function(i, y, a0 = 1, b0 = 1) {
  yi <- y[i]
  s_minus <- sum(y) - yi         # successi tra i rimanenti
  n_minus <- n - 1
  alpha <- a0 + s_minus
  beta  <- b0 + (n_minus - s_minus)
  p1 <- alpha / (alpha + beta)   # P(y=1 | y_-i)
  p  <- if (yi == 1) p1 else (1 - p1)
  log(p)
}

# Log-predittive punto-per-punto
lp_A <- sapply(seq_along(y), loo_log_pred_beta, y = y)
lp_B <- rep(log(0.5), n)

# ELPD-LOO
elpd_A <- sum(lp_A)
elpd_B <- sum(lp_B)

# Differenza e SE (stima semplice da varianza dei contributi puntowise)
diff_pt <- lp_A - lp_B
se_diff <- sqrt(n * var(diff_pt))

# Tabella riassuntiva
res <- data.frame(
  i = 1:n, y = y,
  lp_A = round(lp_A, 6),
  lp_B = round(lp_B, 6),
  diff = round(diff_pt, 6)
)
print(res)
#>   i y   lp_A   lp_B   diff
#> 1 1 1 -0.405 -0.693  0.288
#> 2 2 1 -0.405 -0.693  0.288
#> 3 3 1 -0.405 -0.693  0.288
#> 4 4 0 -1.792 -0.693 -1.099
#> 5 5 1 -0.405 -0.693  0.288
cat(sprintf("\nELPD-LOO Modello A: %.6f\n", elpd_A))
#> 
#> ELPD-LOO Modello A: -3.413620
cat(sprintf("ELPD-LOO Modello B: %.6f\n", elpd_B))
#> ELPD-LOO Modello B: -3.465736
cat(sprintf("Differenza (A-B)  : %.6f\n", elpd_A - elpd_B))
#> Differenza (A-B)  : 0.052116
cat(sprintf("SE differenza     : %.6f\n", se_diff))
#> SE differenza     : 1.386294

Interpretazione: ogni riga è il log-score fuori campione dei due modelli. Con 4/5 successi, il Modello A assegna probabilità \(>0.5\) ai successi (e \(<0.5\) all’unico insuccesso), quindi tende ad avere ELPD-LOO più alto. Con \(n\) piccolo, l’SE è grande: la lezione è didattica, non inferenziale.

Regola pratica: \(|\Delta \mathrm{ELPD}| \ge 2 \times \mathrm{SE}\) → evidenza credibile a favore del modello migliore.

NotaImplementazione pratica: stima e confronto dell’ELPD-LOO

Strumenti computazionali

• Funzioni loo::loo() e loo::loo_compare() nel pacchetto R omonimo; già integrate negli ambienti brms e rstanarm.
• Implementazione Stan personalizzata: memorizzare la log-verosimiglianza puntuale nel blocco generated quantities e fornire la matrice risultante a loo().

Procedura operativa

1. Stimare ciascun modello candidato (brm(), stan_glm() o codice Stan personalizzato).
2. Estrarre la log-verosimiglianza puntuale (tramite log_lik() o variabile definita in Stan).
3. Applicare loo() per ottenere ELPD, errore standard e diagnostiche Pareto-\(k\)
4. Confrontare i modelli con loo_compare(), riportando \(\Delta \mathrm{ELPD}\) ed errore standard.
5. In presenza di valori elevati di Pareto-\(k\): ricorrere a K-fold CV, LOO per cluster, o valutare un refitting più robusto.

Criteri decisionali

• Privilegiare i modelli con valori più elevati di ELPD-LOO.
• Comunicare sistematicamente \(\Delta \mathrm{ELPD}\) accompagnato dall’errore standard.
• Per differenze inferiori a \(2\times \mathrm{SE}\), considerare prioritariamente parsimonia e interpretabilità.
• Nel caso di serie temporali, adottare schemi di validazione leave-future-out strutturati a blocchi temporali.

44.7 Criteri di informazione e confronto tra modelli

La metrica fondamentale per la valutazione predittiva è l’ELPD (Expected Log Predictive Density), che misura quanto la distribuzione predittiva del modello \(q(\tilde{y} \mid y)\) si avvicini alla distribuzione vera dei dati \(p(\tilde{y})\), minimizzando la divergenza di Kullback-Leibler. I criteri di informazione costituiscono diverse approssimazioni - di varia accuratezza - per stimare questa capacità predittiva su nuovi dati, bilanciando in modo ottimale il trade-off tra adattamento ai dati osservati e complessità del modello.

Tra i vari criteri disponibili, LOO e WAIC rappresentano le approssimazioni concettualmente più vicine all’ELPD in un’ottica bayesiana, mentre AIC e BIC si collocano nell’ambito della statistica frequentista, dove funzionano bene per modelli con strutture parametriche semplici e dati indipendenti. Essendo statistiche asintotiche, il loro utilizzo richiede campioni di dati di numerosità elevata per garantire stime affidabili.

44.7.0.1 Tipi di criteri e formule principali

AIC – Per modelli classici
Stima la divergenza KL basandosi sulla verosimiglianza massimizzata.
\(\text{AIC} = -2\log p(y\mid \hat\theta) + 2k\)
Penalità fissa: \(2k\) parametri
Ideale per: campioni grandi, modelli regolari senza gerarchie

BIC – Per modelli parsimoniosi
Approssima la verosimiglianza marginale, favorendo la semplicità.
\(\text{BIC} = -2\log p(y\mid \hat\theta) + k\log n\)
Penalità crescente: \(k\log n\)
Adatto a: selezione del modello “vero”, non predizione

WAIC – Approccio bayesiano
Stima direttamente l’accuratezza predittiva in ottica bayesiana.
\(\text{elpd}_{\text{waic}} = \sum_i \log \big(\tfrac{1}{S}\sum_s p(y_i\mid \theta^{(s)})\big) - p_{\text{waic}}\)
Penalità adattiva: \(p_{\text{waic}} = \sum_i V_s[\log p(y_i\mid \theta^{(s)})]\)
Utilizzo: alternativa pratica a LOO, richiede likelihood puntuale

LOO (PSIS-LOO) – Riferimento predittivo
Stima l’ELPD mediante validazione leave-one-out.
\(\text{elpd}_{\text{loo}} = \sum_i \log q(y_i\mid y_{-i})\)
Penalità stimata: \(p_{\text{loo}} \approx \sum_i (\text{lppd}_i - \text{elpd}_{\text{loo},i})\)
Vantaggi: efficiente con PSIS, include diagnostiche di affidabilità

\(V_s[\cdot] = \operatorname{Var}_s[\cdot]\); \(\text{lppd}_i = \log\big(\tfrac{1}{S}\sum_s p(y_i\mid \theta^{(s)})\big)\)

ConsiglioScala e convenzioni (Stan/loo)

• Scala ELPD: più alto è meglio.
• Scala IC (devianza): \(\text{LOOIC}=-2,\text{elpd}*{\text{loo}}\), \(\text{WAIC}=-2,\text{elpd}*{\text{waic}}\) → più basso è meglio.
• In loo, brms e rstanarm troverai sia \(\text{elpd}\) e SE, sia l’IC corrispondente. Riporta sempre \(\Delta \text{ELPD}\) con il suo SE.

44.7.0.2 Quando usare AIC, BIC, WAIC o LOO

• Predizione bayesiana su nuovi dati → PSIS-LOO costituisce lo strumento privilegiato, con WAIC come utile riscontro, data la frequente convergenza dei risultati.

• Modelli frequentisti regolari → AIC offre una soluzione immediata per campioni di ampiezza medio-elevata e strutture semplici, mentre BIC risulta preferibile quando si intende enfatizzare la parsimonia o la consistenza asintotica.

• Strutture gerarchiche e dati dipendenti → è fondamentale identificare correttamente l’unità di predizione pointwise, generalmente il cluster; in tali contesti, LOO e WAIC mantengono la loro validità, a differenza di AIC e BIC.

• Analisi di serie storiche → si raccomandano protocolli di validazione leave-future-out strutturati a blocchi, in alternativa all’applicazione standard del LOO.

44.7.0.3 Penalizzazione effettiva e complessità reale

Nei modelli bayesiani, la complessità non corrisponde al mero numero di parametri, ma dipende criticamente da quanto i dati siano informativi riguardo ai parametri stessi. I criteri WAIC e LOO quantificano una penalizzazione effettiva (\(p_{\text{waic}}\), \(p_{\text{loo}}\)) basata sulla variabilità della log-verosimiglianza posteriore punto per punto. Questo approccio genera una penalizzazione adattiva, che riflette più fedelmente la complessità reale del modello, particolarmente in contesti gerarchici o non regolari.

44.8 Workflow operativo e buone pratiche

44.8.1 Procedura passo-passo

Preparazione: definisci l’unità di validazione appropriata (singole osservazioni per dati indipendenti, interi cluster per dati raggruppati).

1. Stima i modelli candidati
2. Estrai la verosimiglianza punto per punto dalle simulazioni posteriori
3. Calcola PSIS-LOO (loo()) per ottenere:
   • \(\mathrm{elpd}_{\mathrm{loo}}\) e errore standard
   • Diagnostiche Pareto-\(k\)
4. Confronta i modelli con loo_compare()
5. Interpreta i risultati considerando:
   • Differenze di ELPD e relativa incertezza
   • Semplicità e interpretabilità del modello
   • Adeguatezza al contesto applicativo

44.8.1.1 Errori comuni e raccomandazioni

• Metriche in-sample → Evita MSE sui dati di training (sovrastima le prestazioni). Usa ELPD-LOO o WAIC.
• Scale incoerenti → Non confrontare direttamente AIC/BIC con ELPD. Calcola differenze sulla stessa scala (ΔELPD).
• Validazione impropria → Con dati raggruppati, escludi interi cluster, non singole osservazioni.
• Diagnostiche ignorate → Con PSIS-LOO:
  • \(k \le 0.5\): procedi con fiducia
  • \(0.5 < k \le 0.7\): risultati da interpretare con cautela
  • \(k > 0.7\): passa a K-fold CV per stime più robuste
• Incertezza trascurata → Riporta sempre l’errore standard di ΔELPD.
• Pre-processing scorretto → Trasformazioni e selezione di variabili vanno ripetute per ogni fold di validazione.

44.8.1.2 Verifica incrociata con WAIC

(Raccomandato) Calcola WAIC come verifica aggiuntiva di robustezza. I due criteri (LOO e WAIC) dovrebbero convergere verso conclusioni simili.

Orientamento finale: Tutti i criteri mirano a stimare l’ELPD. AIC e BIC offrono approssimazioni classiche, mentre WAIC e LOO rappresentano le controparti bayesiane più dirette. La scelta dello strumento dipende dalla struttura dei dati, ma la comunicazione dei risultati deve sempre includere sia l’entità delle differenze che la loro incertezza.

Riflessioni conclusive

Il principio guida di questo capitolo è semplice: un buon modello è quello che fa buone previsioni su dati che non ha mai visto. Nel paradigma bayesiano, questo si traduce nell’uso della distribuzione predittiva a posteriori \(p(\tilde y \mid y)\) e nella sua bontà, misurata dall’Expected Log Predictive Density (ELPD)—l’accuratezza predittiva attesa su dati futuri. Dal punto di vista teorico, massimizzare l’ELPD equivale a minimizzare la divergenza di Kullback-Leibler dalla verità, creando un ponte diretto tra teoria dell’informazione e pratica della modellistica.

Sul piano operativo, l’ELPD viene stimato in modo efficiente con PSIS-LOO, che, partendo da un unico adattamento del modello, fornisce anche le fondamentali diagnostiche Pareto-\(k\) per valutarne l’affidabilità. WAIC rappresenta un’alternativa pienamente bayesiana, spesso in accordo con LOO. AIC e BIC, sebbene utili in contesti classici e regolari, possono risultare fuorvianti in scenari bayesiani con piccoli campioni, strutture gerarchiche o dipendenze complesse, poiché si basano su stime puntuali e approssimazioni asintotiche.

Nel confronto tra modelli, non ci si deve limitare a dichiarare un vincitore. È una buona pratica di reporting includere:

• I valori di \(\text{ELPD}_{\text{LOO}}\) (o LOOIC, su scala devianza);
• Le differenze \(\Delta \mathrm{ELPD}\) con i relativi errori standard;
• Un riepilogo delle diagnostiche Pareto-\(k\) (es., numero di osservazioni con \(k>0.7\)) e le eventuali azioni intraprese (K-fold, LOO per cluster, refit mirato);
• L’ELPD per osservazione (\(\mathrm{ELPD}/n\)) quando si confrontano dataset di dimensioni diverse.

È cruciale scegliere l’unità di predizione corretta (ad esempio, per cluster in dati gerarchici) e adottare schemi leave-future-out in presenza di dipendenze temporali. Ricorda che molti software (come l’ecosistema Stan/loo) riportano i risultati anche in scala devianza (\(\text{LOOIC} = -2 \cdot \text{ELPD}_{\text{LOO}}\)): in questo caso, valori più bassi indicano un modello migliore.

Infine, la selezione del modello non è un processo meccanico. Differenze esigue nell’ELPD, specialmente se accompagnate da un errore standard ampio, raramente hanno rilevanza sostanziale. In questi casi, è preferibile privilegiare parsimonia e coerenza teorica. L’obiettivo ultimo non è solo prevedere, ma comprendere. L’ELPD e la LOO forniscono una bussola quantitativa preziosa, che deve essere integrata con la conoscenza del dominio, la struttura dei dati e l’interpretabilità del modello.

Messaggio fondamentale: il modello migliore è quello che predice meglio il futuro senza sacrificare interpretabilità e coerenza teorica.

NotaInformazioni sull’ambiente di sviluppo

sessionInfo()
#> R version 4.5.1 (2025-06-13)
#> Platform: aarch64-apple-darwin20
#> Running under: macOS Tahoe 26.0
#> 
#> Matrix products: default
#> BLAS:   /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib 
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
#> 
#> locale:
#> [1] C.UTF-8/UTF-8/C.UTF-8/C/C.UTF-8/C.UTF-8
#> 
#> time zone: Europe/Rome
#> tzcode source: internal
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#>  [1] gt_1.1.0              ragg_1.5.0            tinytable_0.13.0     
#>  [4] withr_3.0.2           systemfonts_1.2.3     patchwork_1.3.2      
#>  [7] ggdist_3.3.3          tidybayes_3.0.7       bayesplot_1.14.0     
#> [10] ggplot2_4.0.0         reliabilitydiag_0.2.1 priorsense_1.1.1     
#> [13] posterior_1.6.1       loo_2.8.0             rstan_2.32.7         
#> [16] StanHeaders_2.32.10   brms_2.23.0           Rcpp_1.1.0           
#> [19] sessioninfo_1.2.3     conflicted_1.2.0      janitor_2.2.1        
#> [22] matrixStats_1.5.0     modelr_0.1.11         tibble_3.3.0         
#> [25] dplyr_1.1.4           tidyr_1.3.1           rio_1.2.4            
#> [28] here_1.0.2           
#> 
#> loaded via a namespace (and not attached):
#>  [1] gridExtra_2.3         inline_0.3.21         sandwich_3.1-1       
#>  [4] rlang_1.1.6           magrittr_2.0.4        multcomp_1.4-28      
#>  [7] snakecase_0.11.1      compiler_4.5.1        vctrs_0.6.5          
#> [10] stringr_1.5.2         pkgconfig_2.0.3       arrayhelpers_1.1-0   
#> [13] fastmap_1.2.0         backports_1.5.0       rmarkdown_2.30       
#> [16] purrr_1.1.0           xfun_0.53             cachem_1.1.0         
#> [19] jsonlite_2.0.0        broom_1.0.10          parallel_4.5.1       
#> [22] R6_2.6.1              stringi_1.8.7         RColorBrewer_1.1-3   
#> [25] lubridate_1.9.4       estimability_1.5.1    knitr_1.50           
#> [28] zoo_1.8-14            Matrix_1.7-4          splines_4.5.1        
#> [31] timechange_0.3.0      tidyselect_1.2.1      abind_1.4-8          
#> [34] codetools_0.2-20      curl_7.0.0            pkgbuild_1.4.8       
#> [37] lattice_0.22-7        bridgesampling_1.1-2  S7_0.2.0             
#> [40] coda_0.19-4.1         evaluate_1.0.5        survival_3.8-3       
#> [43] RcppParallel_5.1.11-1 xml2_1.4.0            pillar_1.11.1        
#> [46] tensorA_0.36.2.1      checkmate_2.3.3       stats4_4.5.1         
#> [49] distributional_0.5.0  generics_0.1.4        rprojroot_2.1.1      
#> [52] rstantools_2.5.0      scales_1.4.0          xtable_1.8-4         
#> [55] glue_1.8.0            emmeans_1.11.2-8      tools_4.5.1          
#> [58] fs_1.6.6              mvtnorm_1.3-3         grid_4.5.1           
#> [61] QuickJSR_1.8.1        colorspace_2.1-2      nlme_3.1-168         
#> [64] cli_3.6.5             textshaping_1.0.3     svUnit_1.0.8         
#> [67] Brobdingnag_1.2-9     V8_8.0.0              gtable_0.3.6         
#> [70] digest_0.6.37         TH.data_1.1-4         htmlwidgets_1.6.4    
#> [73] farver_2.1.2          memoise_2.0.1         htmltools_0.5.8.1    
#> [76] lifecycle_1.0.4       MASS_7.3-65

Bibliografia

McElreath, R. (2020). Statistical rethinking: A Bayesian course with examples in R and Stan (2nd Edition). CRC Press.

---

1. L’errore quadratico medio (MSE) misura la discrepanza tra valori previsti e osservati: \(\text{MSE} = \frac{1}{n}\sum_{i=1}^n (y_i - \hat{y}_i)^2\). Sebbene utile, non cattura l’incertezza probabilistica delle previsioni.