46  Valutare i modelli bayesiani

Introduzione

Nei capitoli precedenti abbiamo visto due concetti fondamentali: l’entropia, che misura l’incertezza insita in una distribuzione, e la divergenza di Kullback–Leibler (\(D_{\text{KL}}\)), che quantifica la distanza tra due distribuzioni di probabilità. Ora possiamo fare un passo ulteriore: usare queste idee per valutare e confrontare modelli statistici nel contesto bayesiano.

Il punto di partenza è una domanda cruciale: quanto bene il modello riesce a prevedere nuovi dati? Un buon modello non deve solo adattarsi ai dati osservati, ma deve anche saper generalizzare a situazioni future. Questa distinzione – adattamento vs. generalizzazione – è il cuore della valutazione predittiva.

Immaginiamo, ad esempio, di sviluppare un test psicologico per stimare l’ansia degli studenti prima di un esame. Non basta sapere che il modello descrive bene il campione usato per costruirlo: vogliamo essere ragionevolmente sicuri che le stesse previsioni funzionino anche per studenti che non hanno partecipato allo studio. In psicologia, scegliere tra due modelli è simile a decidere quale test usare: in entrambi i casi si cerca lo strumento che fornisce previsioni più affidabili.

Panoramica del capitolo

Per rispondere alla domanda fondamentale sulla qualità predittiva, seguiremo un percorso logico che ci condurrà dagli strumenti teorici ai metodi pratici:

  • Prima costruiremo la base teorica: vedremo come la distribuzione predittiva posteriore incorpora l’incertezza sui parametri nelle nostre previsioni
  • Poi definiremo le misure di accuratezza: il log-score per valutare la bontà delle previsioni punto per punto
  • Distingueremo tra valutazione in-sample e out-of-sample: LPPD (sui dati osservati) vs. ELPD (capacità di generalizzazione)
  • Collegheremo tutto alla divergenza KL: per capire perché massimizzare l’ELPD equivale a trovare il modello più vicino alla realtà
  • Implementeremo metodi pratici: LOO-CV per stimare l’ELPD senza conoscere la vera distribuzione dei dati
  • Confronteremo con altri criteri: AIC, BIC, WAIC e i loro ambiti di applicazione

L’obiettivo è fornire strumenti pratici e un quadro concettuale chiaro per guidare la scelta del modello più adatto al problema in esame.

  • Per comprendere appieno questo capitolo è utile leggere il capitolo 7 Ulysses’ Compass di Statistical Rethinking (McElreath (2020)).
here::here("code", "_common.R") |> 
  source()

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

46.1 Il punto di partenza: dalle previsioni deterministiche a quelle probabilistiche

Prima di addentrarci nei dettagli tecnici, facciamo un passo indietro per capire 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 “veri” e fissi. Se \(\hat{\theta}\) è la nostra stima, la previsione per un nuovo dato \(\tilde{y}\) è semplicemente:

\[ p(\tilde{y} \mid \hat{\theta}). \] Questo approccio ignora completamente l’incertezza sulla stima dei parametri.

Nell’approccio bayesiano, invece, riconosciamo che i parametri sono incerti. Anche dopo aver osservato i dati, la nostra conoscenza di \(\theta\) è descritta da un’intera distribuzione posteriore \(p(\theta \mid y)\), non da un singolo valore. Le previsioni devono quindi riflettere questa incertezza.

NotaAnalogia didattica

Immaginiamo di dover prevedere il tempo domani. L’approccio “classico” è come consultare un solo meteorologo e fidarsi completamente della sua previsione. L’approccio bayesiano è come consultare tutti i meteorologi disponibili, pesare le loro opinioni in base alla loro affidabilità passata, e costruire una previsione che tenga conto di tutti i punti di vista plausibili.

46.2 Distribuzione predittiva posteriore: il cuore delle previsioni bayesiane

Nel capitolo sul modello beta–binomiale l’abbiamo già incontrata: è 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 quantifica 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.

46.2.1 La formula fondamentale

Per prevedere un nuovo dato \(\tilde{y}\), non fissiamo \(\theta\). Mediamo invece tutte le previsioni condizionate \(p(\tilde{y} \mid \theta)\) pesandole con la densità posteriore \(p(\theta\mid y)\):

\[ p(\tilde{y} \mid y) = \int p(\tilde{y} \mid \theta)\, p(\theta \mid y)\, d\theta \tag{46.1}\] Questa è la distribuzione predittiva posteriore, e rappresenta la nostra migliore previsione per dati futuri dato quello che abbiamo osservato.

46.2.2 Decomposizione dell’integrale: cosa significa realmente

L’Equazione 46.1 può sembrare astratta, ma ha un’interpretazione intuitiva molto chiara:

  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)\): questa è la nostra incertezza su quale sia il vero valore di \(\theta\);
  3. L’integrale: combina le previsioni per tutti i possibili valori di \(\theta\), pesandole secondo quanto crediamo che ciascun valore sia plausibile.

Se conoscessimo il valore vero di \(\theta\), potremmo prevedere i dati futuri usando la distribuzione predittiva condizionata:

\[ p(\tilde y \mid \theta). \] Il problema è che \(\theta\) non lo conosciamo: abbiamo soltanto la distribuzione a posteriori \(p(\theta\mid y)\). Perciò, la distribuzione predittiva posteriore si costruisce combinando le previsioni condizionate per ogni valore possibile di \(\theta\), pesandole con quanto ciascun valore è plausibile a posteriori:

\[ p(\tilde y\mid y) = \int p(\tilde y\mid \theta)\,p(\theta\mid y)\,d\theta. \] Per fare un esempio concreto, consideriamo il caso binomiale. Supponiamo che i dati futuri siano 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}. \] La distribuzione predittiva posteriore diventa:

\[ p(\tilde y = x \mid y) = \int \binom{m}{x}\,\theta^x(1-\theta)^{m-x}\,p(\theta\mid y)\,d\theta. \] In alcuni casi particolari (per esempio con prior Beta e dati binomiali) questo integrale si può risolvere analiticamente, ottenendo la Beta–Binomiale. Ma in generale non c’è una formula chiusa e serve un’approssimazione numerica.

Approssimazione numerica con il metodo su griglia: L’idea è semplice: sostituire l’integrale con una somma pesata su una griglia di valori possibili di \(\theta\). I passaggi algoritmici sono i seguenti.

  1. Definire una griglia di valori di \(\theta\), ad esempio 1000 punti equispaziati tra 0 e 1:

    \[ \theta_1, \theta_2, \dots, \theta_J. \]

  2. Calcolare la posteriore su ciascun punto della griglia. Nel caso Beta–Binomiale:

    \[ p(\theta_j \mid y) \propto \theta_j^{\,k+a-1}(1-\theta_j)^{n-k+b-1}. \]

    Poi normalizzare per avere somme che valgono 1:

    \[ w_j = \frac{p(\theta_j \mid y)}{\sum_{\ell=1}^J p(\theta_\ell \mid y)}. \]

  3. Combinare le predizioni condizionate. Per ogni valore futuro \(x=0,\dots,m\), si calcola:

    \[ p(\tilde y = x \mid y) \approx \sum_{j=1}^J w_j \, \binom{m}{x}\theta_j^x(1-\theta_j)^{m-x}. \]

  4. Interpretazione: la pmf ottenuta è la nostra approssimazione numerica della distribuzione predittiva posteriore. Da essa possiamo:

    • calcolare probabilità,
    • generare campioni di \(\tilde y\),
    • confrontare la predizione con i dati osservati.

Da ricordare:

  • La predittiva non si ottiene facendo la media dei valori di \(\tilde y\), ma costruendo un’intera distribuzione di probabilità.
  • Il metodo su griglia è il più semplice: discretizza \(\theta\), pesa ogni valore con la sua plausibilità a posteriori, e combina le predizioni condizionate.
  • In problemi più complessi, la stessa logica viene implementata tramite MCMC: invece di usare una griglia fissa, si usano campioni \(\theta^{(s)}\) dalla posteriore.

Esaminiamo ora uno script in R che implementa passo per passo l’approssimazione della distribuzione predittiva posteriore binomiale con il metodo su griglia.

# ESEMPIO DIDATTICO: predittiva posteriore per Binomiale con metodo su griglia
# Dati e prior
k <- 10     # successi osservati
n <- 50     # prove osservate
a <- 1      # prior Beta(a, b)
b <- 1
m <- 10     # numero di prove future per la predizione (scelta didattica)
J <- 2000   # numero di punti griglia su theta in [0,1]

# -------------------------------------------------------------
# PASSAGGIO 1: Griglia su theta
# -------------------------------------------------------------

theta <- seq(0, 1, length.out = J)

# -------------------------------------------------------------
# PASSAGGIO 2: 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)

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

# -------------------------------------------------------------
# PASSAGGIO 3: combinare le predittive condizionate p(tilde y | theta)
# -------------------------------------------------------------
# Obiettivo: costruire la pmf predittiva p(tilde y = x | y) 
# per ogni x = 0,...,m come media pesata (sulla griglia di θ) 
# delle pmf condizionate binomiali.

# 1) Definiamo i valori futuri possibili di tilde y
x_vals <- 0:m

# 2) Inizializziamo una matrice vuota: 
#    - J righe (una per ciascun θ_j della griglia)
#    - (m+1) colonne (una per ogni valore possibile di x)
px_given_theta <- matrix(NA_real_, nrow = J, ncol = m + 1)

# 3) Riempiamo la matrice: per ogni θ_j (riga j) e per ogni x (colonna i)
#    calcoliamo P(tilde y = x | θ_j) = Binomiale(x | m, θ_j)
for (j in 1:J) {
  for (i in 1:(m + 1)) {
    x <- x_vals[i]
    px_given_theta[j, i] <- dbinom(x, size = m, prob = theta[j])
  }
}

# 4) Combinazione pesata:
#    p(tilde y = x | y) ≈ somma_j w_j * P(tilde y = x | θ_j).
#    Per ciascun valore di x (colonna i), facciamo la somma pesata.
pred_pmf <- numeric(m + 1)
for (i in 1:(m + 1)) {
  pred_pmf[i] <- sum(w * px_given_theta[, i])
}

# Nota didattica:
# - Ogni colonna della matrice px_given_theta contiene le probabilità condizionate 
#   P(tilde y = x_i | θ_j) per tutti i valori di griglia θ_j.
# - Moltiplicando riga per riga queste probabilità per i pesi posteriori w_j 
#   e sommando, otteniamo la probabilità predittiva p(tilde y = x_i | y).
# - In questo modo l'integrale viene approssimato da una somma pesata.

# -------------------------------------------------------------
# PASSAGGIO 4: Risultato: una pmf su {0,1,...,m}
# -------------------------------------------------------------

pred_df <- data.frame(x = x_vals, p = pred_pmf)
print(pred_df)
#>     x          p
#> 1   0 0.11391218
#> 2   1 0.25060680
#> 3   2 0.27617892
#> 4   3 0.19946255
#> 5   4 0.10397516
#> 6   5 0.04068593
#> 7   6 0.01205509
#> 8   7 0.00266151
#> 9   8 0.00041780
#> 10  9 0.00004200
#> 11 10 0.00000205
sum(pred_df$p)  # dovrebbe essere ~1
#> [1] 1
# (Opzionale) campionamento dalla predittiva posteriore approssimata
# Estrae N valori da {0,...,m} con le probabilità 'p'
set.seed(123)
N <- 5000
tilde_y_samples <- sample(pred_df$x, size = N, replace = TRUE, prob = pred_df$p)

# Controllo: istogramma delle simulazioni vs pmf teorica approssimata
ggplot() +
  geom_histogram(
    data = data.frame(x = tilde_y_samples), aes(x = x, y = after_stat(density)),
    binwidth = 1, breaks = seq(-0.5, m + 0.5, by = 1), fill = "skyblue", 
    color = "black") +
  geom_point(data = pred_df, aes(x = x, y = p), pch = 19, cex = 3) + 
  geom_line(data = pred_df, aes(x = x, y = p), lwd = 1.5) + 
  ylim(0, max(pred_df$p) * 1.1) +
  labs(
    title = "Posterior Predictive (grid) – m=10",
    x = expression(tilde(y)),
    y = "Density")

Nota: il vettore pred_df$p è la pmf della predittiva posteriore approssimata; da qui si leggono probabilità, si calcolano quantità riassuntive e si può estrarre \(\tilde y\).

Verifica quando esiste la formula chiusa: Quando prior e likelihood sono coniugate (Beta + Binomiale), la predittiva è Beta–Binomiale. Possiamo usarla solo come verifica didattica:

# Confronto con Beta-Binomiale (se applicabile)
a_post <- a + k
b_post <- b + n - k

# pmf beta-binomial (con funzione base: dbetabinom in VGAM, altrimenti la implementiamo)
dbetabinom <- function(x, m, a, b) {
  # Beta-Binomiale: choose(m, x) * Beta(x+a, m-x+b) / Beta(a, b)
  choose(m, x) * beta(x + a, m - x + b) / beta(a, b)
}

bb_pmf <- sapply(0:m, function(x) dbetabinom(x, m, a_post, b_post))
cbind(grid = pred_df$p, beta_binom = bb_pmf)[1:6, ]  # prime 6 righe a confronto
#>        grid beta_binom
#> [1,] 0.1139     0.1139
#> [2,] 0.2506     0.2506
#> [3,] 0.2762     0.2762
#> [4,] 0.1995     0.1995
#> [5,] 0.1040     0.1040
#> [6,] 0.0407     0.0407
max(abs(pred_df$p - bb_pmf))   # lo scarto massimo (dovrebbe essere ~0 con J grande)
#> [1] 1.68e-15

46.2.3 Notazione e terminologia

Notazione: Useremo talvolta la forma compatta \(p(\cdot \mid y)\) per indicare la predittiva posteriore del modello. Quando ci servirà evidenziare la previsione marginale per una singola osservazione \(y_i\), scriveremo:

\[ p(y_i \mid y) = \int p(y_i \mid \theta)\, p(\theta \mid y)\, d\theta, \] cioè la verosimiglianza \(p(y_i\mid\theta)\) integrata rispetto alla posteriore \(p(\theta\mid y)\).

Idea chiave: La predittiva posteriore propaga l’incertezza sui parametri alle previsioni. È questo passaggio a rendere le valutazioni predittive coerenti con il principio bayesiano, e quindi utilizzabili nel confronto tra modelli e nella stima di quantità legate alla “distanza” dal generatore dei dati.

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 bayesiana media Log-score, LPPD
\(p(y_i \mid y_{-i})\) Predizione LOO (leave-one-out) ELPD
\(p(\tilde{y} \mid y)\) Distribuzione predittiva complessiva Divergenza KL, confronto modelli

46.3 Il problema fondamentale della valutazione predittiva

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

46.3.1 Il dilemma teorico

Quando costruiamo un modello, vogliamo sapere quanto bene riesce a predire dati futuri. In altre parole: se ripetessimo l’esperimento o raccogliessimo nuovi dati, quanto sarebbero vicine le predizioni del modello a quei dati?

Per formalizzare questa idea, distinguiamo tra due distribuzioni:

  • la distribuzione vera dei dati futuri \(p(\tilde{y})\), che purtroppo non conosciamo,
  • la distribuzione predittiva del modello \(p(\tilde{y} \mid y)\), cioè le predizioni basate sui dati osservati \(y\).

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

46.3.2 Una misura di distanza: la divergenza di Kullback–Leibler

La divergenza di Kullback–Leibler (KL) fornisce una misura di questa distanza:

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

  • se \(q\) (il nostro modello) assegna probabilità alte agli stessi eventi che sono probabili in \(p\) (la realtà), la distanza sarà piccola,
  • se invece \(q\) “sbaglia bersaglio”, assegnando probabilità alte a eventi che in realtà sono rari, la distanza sarà grande.

Idea chiave. La KL divergence misura quanta informazione perdiamo se usiamo le predizioni del modello \(q\) al posto della distribuzione vera \(p\).

46.3.3 Un ostacolo pratico insormontabile

Il problema è che \(p(\tilde{y})\) non lo conosciamo mai: non abbiamo accesso alla “vera” distribuzione dei dati futuri. Questa è la sfida fondamentale della validazione predittiva: come possiamo valutare la qualità delle nostre previsioni senza conoscere la verità?

Per questo dobbiamo ricorrere a strategie di approssimazione, come la validazione incrociata e criteri predittivi come ELPD, che vedremo nelle prossime sezioni.

46.3.3.1 Mini-esempio intuitivo

Immagina una moneta truccata che dà testa il 70% delle volte.

Vogliamo confrontare due modelli:

  • Modello A: ipotizza una moneta equa (\(p = 0.5\)),
  • Modello B: ipotizza una probabilità leggermente sbilanciata (\(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 serve proprio a quantificare quanta informazione perdiamo quando ci affidiamo a un modello meno accurato (come A) invece che a uno più vicino alla verità (come B).

Il punto cruciale è che nella pratica non conosciamo mai la probabilità vera della moneta. Abbiamo soltanto i dati osservati, cioè gli esiti dei lanci. Per decidere quale modello predice meglio dobbiamo quindi basarci sui dati disponibili: è qui che entrano in gioco i metodi di confronto predittivo che studieremo.

46.4 Il log-score: misurare l’accuratezza predittiva punto per punto

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

46.4.1 La domanda di base

Per ogni osservazione nei nostri dati, vogliamo sapere: quanto il nostro modello considerava plausibile questo specifico risultato? Il log-score risponde proprio a questa domanda.

46.4.2 Definizione formale

Il log-score per un’osservazione \(y_i\) è semplicemente il logaritmo della probabilità che il modello assegnava a quell’osservazione:

\[ \text{Log-score}(y_i) = \log p(y_i \mid y) , \tag{46.3}\] dove \(p(y_i \mid y)\) è la distribuzione predittiva posteriore che abbiamo appena imparato a calcolare:

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

46.4.3 Interpretazione: “scommettere” sui dati

Il log-score può essere interpretato come quanto il modello avrebbe “scommesso” su quel particolare risultato:

  • se il modello assegna alta probabilità* a \(y_i\)*: \(p(y_i \mid y) \approx 1\), quindi \(\log p(y_i \mid y) \approx 0\) (buono);
  • se il modello assegna bassa probabilità* a \(y_i\)*: \(p(y_i \mid y) \approx 0\), quindi \(\log p(y_i \mid y) \ll 0\) (molto negativo, scarso).

Perché il logaritmo?
Il logaritmo trasforma prodotti di probabilità in somme. Così possiamo sommare contributi punto per punto dei dati invece di moltiplicarli; inoltre stabilizza i numeri molto piccoli tipici delle verosimiglianze.

46.4.4 Dal singolo dato al punteggio totale

Per avere una visione complessiva della performance del modello, sommiamo i log-score su tutte le osservazioni:

\[ S = \sum_{i=1}^n \log p(y_i \mid y) . \tag{46.4}\] Più \(S\) è alto (meno negativo), più il modello “scommette” bene sui dati osservati (in-sample).

46.4.5 Due filosofie a confronto: parametri fissi vs. incerti

Il calcolo del log-score può seguire due approcci concettualmente diversi, che riflettono due filosofie statistiche diverse.

46.4.5.1 Approccio classico: parametri fissi

Nell’impostazione frequentista classica, usiamo una stima puntuale dei parametri (ad es. Massima Verosimiglianza o MAP) e ignoriamo l’incertezza:

\[ \text{Log-score classico} = \log p(y_i \mid \hat{\theta}). \]

46.4.5.2 Approccio bayesiano: gestione dell’incertezza sui parametri

Nell’approccio bayesiano non fissiamo i parametri a un singolo valore stimato, ma li trattiamo come incerti. Questo significa che, invece di calcolare la verosimiglianza con un $$, “mescoliamo” tutte le verosimiglianze possibili, pesandole in base alla loro plausibilità a posteriori:

\[ \begin{align} \text{Log-score bayesiano} &= \log p(y_i \mid y) \\ &= \log \int p(y_i \mid \theta)\, p(\theta \mid y)\, d\theta . \end{align} \tag{46.5}\]

In altre parole, la probabilità predittiva di \(y_i\) non dipende da un solo \(\theta\), ma dalla media delle predizioni condizionate su tutti i valori plausibili dei parametri.

Differenza chiave
- L’approccio frequentista chiede: “Quanto sono plausibili i dati se i parametri valgono esattamente \(\hat{\theta}\)?”
- L’approccio bayesiano chiede: “Quanto sono plausibili i dati, in media, considerando tutti i valori di \(\theta\) compatibili con i dati osservati?”

La seconda prospettiva è più onesta perché riconosce l’incertezza sui parametri.

46.4.6 Implementazione pratica con il metodo MCMC

L’integrale nell’nell’Equazione 46.5 raramente ha una soluzione analitica (cioè non si può calcolare con una formula esatta). Possiamo però stimarlo in modo pratico ed efficiente utilizzando i campioni generati da un algoritmo MCMC (Markov Chain Monte Carlo).

Supponiamo di avere una serie di campioni di parametri, \(\theta^{(1)},\dots,\theta^{(S)}\), estratti dalla distribuzione a posteriori \(p(\theta\mid y)\). Il procedimento per approssimare la probabilità predittiva si articola in due semplici passi:

46.4.6.1 Passo 1: calcolare la verosimiglianza per ogni campione

Per ogni set di parametri \(\theta^{(s)}\) che abbiamo campionato, calcoliamo la probabilità (verosimiglianza) di osservare il dato \(y_i\) sotto quei parametri:

\[ p\bigl(y\_i \mid \theta^{(s)}\bigr). \]

Fare questo per tutti i campioni \(S\) ci fornisce un insieme di valori:

\[ \bigl\{\, p(y\_i \mid \theta^{(1)}),\; p(y\_i \mid \theta^{(2)}),\; \dots,\; p(y\_i \mid \theta^{(S)}) \,\bigr\} \] Questa collezione rappresenta come la plausibilità del dato \(y_i\) cambi al variare dei parametri, ponderata per la loro probabilità a posteriori.

46.4.6.2 Passo 2: calcolare la media dei valori ottenuti

La probabilità predittiva per \(y_i\) (che tiene conto di tutta l’incertezza sui parametri) è approssimata semplicemente calcolando la media aritmetica dell’insieme di valori ottenuti nel passo precedente:

\[ p(y_i \mid y) \;\approx\; \frac{1}{S}\sum_{s=1}^S p\bigl(y_i \mid \theta^{(s)}\bigr). \tag{46.6}\]

Il risultato è un singolo valore numerico (uno scalare) che sintetizza in una previsione probabilistica tutto ciò che abbiamo appreso sull’incertezza dei parametri del modello.

46.4.7 La LPPD: una misura bayesiana di bontà della previsione

Per valutare la capacità predittiva dell’intero modello su tutti i dati, ripetiamo il procedimento per ogni osservazione \(y_i\) e procediamo così:

  1. per ogni osservazione \(y_i\), calcoliamo la sua probabilità predittiva media \(p(y_i \mid y)\);
  2. prendiamo il logaritmo naturale di questa probabilità. (Usiamo il logaritmo per motivi computazionali e perché trasforma prodotti in somme);
  3. sommiamo i logaritmi di tutte le \(n\) osservazioni.

Il risultato di questo processo è la Log Pointwise Predictive Density (LPPD):

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

Confronto e Sintesi:

  • Il log-score classico (usato nella statistica frequentista) valuta la previsione utilizzando un unico valore stimato dei parametri (ad esempio, la stima di massima verosimiglianza \(\hat{\theta}\)). Questo ignora completamente l’incertezza esistente sulla stima dei parametri.
  • La LPPD bayesiana compie la stessa operazione fondamentale, ma in modo più robusto: invece di usare un singolo valore dei parametri, media le previsioni su tutte le migliaia di valori plausibili dei parametri campionati dalla distribuzione a posteriori. In questo modo, la misura di bontà predittiva incorpora in modo naturale tutta l’incertezza del modello.

46.4.8 Il problema nascosto: overfitting in-sample

La LPPD è calcolata sugli stessi dati usati per stimare il modello: modelli molto flessibili possono “scommettere bene” anche sul rumore, gonfiando la LPPD in-sample.

Analogia: È come valutare uno studente facendogli ripetere gli stessi esercizi che ha già risolto durante lo studio. Otterrà un punteggio alto, ma non sappiamo se sarebbe altrettanto bravo con problemi nuovi.

Per valutare la capacità di generalizzazione, serve una stima out-of-sample. Nelle prossime sezioni introdurremo la validazione incrociata leave-one-out (LOO-CV) e l’ELPD (Expected Log Pointwise Predictive Density), che forniscono una versione “fuori campione” della LPPD per il confronto predittivo tra modelli.

Consideriamo un singolo dato \(y_i = 3\) successi su \(n=5\) tentativi (Binomiale). Abbiamo tre valori plausibili per \(\theta\) dalla posterior, con pesi didattici \(w^{(s)}\) (nella pratica MCMC i pesi sono uguali):

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

Per ogni campione \(\theta^{(s)}\) calcoliamo \(p(y_i \mid \theta^{(s)})\), otteniamo la collezione di likelihood \(\{p(y_i\mid \theta^{(s)})\}_{s=1}^S\), poi facciamo la media pesata (eq. Equazione 46.6) per ottenere \(p(y_i\mid y)\), e infine il log-score \(\log p(y_i\mid y)\) (eq. Equazione 46.3).

# Dato osservato (un solo punto)
y_i  <- 3
n_i  <- 5

# "Campioni" posteriori (qui pochi e con pesi espliciti per chiarezza didattica)
theta_vals        <- c(0.3, 0.5, 0.7)      # θ^(1), θ^(2), θ^(3)
posterior_weights <- c(0.2, 0.5, 0.3)      # w^(1), w^(2), w^(3); in MCMC tipicamente uguali

# (1) Likelihood punto-per-punto: p(y_i | θ^(s))
likelihoods <- dbinom(y_i, size = n_i, prob = theta_vals)
likelihoods  # questa è la collezione { p(y_i | θ^(s)) }_s
#> [1] 0.132 0.312 0.309

# (2) Media (pesata) sulle likelihood ⇒ p(y_i | y) ≈ Σ_s w^(s) p(y_i | θ^(s))
p_yi_given_y <- sum(posterior_weights * likelihoods)

# (3) Log-score (per un solo dato coincide con la LPPD del singolo punto)
log_score_i <- log(p_yi_given_y)

# Stampa riassuntiva con notazione coerente
cat("Campioni θ^{(s)}:        ", theta_vals, "\n")
#> Campioni θ^{(s)}:         0.3 0.5 0.7
cat("Likelihood p(y_i|θ^{(s)}):", round(likelihoods, 4), "\n")
#> Likelihood p(y_i|θ^{(s)}): 0.132 0.312 0.309
cat("p(y_i|y) (media pesata):  ", round(p_yi_given_y, 4), "\n")
#> p(y_i|y) (media pesata):   0.275
cat("log p(y_i|y):             ", round(log_score_i, 4), "\n")
#> log p(y_i|y):              -1.29

Nota didattica: Nella pratica con MCMC i pesi sono uguali, \(w^{(s)}=\tfrac{1}{S}\), quindi \(p(y_i\mid y) \approx \tfrac{1}{S}\sum_{s=1}^S p(y_i\mid \theta^{(s)})\) (Equazione 46.6). Con più osservazioni \(\{y_i\}_{i=1}^n\), la LPPD è la somma dei log-score punto-per-punto (Equazione 46.7).

46.5 La svolta: dall’adattamento alla generalizzazione

46.5.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 che ha già visto molte volte durante l’allenamento.
  2. Scenario B: lo testi con nuove fotografie di persone mai viste prima.

Nel primo caso, lo studente probabilmente avrà un punteggio molto alto, ma non sapremo se ha davvero imparato a riconoscere le emozioni o se si limita a ricordare quelle immagini specifiche. Il secondo scenario, invece, è più onesto: misura la capacità di generalizzare la competenza a stimoli nuovi.

Lo stesso accade con i modelli statistici.

  • La LPPD corrisponde allo Scenario A: valuta il modello sugli stessi dati usati per adattarlo, rischiando di dare un’illusione di performance eccellente.
  • Per sapere se il modello sa davvero “generalizzare”, serve testarlo come nello Scenario B: con dati nuovi o tramite tecniche di validazione incrociata.

46.5.2 Guardare oltre i dati osservati

Quando valutiamo un modello, non ci basta sapere quanto bene spiega i dati che ha già visto. La vera domanda è: quanto bene predirebbe dati nuovi, mai osservati?

La Expected Log Predictive Density (ELPD) risponde a questa domanda. La logica è la stessa della LPPD, ma con una differenza cruciale: la previsione di ogni osservazione \(y_i\) viene fatta senza usare \(y_i\) per stimare il modello. Questa tecnica si chiama Leave-One-Out (LOO):

\[ \text{ELPD} = \sum_{i=1}^n \log p(y_i \mid y_{-i}), \tag{46.8}\] dove \(y_{-i}\) indica il dataset a cui è stata tolta l’osservazione \(i\).

46.5.2.1 Un esempio concreto per chiarire la differenza

Immagina di voler costruire un modello che predice i punteggi di memoria a breve termine degli studenti a partire dal loro livello di concentrazione.

  • Con la LPPD, il modello viene valutato sugli stessi studenti che sono serviti per stimarlo. È come dire: “quanto bene il modello spiega questi dati noti?”.
  • Con la ELPD, invece, ogni volta togliamo uno studente dal campione, stimiamo il modello sugli altri e proviamo a predire il punteggio di quello escluso. È come chiedere: “quanto bene il modello predirebbe un nuovo studente, mai visto prima?”.

46.5.2.2 Procedura passo per passo

  1. Prendiamo il primo studente ed escludiamolo dal dataset.
  2. Adattiamo il modello usando i dati dei rimanenti studenti.
  3. Prediciamo il punteggio di memoria dello studente escluso.
  4. Ripetiamo lo stesso procedimento per ogni studente, uno alla volta.
  5. Sommiamo tutti i log-score ottenuti: questo è l’ELPD.

46.6 Il collegamento con la divergenza KL

46.6.1 La teoria che unifica tutto

Abbiamo visto che l’ELPD fornisce una misura empirica della capacità predittiva di un modello. Esiste, tuttavia, una giustificazione teorica profonda e unificante che spiega il motivo per cui massimizzare l’ELPD è il principio corretto per la selezione dei modelli. Questa giustificazione poggia sul concetto di divergenza di Kullback-Leibler (KL).

46.6.1.1 Cosa misura la divergenza KL?

La divergenza KL, indicata come \(D_{\text{KL}}\), misura la “distanza” informazionale tra la distribuzione vera dei dati, \(p(\tilde{y})\) (la realtà che vogliamo catturare), e la distribuzione predittiva del nostro modello, \(q(\tilde{y} \mid y)\) (la nostra approssimazione). È definita come:

\[ 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\) è calcolata rispetto alla distribuzione vera \(p(\tilde{y})\).

46.6.1.2 Scomponiamo la divergenza KL

Per comprendere a fondo, espandiamo la definizione:

\[ D_{\text{KL}}(p \parallel q) = \underbrace{\mathbb{E}_p[\log p(\tilde{y})]}_{\text{(1) Entropia}} - \underbrace{\mathbb{E}_p[\log q(\tilde{y} \mid y)]}_{\text{(2) Accuratezza predittiva}}. \tag{46.9}\]

Analizziamo i due termini:

  1. \(\mathbb{E}_p[\log p(\tilde{y})]\) (Entropia): Rappresenta il contenuto informativo intrinseco della distribuzione vera. È una quantità fissa, immutabile e, soprattutto, indipendente dal modello che stiamo considerando. È una costante.
  2. \(-\mathbb{E}_p[\log q(\tilde{y} \mid y)]\) (Log-verosimiglianza attesa): Questo è il termine cruciale. Misura quanto è buona la nostra distribuzione predittiva \(q\) nel prevedere nuovi dati provenienti dalla vera distribuzione \(p\). Nota: questo è esattamente l’opposto della quantità che stimiamo con l’ELPD \((\sum \log q(\tilde{y} \mid y))\).

46.6.1.3 Il collegamento fondamentale

Poiché il primo termine dell’Equazione 46.9 è una costante, minimizzare la divergenza KL \(D_{\text{KL}}(p \parallel q)\) equivale esattamente a massimizzare il secondo termine, ovvero l’accuratezza predittiva attesa. Questo risultato si traduce in una regola pratica potentissima per il confronto tra modelli. Date due distribuzioni predittive, \(q_A\) e \(q_B\), la differenza nelle loro divergenze KL è:

\[ \begin{aligned} D_{\text{KL}}(p \parallel q_A) - D_{\text{KL}}(p \parallel q_B) &= \left( \cancel{\mathbb{E}_p[\log p(\tilde{y})]} - \mathbb{E}_p[\log q_A(\tilde{y} \mid y)] \right) \notag\\ &\qquad - \left( \cancel{\mathbb{E}_p[\log p(\tilde{y})]} - \mathbb{E}_p[\log q_B(\tilde{y} \mid y)] \right) \\ &= \mathbb{E}_p[\log q_B(\tilde{y} \mid y)] - \mathbb{E}_p[\log q_A(\tilde{y} \mid y)] \\ &= \text{ELPD}(q_B) - \text{ELPD}(q_A) \end{aligned} \] dove abbiamo cancellato il termine entropia costante.

46.6.1.4 Conclusione teorica fondamentale

Il risultato precedente ci porta alla conclusione chiave di tutta la teoria:

\[ \text{Massimizzare l'ELPD} \;\; \equiv \;\; \text{Minimizzare la divergenza KL dalla verità}. \] In altre parole, quando preferiamo il modello con l’ELPD più alto, non stiamo solo seguendo un criterio empirico. Stiamo scegliendo consapevolmente il modello la cui distribuzione predittiva è, in media, più vicina alla realtà sottostante in senso informazionale. Questo principio unifica la teoria dell’informazione con la pratica della valutazione e selezione dei modelli predittivi.

Vogliamo confrontare due modelli predittivi per il numero di “teste” in \(n=10\) lanci. Supponiamo che

  • la distribuzione vera è \(p(y)=\text{Binom}(n=10,\;p=0.6)\),
  • il modello candidato prevede \(q(y)=\text{Binom}(n=10,\;q=0.5)\).

L’ELPD di un modello è l’aspettativa, rispetto alla distribuzione vera \(p\), del log-score del modello: \(\mathrm{ELPD}(q)=\mathbb{E}_{p}[\log q(Y)]\). Nel caso discreto, l’aspettativa diventa una somma su tutti i possibili valori \(y=0,\dots,n\).

# Parametri del problema
n <- 10          # numero di lanci
p <- 0.6         # probabilità vera di "testa"
q <- 0.5         # probabilità ipotizzata dal modello candidato

# 1) Supporto dei possibili esiti
y_vals <- 0:n

# 2) Distribuzione vera p(y) su tutto il supporto
p_y <- dbinom(y_vals, size = n, prob = p)

# 3) Log-predittiva del modello candidato q su tutto il supporto
log_q_y <- log(dbinom(y_vals, size = n, prob = q))

# 4) ELPD del modello candidato: somma dei log q(y) pesati da p(y)
elpd_q <- sum(p_y * log_q_y)

# 5) "Modello vero": usa q = p. Log-predittiva del modello vero
log_p_y <- log(dbinom(y_vals, size = n, prob = p))

# 6) ELPD del modello vero: somma dei log p(y) pesati da p(y)
elpd_p <- sum(p_y * log_p_y)

# 7) Divergenza KL tra p e q: somma p(y) * log [p(y)/q(y)]
kl_pq <- sum(p_y * (log_p_y - log_q_y))

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

Cosa stiamo verificando?

  1. \(\mathrm{ELPD}(q)=\sum_y p(y)\log q(y)\) è più basso (più negativo) del valore ottenuto dal modello vero \(\mathrm{ELPD}(p)=\sum_y p(y)\log p(y)\). → Il modello con \(q=0.6\) è più predittivo di quello con \(q=0.5\).

  2. La differenza tra i due ELPD è uguale (vicina numericamente) alla divergenza di Kullback–Leibler:

\[ \begin{align} \mathrm{ELPD}(p)-\mathrm{ELPD}(q) &= \sum_y p(y)\big[\log p(y)-\log q(y)\big] \notag\\ &= D_{\mathrm{KL}}(p\|q)\;>\;0. \end{align} \] Questo mostra algebricamente e numericamente il legame: massimizzare l’ELPD equivale a minimizzare la divergenza KL.

In pratica: In questo esempio abbiamo potuto calcolare l’ELPD vero perché conoscevamo l’intera distribuzione generatrice \(p(y)\) e potevamo integrare esattamente. Nella realtà, \(p(y)\) è sconosciuta: disponiamo solo di un campione osservato. In questi casi stimiamo l’ELPD empiricamente, ad esempio con la Leave-One-Out Cross-Validation (LOO-CV), che sostituisce l’aspettativa rispetto a \(p\) con una media sui dati raccolti, lasciando fuori una osservazione alla volta.

Collegamento chiave.

L’ELPD è una stima empirica (con segno cambiato) della divergenza di Kullback–Leibler. Più alto è l’ELPD, migliore è la capacità predittiva del modello.

46.7 Stimare l’ELPD nella pratica: la Leave-One-Out Cross-Validation

Abbiamo chiarito che l’ELPD è la misura ideale della bontà predittiva di un modello, perché è direttamente collegata alla divergenza KL. Il problema è che, per definizione, richiede un’aspettativa rispetto alla vera distribuzione dei dati futuri \(p(\tilde{y})\), che non conosciamo mai.

Come possiamo allora stimarlo? La soluzione più usata è la Leave-One-Out Cross-Validation (LOO-CV), che ci permette di avvicinarci all’ELPD teorico usando soltanto i dati osservati.

46.7.1 L’idea alla base della LOO-CV

Il principio è semplice: trattare ogni osservazione del dataset come se fosse “nuova” e verificare se il modello, addestrato sui dati rimanenti, riesce a prevederla.

Il procedimento è questo:

  1. Si prende un’osservazione \(y_i\).
  2. La si rimuove temporaneamente dal dataset.
  3. Si stima il modello sui dati rimanenti \(y_{-i}\).
  4. Si calcola la densità predittiva che il modello assegna al dato escluso: \(p(y\_i \mid y_{-i})\).
  5. Si ripete per tutte le osservazioni e si sommano i logaritmi.

In formula:

\[ \text{ELPD}_{\text{LOO}} = \sum_{i=1}^n \log p(y_i \mid y_{-i}) . \tag{46.10}\]

Così otteniamo una stima out-of-sample: il modello viene valutato su dati che non ha mai visto.

46.7.2 Perché funziona

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

\[ \text{ELPD}_{\text{LOO}} \approx \mathbb{E}_p[\log q(\tilde{y}\mid y)] . \tag{46.11}\]

46.7.3 Confrontare i modelli con LOO-CV

Una volta stimato l’ELPD-LOO, possiamo confrontare due modelli calcolando la differenza:

\[ \Delta \text{ELPD} = \text{ELPD}_{\text{LOO}}(M_1) - \text{ELPD}_{\text{LOO}}(M_2). \tag{46.12}\]

Se la differenza è positiva, il modello \(M_1\) ha una distribuzione predittiva più vicina alla realtà di quella di \(M_2\).

È utile stimare anche un errore standard della differenza. Come regola empirica, una differenza di almeno due volte l’SE indica un vantaggio credibile di un modello sull’altro.

46.7.4 Overfitting e vantaggio della LOO-CV

Se valutassimo un modello sugli stessi dati usati per addestrarlo, la sua performance apparirebbe gonfiata (overfitting). La LOO-CV aggira questo problema: ogni osservazione viene valutata solo con modelli che non l’hanno vista. Il punteggio ottenuto è quindi una misura più realistica della capacità di generalizzare a nuovi dati.

46.7.5 PSIS-LOO: la scorciatoia pratica

Un limite della LOO tradizionale è che richiederebbe di riadattare il modello \(n\) volte, cosa spesso impraticabile. Per questo oggi si usa il metodo Pareto-Smoothed Importance Sampling (PSIS-LOO), che consente di stimare l’ELPD-LOO a partire da un unico adattamento del modello, sfruttando i campioni MCMC.

In R, tutto ciò è implementato nel pacchetto loo, già integrato in brms e rstanarm, attraverso funzioni come loo() e loo_compare(). Oltre ai valori di ELPD, queste funzioni forniscono anche diagnostiche (le Pareto k) per capire se la stima è affidabile.

46.7.5.1 In sintesi

  • L’ELPD misura la capacità predittiva del modello su dati futuri.
  • Non conoscendo la distribuzione vera, usiamo la LOO-CV per stimarlo.
  • La differenza di ELPD-LOO tra modelli approssima la differenza nelle loro divergenze KL.
  • PSIS-LOO rende il calcolo efficiente anche per modelli complessi.
  • La regola pratica: preferire il modello con ELPD-LOO più alto, tenendo conto anche della sua semplicità e interpretabilità.

Questo esempio mostra come passare dalla definizione teorica dell’ELPD alla stima pratica via Leave-One-Out, usando un caso elementare Beta–Bernoulli.

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

Modello A (Bayesiano adattato ai dati): Bernoulli\((\theta)\) con prior \(\theta\sim \text{Beta}(1,1)\) (uninformativa). Per LOO: * per ogni \(i\), escludiamo \(y_i\); * calcoliamo la posteriore \(\theta \mid y_{-i} \sim \text{Beta}(1+s_{-i},\,1+n_{-i}-s_{-i})\), dove \(s_{-i}\) è il numero di successi tra i \(n-1\) rimanenti; * calcoliamo la probabilità predittiva per \(y_i\).

Modello B (di confronto): Moneta equa fissa (\(q=0.5\)): la predittiva è sempre \(0.5\), indipendentemente dai dati.

# 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
  n_minus <- n - 1
  alpha <- a0 + s_minus
  beta  <- b0 + (n_minus - s_minus)
  p1 <- alpha / (alpha + beta)
  p  <- if (yi == 1) p1 else (1 - p1)
  log(p)
}

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

# ELPD-LOO
elpd_beta  <- sum(lp_beta)
elpd_fixed <- sum(lp_fixed)

# Differenza e SE
diff_pt <- lp_beta - lp_fixed
se_diff <- sqrt(n * var(diff_pt))

# Tabella riassuntiva
res <- data.frame(
  i = 1:n, y = y,
  lp_beta = round(lp_beta, 6),
  lp_fixed = round(lp_fixed, 6),
  diff = round(diff_pt, 6)
)
print(res)
#>   i y lp_beta lp_fixed   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_beta))
#> 
#> ELPD-LOO Modello A: -3.413620
cat(sprintf("ELPD-LOO Modello B: %.6f\n", elpd_fixed))
#> ELPD-LOO Modello B: -3.465736
cat(sprintf("Differenza (A-B)  : %.6f\n", elpd_beta - elpd_fixed))
#> Differenza (A-B)  : 0.052116
cat(sprintf("SE differenza     : %.6f\n", se_diff))
#> SE differenza     : 1.386294

Interpretazione:

  • Ogni riga della tabella mostra la log-predittiva fuori campione per entrambi i modelli
  • In un campione con 4 successi su 5, il Modello A assegna più di 0.5 di probabilità ai successi, e meno di 0.5 all’unico insuccesso
  • L’ELPD-LOO di A può risultare leggermente più alto di quello di B, ma l’errore standard è grande perché \(n\) è piccolo

Regola pratica: una differenza \(|\Delta \text{ELPD}|\) di almeno 2 volte l’SE fornisce un’indicazione più affidabile di superiorità del modello. In esempi così piccoli l’obiettivo è puramente didattico: capire come si calcola e cosa significa.

Concetto chiave

  • L’ELPD valuta la capacità predittiva su dati non visti.
  • La LOO-CV lo stima in modo efficiente con PSIS-LOO.

Strumenti

  • Funzione loo() del pacchetto loo, integrata in brms e rstanarm.
  • Diagnostica con Pareto k, confronto con loo_compare().

Workflow tipico in R

  1. Adattare ogni modello (brm() o stan_glm()).
  2. Estrarre log_lik() e calcolare loo().
  3. Confrontare modelli con loo_compare().

Decisione

  • Preferire l’ELPD-LOO più alto.
  • Differenza ≥ 2×SE → indicazione di vantaggio sostanziale.
  • Valutare anche semplicità e interpretabilità.

46.8 Criteri di informazione

Oltre alla convalida incrociata Leave-One-Out, la statistica offre altri strumenti per stimare la qualità predittiva di un modello senza conoscere la distribuzione vera dei dati. Molti di questi metodi affondano le loro radici teoriche nel concetto di divergenza di Kullback-Leibler, che misura la distanza tra la distribuzione generatrice dei dati e quella stimata dal nostro modello.

L’obiettivo comune è valutare la capacità di un modello di generalizzare, ovvero di fare buone previsioni su dati non osservati, senza farsi trarre in inganno dall’overfitting. Tutti i criteri seguono una logica simile, bilanciando due componenti: una misura della bontà di adattamento ai dati e una penalità per la complessità del modello stesso. I vari criteri si distinguono proprio per come definiscono queste due componenti e per le assunzioni su cui si basano.

46.8.1 Una panoramica dei criteri principali

L’AIC (Akaike Information Criterion) approssima la distanza di Kullback-Leibler utilizzando la verosimiglianza massimizzata e applica una penalità semplice, proporzionale al numero di parametri. È un criterio veloce e ampiamente utilizzato, particolarmente utile per modelli regolari con campioni non troppo piccoli e privi di struttura gerarchica.

Il BIC (Bayesian Information Criterion) segue una logica simile all’AIC, ma introduce una penalità per la complessità che cresce all’aumentare della dimensione del campione. Questo lo porta tendenzialmente a preferire modelli più parsimoniosi quando il numero di osservazioni è grande e, sotto specifiche ipotesi, può essere collegato alla verosimiglianza marginale.

Il WAIC (Widely Applicable Information Criterion) rappresenta una versione pienamente bayesiana. Utilizza l’intera distribuzione predittiva a posteriori per valutare il fit e stima una penalità per la complessità effettiva del modello, che può differire dal semplice numero di parametri. È particolarmente adatto per modelli complessi o non regolari ed è concettualmente molto vicino alla stima LOO.

Infine, il LOO-CV (Leave-One-Out Cross-Validation), specialmente nella sua efficiente implementazione PSIS-LOO, stima direttamente l’Expected Log Predictive Density (ELPD) escludendo un dato alla volta. È spesso considerato il gold standard per il confronto predittivo nell’ambito della modellazione bayesiana, grazie alla sua robustezza e alle utili diagnostiche che fornisce.

46.8.1.1 Come orientarsi nella scelta

Una regola pratica è che se l’obiettivo principale è la previsione fuori campione in un contesto bayesiano, il PSIS-LOO o il WAIC sono generalmente da preferire ad AIC e BIC. In un approccio frequentista classico, con modelli regolari e campioni di dimensioni medio-grandi, l’AIC rimane un buon compromesso, mentre il BIC può essere più appropriato quando si desidera enfatizzare la parsimonia.

Per modelli bayesiani con obiettivo predittivo e dati reali (spesso non iid o gerarchici), il PSIS-LOO è la prima scelta, con il WAIC utile come riscontro. Con campioni piccoli, strutture complesse o unità dipendenti, è bene evitare criteri puramente asintotici come AIC e BIC, preferendo invece LOO o definendo con attenzione l’unità di esclusione (ad esempio, per soggetto o per gruppo). Nei modelli gerarchici o multilivello, LOO e WAIC possono essere applicati in modo coerente, prestando attenzione a non escludere singole osservazioni se queste non sono indipendenti, ma piuttosto interi cluster.

46.8.1.2 Errori comuni e best practice

Un errore frequente è utilizzare il Mean Squared Error (MSE) sul campione di addestramento come metro di giudizio, poiché questo valore non penalizza la complessità e tende quindi a favorire modelli eccessivamente flessibili e soggetti a overfitting. Allo stesso modo, è importante ricordare che AIC e BIC si basano su stime puntuali (MLE o MAP) e non catturano l’incertezza completa sui parametri, il che li rende meno ideali in un contesto bayesiano puro. WAIC e LOOCV, al contrario, sono espressamente concepiti per stimare la performance predittiva su dati nuovi.

Quando si riporta un confronto tra modelli, è buona norma includere non solo il modello “vincente”, ma anche la differenza di ELPD con il suo errore standard, le diagnostiche sui parametri di Pareto-k, una stima della complessità effettiva e un commento sostantivo che spieghi il motivo della preferenza, che potrebbe risiedere nella parsimonia, nell’interpretabilità dei parametri o nella robustezza.

46.8.1.3 In sintesi: il workflow essenziale

Un mini-workflow consigliato per un approccio bayesiano prevede di: adattare i modelli; calcolare il LOO per ciascuno di essi e controllare i parametri di Pareto-k; se si riscontrano valori di k elevati, considerare una convalida incrociata K-fold o una LOO per cluster; confrontare i modelli con appositi strumenti e riportare le differenze di ELPD; opzionalmente, calcolare il WAIC come controllo incrociato; argomentare infine la scelta finale anche in base a parsimonia e interpretabilità.

La selezione del modello, in definitiva, ruota attorno a una domanda essenziale: quanto bene il modello predice dati che non ha mai visto? Il riferimento teorico è l’Expected Log Predictive Density (ELPD), che misura quanto la distribuzione predittiva del modello si avvicina alla vera distribuzione dei dati. Poiché quest’ultima è sconosciuta, l’ELPD va stimato con strumenti come LOO-CV e WAIC, che oggi rappresentano gli standard più affidabili per guidare una scelta consapevole, equilibrata e focalizzata sulla capacità di generalizzazione.

Riflessioni conclusive

Il principio fondamentale della modellazione bayesiana risiede nella valutazione della qualità di un modello attraverso la sua capacità di produrre previsioni probabilistiche robuste, rappresentate dalla distribuzione predittiva a posteriori \(p(\tilde{y} \mid y)\).

La misura che guida questa valutazione è l’Expected Log Predictive Density (ELPD), che quantifica la capacità predittiva del modello su dati non osservati. A differenza delle metriche in-sample, soggette a sovradattamento, l’ELPD fornisce una stima imparziale della capacità di generalizzazione. Teoricamente, massimizzare l’ELPD equivale a minimizzare la divergenza di Kullback-Leibler rispetto alla vera distribuzione generatrice dei dati.

Operativamente, l’ELPD viene stimato mediante PSIS-LOO, integrato con i diagnostici Pareto-k. Il WAIC costituisce un’alternativa bayesiana solida, spesso coerente con LOO. Al contrario, criteri come AIC e BIC, sebbene computazionalmente efficienti, si basano su stime puntuali e approssimazioni asintotiche, risultando meno affidabili in contesti di campioni piccoli o modelli gerarchici.

Nel confronto tra modelli, è essenziale riportare non solo l’ELPD-LOO, ma anche le differenze ΔELPD e i relativi errori standard. Tuttavia, la selezione del modello non dovrebbe ridursi a un esercizio meccanico: differenze trascurabili nell’ELPD, specialmente se associate ad alta incertezza, possono essere irrilevanti sul piano sostanziale. Modelli meno performanti ma più parsimoniosi o teoricamente fondati possono rappresentare scelte migliori.

L’obiettivo finale è bilanciare capacità predittiva e coerenza teorica, ricordando che lo scopo della modellazione non è solo prevedere, ma comprendere. La valutazione deve quindi integrare strumenti come il PSIS-LOO con considerazioni sull’incertezza statistica, la struttura dei dati e il contesto teorico di riferimento.

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.0.0              pillar_1.11.1         tinytable_0.13.0     
#>  [4] patchwork_1.3.2       ggdist_3.3.3          tidybayes_3.0.7      
#>  [7] bayesplot_1.14.0      ggplot2_4.0.0         reliabilitydiag_0.2.1
#> [10] priorsense_1.1.1      posterior_1.6.1       loo_2.8.0            
#> [13] rstan_2.32.7          StanHeaders_2.32.10   brms_2.23.0          
#> [16] Rcpp_1.1.0            sessioninfo_1.2.3     conflicted_1.2.0     
#> [19] janitor_2.2.1         matrixStats_1.5.0     modelr_0.1.11        
#> [22] tibble_3.3.0          dplyr_1.1.4           tidyr_1.3.1          
#> [25] rio_1.2.3             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        systemfonts_1.2.3    
#> [10] vctrs_0.6.5           stringr_1.5.2         pkgconfig_2.0.3      
#> [13] arrayhelpers_1.1-0    fastmap_1.2.0         backports_1.5.0      
#> [16] labeling_0.4.3        rmarkdown_2.29        ragg_1.5.0           
#> [19] purrr_1.1.0           xfun_0.53             cachem_1.1.0         
#> [22] jsonlite_2.0.0        broom_1.0.10          parallel_4.5.1       
#> [25] R6_2.6.1              stringi_1.8.7         RColorBrewer_1.1-3   
#> [28] lubridate_1.9.4       estimability_1.5.1    knitr_1.50           
#> [31] zoo_1.8-14            Matrix_1.7-4          splines_4.5.1        
#> [34] timechange_0.3.0      tidyselect_1.2.1      abind_1.4-8          
#> [37] codetools_0.2-20      curl_7.0.0            pkgbuild_1.4.8       
#> [40] lattice_0.22-7        withr_3.0.2           bridgesampling_1.1-2 
#> [43] S7_0.2.0              coda_0.19-4.1         evaluate_1.0.5       
#> [46] survival_3.8-3        RcppParallel_5.1.11-1 xml2_1.4.0           
#> [49] tensorA_0.36.2.1      checkmate_2.3.3       stats4_4.5.1         
#> [52] distributional_0.5.0  generics_0.1.4        rprojroot_2.1.1      
#> [55] rstantools_2.5.0      scales_1.4.0          xtable_1.8-4         
#> [58] glue_1.8.0            emmeans_1.11.2-8      tools_4.5.1          
#> [61] mvtnorm_1.3-3         grid_4.5.1            QuickJSR_1.8.0       
#> [64] colorspace_2.1-1      nlme_3.1-168          cli_3.6.5            
#> [67] textshaping_1.0.3     svUnit_1.0.8          Brobdingnag_1.2-9    
#> [70] V8_7.0.0              gtable_0.3.6          digest_0.6.37        
#> [73] TH.data_1.1-4         htmlwidgets_1.6.4     farver_2.1.2         
#> [76] memoise_2.0.1         htmltools_0.5.8.1     lifecycle_1.0.4      
#> [79] MASS_7.3-65

Bibliografia

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