77  Valutare i modelli bayesiani: Log-Score, LPPD, ELPD e LOO-CV

“Se uno scienziato è costretto a scegliere fra due ipotesi, il suo istinto sarà quello di scegliere la più semplice. Questa inclinazione non è un capriccio, ma una necessità logica.”

– E.T. Jaynes, Fisico e Statistico

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 bene ai dati già osservati, ma anche saper generalizzare a situazioni future o a campioni mai visti. Questa distinzione — adattamento vs. generalizzazione — è il cuore della valutazione predittiva.

Per rendere concreta questa idea, immaginiamo di aver sviluppato un test psicologico per prevedere il livello di ansia degli studenti alla vigilia di un esame. Non basta sapere che il modello descrive bene i dati del campione che abbiamo usato per costruirlo: vogliamo anche essere ragionevolmente sicuri che le stesse previsioni funzionino per studenti che non hanno partecipato allo studio. In psicologia, scegliere tra due modelli non è diverso dal decidere quale test usare per prevedere un disturbo: entrambi mirano a capire quale strumento fornisce previsioni più affidabili sui dati futuri.

In questo capitolo esploreremo gli strumenti fondamentali per la valutazione e il confronto di modelli nell’ambito dell’inferenza bayesiana.

1. La distribuzione predittiva posteriore
Introdurremo la distribuzione predittiva posteriore, che incorpora l’incertezza sui parametri per generare previsioni coerenti con lo stato di conoscenza del modello. Questo strumento rappresenta il ponte naturale tra stima e previsione, garantendo una quantificazione probabilistica completa dell’incertezza.

2. Misure di accuratezza predittiva
Discuteremo il log-score, una metrica punto per punto che valuta la qualità delle previsioni, e due sue sintesi fondamentali:

• la LPPD (Log Pointwise Predictive Density), che misura la bontà di adattamento su dati osservati;
  
• l’ELPD (Expected Log Predictive Density), che stima l’abilità predittiva attesa su nuove osservazioni.

3. Validazione empirica e confronto tra modelli
Presenteremo la tecnica Leave-One-Out Cross-Validation (LOO-CV), un approccio efficiente per stimare l’ELPD senza bisogno di nuovi dati, dimostrando come questa metodologia fornisca una valutazione robusta delle prestazioni predittive.

4. Fondamenti teorici e interpretazione
Approfondiremo il legame tra ELPD e divergenza di Kullback-Leibler, che consente di interpretare il confronto tra modelli come una ricerca del modello più vicino alla vera distribuzione generatrice dei dati. Questa connessione teorica fornisce una solida giustificazione informazionale per le procedure di selezione bayesiana.

L’obiettivo del capitolo è offrire una panoramica completa e operativa, che unisca principi teorici a strumenti applicativi, guidando il lettore nella scelta razionale del modello più adatto al problema in esame.

Panoramica del capitolo

• Cos’è la distribuzione predittiva posteriore e come si costruisce.
• Cosa misura il log-score e come si calcola nella pratica.
• Distinzione tra LPPD ed ELPD e il loro significato;
• Come LOO-CV fornisce una stima dell’ELPD;
• Il confronto tra modelli alla divergenza di Kullback-Leibler.

---

Prerequisiti

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

Preparazione del Notebook

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

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

77.1 Distribuzione predittiva posteriore

Nel capitolo precedente abbiamo usato la divergenza di Kullback–Leibler (KL) come misura teorica della distanza tra realtà e modello. Qui ci chiediamo: come stimiamo questa distanza quando la “vera” distribuzione generatrice è ignota? Un tassello fondamentale è la distribuzione predittiva posteriore.

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.

In sintesi: 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. 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à posteriore.

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

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

Intuizione

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.

Esempio numerico

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)
pred_df
#>     x         p
#> 1   0 1.139e-01
#> 2   1 2.506e-01
#> 3   2 2.762e-01
#> 4   3 1.995e-01
#> 5   4 1.040e-01
#> 6   5 4.069e-02
#> 7   6 1.206e-02
#> 8   7 2.662e-03
#> 9   8 4.178e-04
#> 10  9 4.200e-05
#> 11 10 2.049e-06

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.11391    0.11391
#> [2,] 0.25061    0.25061
#> [3,] 0.27618    0.27618
#> [4,] 0.19946    0.19946
#> [5,] 0.10398    0.10398
#> [6,] 0.04069    0.04069

max(abs(pred_df$p - bb_pmf))   # lo scarto massimo (dovrebbe essere ~0 con J grande)
#> [1] 1.679e-15

Notazione. Useremo talvolta la forma compatta \(q(\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.

77.1.1 Il problema della valutazione predittiva

Il nostro obiettivo è capire quanto la distribuzione predittiva posteriore \(q(\tilde{y} \mid y)\) si avvicini alla vera distribuzione generatrice dei dati futuri, \(p(\tilde{y})\). In teoria, questa distanza si misura con la divergenza di Kullback–Leibler (KL):

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

Qui però incontriamo subito un problema concettuale: non conosciamo \(p(\tilde{y})\). Per superare questo ostacolo, possiamo ricorrere a misure surrogate che, pur non avendo accesso diretto a \(p(\tilde{y})\), permettono di stimare la qualità predittiva del modello utilizzando in modo ingegnoso i dati osservati. Tra queste, vedremo il log-score, la LPPD e l’ELPD, che forniscono stime indirette della bontà predittiva.

Mappa 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
\(q(\tilde{y} \mid y)\)	Distribuzione predittiva complessiva	Divergenza KL, confronto modelli

77.2 Il log-score: accuratezza predittiva punto per punto

Abbiamo definito la distribuzione predittiva posteriore. Ora chiediamoci: quanto bene il modello ha previsto ciascun dato osservato? Il log-score risponde proprio a questa domanda: per ogni osservazione \(y_i\) misura quanto il modello la considerava plausibile, cioè quanto avrebbe scommesso su quel dato.

Formalmente,

\[ \log p(y_i \mid y) \;=\; \log \int p(y_i \mid \theta)\, p(\theta \mid y)\, d\theta . \tag{77.1}\]

Se il modello assegna alta probabilità a \(y_i\), \(\log p(y_i \mid y)\) è vicino a 0 (buono). Se assegna bassa probabilità, il log-score è molto negativo (scarso).

Perché il logaritmo?
Il log 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.

77.2.1 Dal singolo dato al punteggio totale

Per avere una visione complessiva, sommiamo i contributi su tutte le osservazioni:

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

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

77.2.2 Parametri fissati vs. parametri incerti

Ci sono due modi concettualmente distinti per valutare il log-score.

Parametri fissati (impostazione classica). Usiamo una stima puntuale dei parametri (ad es. Massima Verosimiglianza o MAP) e ignoriamo l’incertezza:

\[ \log p(y_i \mid \hat{\theta}) . \]

Parametri incerti (impostazione bayesiana). Non fissiamo \(\theta\), ma lo trattiamo come incerto e “mescoliamo” le verosimiglianze pesandole per la plausibilità a posteriori:

\[ p(y_i \mid y) \;=\; \int p(y_i \mid \theta)\, p(\theta \mid y)\, d\theta . \tag{77.3}\]

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

77.2.3 Come stimare l’integrale in pratica: campioni MCMC

L’integrale nell’Equazione 77.3 raramente è calcolabile in forma chiusa. Con i campioni MCMC \(\theta^{(1)},\dots,\theta^{(S)} \sim p(\theta\mid y)\) possiamo approssimarlo così:

1. Per ciascun campione \(\theta^{(s)}\) calcoliamo la verosimiglianza del dato \(y_i\): \[ p\bigl(y_i \mid \theta^{(s)}\bigr). \] Questo produce una collezione di valori \[ \bigl\{\, p(y_i \mid \theta^{(1)}),\; p(y_i \mid \theta^{(2)}),\; \dots,\; p(y_i \mid \theta^{(S)}) \,\bigr\}, \] che rappresenta come la plausibilità di \(y_i\) varia al variare dei parametri plausibili.

2. Media sui campioni (mixing).
   La probabilità predittiva puntuale di \(y_i\), che useremo nel log-score, è la media di quella collezione: \[ p(y_i \mid y) \;\approx\; \frac{1}{S}\sum_{s=1}^S p\bigl(y_i \mid \theta^{(s)}\bigr). \tag{77.4}\] Questa media è uno scalare: condensa l’incertezza sui parametri in un’unica previsione probabilistica per \(y_i\).

Mini-illustrazione: se per tre campioni otteniamo \(\{0.40, 0.50, 0.60\}\), la media è \(0.50\). Questo numero è \(p(y_i \mid y)\) da inserire nel log.

77.2.4 La LPPD: il log-score bayesiano complessivo

Ripetiamo i passi precedenti per ogni osservazione \(y_i\):

1. calcoliamo la probabilità predittiva media \(p(y_i\mid y)\) con l’Equazione 77.4;
   
2. ne prendiamo il logaritmo;
   
3. sommiamo su tutte le osservazioni.

Il risultato è 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{77.5}\]

In sintesi, il log-score classico usa un solo valore dei parametri \((\hat{\theta})\); la LPPD compie lo stesso calcolo ma tiene conto dell’incertezza, mediando su tutti i valori plausibili secondo la posterior.

77.2.5 Attenzione all’overfitting

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. 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.

Esempio.

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 (Equazione 77.4) per ottenere \(p(y_i\mid y)\), e infine il log-score \(\log p(y_i\mid y)\) (Equazione 77.1).

# 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.1323 0.3125 0.3087

# (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.1323 0.3125 0.3087
cat("p(y_i|y) (media pesata):  ", round(p_yi_given_y, 4), "\n")
#> p(y_i|y) (media pesata):   0.2753
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 77.4). Con più osservazioni \(\{y_i\}_{i=1}^n\), la LPPD è la somma dei log-score punto-per-punto (Equazione 77.5).

77.2.6 Expected Log Predictive Density (ELPD): guardare oltre i dati osservati

Se vogliamo valutare la capacità di generalizzazione di un modello, la domanda chiave è: quanto bene predirebbe dati che non ha mai visto? L’ELPD (Expected Log Predictive Density) risponde a questa domanda con la stessa logica della LPPD, ma introduce una differenza fondamentale: la previsione di \(y_i\) viene calcolata escludendo \(y_i\) dall’adattamento del modello (Leave-One-Out, LOO):

\[ \text{ELPD} \;=\; \sum_{i=1}^n \log p(y_i \mid y_{-i}), \tag{77.6}\]

dove \(y_{-i}\) indica il dataset a cui è stata rimossa l’osservazione \(i\).

Esempio Nel caso di un test sull’ansia:

• LPPD → misura quanto bene il modello predice i punteggi di ansia degli studenti già presenti nel campione osservato.
• ELPD → misura quanto bene predirebbe il punteggio di un nuovo studente, usando solo i dati degli altri.

In sostanza, l’ELPD è una stima empirica (con segno cambiato) della divergenza di Kullback–Leibler tra la vera distribuzione dei dati futuri e la distribuzione predittiva del modello. Ci fornisce quindi un indicatore diretto di quanto le previsioni del modello si avvicinano a ciò che accadrà davvero, senza richiedere di conoscere la distribuzione reale.

Interpretazione: l’ELPD è un log-score out-of-sample: per ogni \(y_i\), lo escludiamo, adattiamo il modello agli altri dati, e valutiamo la probabilità predittiva di \(y_i\). Più alto è l’ELPD, migliore è la capacità del modello di generalizzare a dati nuovi.

Esempio.

Supponiamo di avere tre osservazioni \(y_1, y_2, y_3\) e che il modello stimi:

\[ p(y_1 \mid y_2,y_3)=0.6,\quad p(y_2 \mid y_1,y_3)=0.7,\quad p(y_3 \mid y_1,y_2)=0.5. \]

L’ELPD è:

\[ \log 0.6 + \log 0.7 + \log 0.5 \; \approx\; -0.5108 -0.3567 -0.6931 = -1.5606. \]

Un valore meno negativo indica maggiore capacità predittiva fuori campione.

77.2.7 LPPD vs. ELPD in sintesi

Misura	Dati usati per predire \(y_i\)	Valuta	Limite principale
LPPD	Tutti i dati, incluso \(y_i\)	Adattamento in-sample	Rischio di overfitting
ELPD	Tutti i dati tranne \(y_i\) (LOO)	Generalizzazione	—

Metafora In un esperimento di riconoscimento di volti, mostriamo a un partecipante 100 fotografie e lo alleniamo a riconoscerle:

• LPPD → misura quanto bene riconosce quelle stesse foto, già viste in fase di addestramento (in-sample).
• ELPD → misura quanto bene riconosce nuove foto, mai viste prima, cioè immagini fuori dall’insieme di addestramento (out-of-sample).

Se il punteggio LPPD è alto ma l’ELPD è basso, significa che il partecipante — o il modello — ha memorizzato i casi specifici, senza aver appreso regole generali utili per nuovi dati.

77.2.8 Il collegamento con la divergenza KL

La divergenza di Kullback–Leibler \(D_{\text{KL}}\) misura teoricamente la distanza tra la distribuzione vera dei dati, \(p(\tilde{y})\), e la distribuzione predittiva del modello, \(q(\tilde{y} \mid y)\).

Nel confronto tra due modelli \(A\) e \(B\), la differenza nelle loro \(D_{\text{KL}}\) equivale alla differenza nelle rispettive accuratezze predittive medie rispetto a \(p(\tilde{y})\).

Poiché \(p(\tilde{y})\) è sconosciuta, non possiamo calcolare direttamente la KL. L’ELPD fornisce una stima empirica di questa accuratezza predittiva: un valore più alto implica un modello più “vicino” alla distribuzione vera.

\[ \text{Massimizzare ELPD} \;\; \approx \;\; \text{Minimizzare la divergenza KL}. \]

Perché ELPD ≈ - KL?

Per definizione:

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

• Il primo termine \(\mathbb{E}_{p}[\log p(\tilde{y})]\) non dipende dal modello (è fisso per tutti).
• Confrontare due modelli equivale quindi a confrontare solo il secondo termine, che è \(-\)ELPD.

Ecco perché massimizzare l’ELPD equivale a minimizzare la divergenza KL: si sta massimizzando la media log-predittiva che il modello assegna ai dati futuri.

Esempio.

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

• 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:

\[ \mathrm{ELPD}(p)-\mathrm{ELPD}(q) = \sum_y p(y)\big[\log p(y)-\log q(y)\big] = D_{\mathrm{KL}}(p\|q)\;>\;0. \]

→ Questo mostra algebricamente e numericamente il legame: massimizzare l’ELPD equivale a minimizzare la KL.

Nota sul log: nel codice usiamo il log naturale (unità in nat). Se si preferisce il log in base 2 (unità in bit), basta sostituire log() con log2(); tutte le quantità cambiano di una costante di scala, ma i confronti tra modelli restano identici.

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. Questa procedura ci consente di avvicinarci al calcolo ideale della KL, anche senza conoscere \(p(y)\).

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.

77.3 Leave-One-Out Cross-Validation (LOO-CV): stimare l’ELPD nella pratica

Poiché la distribuzione vera dei dati futuri è inaccessibile, dobbiamo usare metodi indiretti per stimare quanto bene il nostro modello prevede nuove osservazioni. La validazione incrociata Leave-One-Out (LOO-CV) è uno di questi metodi e, se combinata con l’uso dell’Expected Log Predictive Density (ELPD), diventa uno strumento potente per il confronto tra modelli.

Abbiamo visto che l’ELPD è la misura ideale della capacità predittiva di un modello su dati futuri. Il problema è che, per definizione, richiede di calcolare un’aspettativa rispetto alla vera distribuzione generatrice \(p(\tilde{y})\), che non conosciamo.

Come possiamo stimarla in pratica? Usando la LOO-CV, che simula la previsione di nuovi dati sfruttando solo le informazioni presenti nei dati osservati.

77.3.1 Cos’è la LOO-CV

La LOO-CV è un esperimento concettuale semplice:

1. Scegli un’osservazione \(y_i\) dal dataset.
2. Escludila dal set di addestramento.
3. Adatta il modello ai dati rimanenti \(y_{-i}\).
4. Calcola la densità predittiva del modello per l’osservazione esclusa: \(p(y_i \mid y_{-i})\).
5. Ripeti per ogni osservazione e somma i logaritmi ottenuti.

Formalmente:

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

dove \(y_{-i}\) indica il dataset senza l’osservazione \(i\).
La struttura è identica a quella dell’ELPD “ideale”, ma ogni termine è calcolato fuori campione, escludendo il dato che viene valutato.

Un’analogia: è come escludere uno studente dall’allenamento e verificare se il modello riesce a predire il suo punteggio d’esame; ripetendo questo processo per tutti gli studenti otteniamo una misura diretta della capacità di generalizzazione.

77.3.2 Perché LOO-CV funziona

L’ELPD può essere scritto come:

\[ \mathbb{E}_p[\log q(\tilde{y} \mid y)], \tag{77.8}\]

dove \(q(\tilde{y} \mid y)\) è la distribuzione predittiva del modello.
Non possiamo calcolare l’aspettativa rispetto a \(p(\tilde{y})\), ma possiamo trattare ogni osservazione \(y_i\) come “nuovo dato” generato da \(p\) e usare la media empirica sulle osservazioni reali come stima dell’aspettativa:

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

In altre parole: LOO-CV misura quanto bene il modello predirebbe ciascun dato se non lo avesse mai visto.

77.3.3 Legame con la divergenza KL

La divergenza di Kullback–Leibler è definita come:

\[ D_{\text{KL}}(p \parallel q) = \mathbb{E}_p[\log p(\tilde{y})] - \mathbb{E}_p[\log q(\tilde{y} \mid y)]. \]

Il primo termine, l’entropia di \(p\), è lo stesso per tutti i modelli e scompare nel confronto.
Ne segue che, per due modelli \(q_1\) e \(q_2\):

\[ D_{\text{KL}}(p \parallel q_1) - D_{\text{KL}}(p \parallel q_2) = \mathbb{E}_p[\log q\_2(\tilde{y} \mid y)] - \mathbb{E}_p[\log q_1(\tilde{y} \mid y)]. \]

Vince il modello con ELPD più alto, perché corrisponde alla minore divergenza KL dalla distribuzione vera.

77.3.4 Confrontare i modelli con LOO-CV

Poiché \(p(\tilde{y})\) è sconosciuta, sostituiamo l’aspettativa teorica con la stima empirica via LOO:

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

\(\Delta\text{ELPD}\) approssima la differenza tra le divergenze KL dei modelli.
Oltre alla differenza, possiamo stimare un errore standard per capire se la superiorità di un modello è robusta o dovuta al caso.

77.3.5 Punti chiave

• Problema: L’ELPD teorico richiede \(p(\tilde{y})\), che è sconosciuta.
  
• Soluzione: LOO-CV fornisce una stima empirica out-of-sample.
  
• Teoria: L’ELPD è direttamente collegato alla parte “accuratezza” della KL-divergence.
  
• Pratica: Massimizzare l’ELPD stimato equivale a scegliere il modello più vicino alla distribuzione vera.

Direi che l’esempio che hai scritto è già molto chiaro e in linea con il testo precedente, ma per integrarlo meglio nel capitolo e mantenere continuità con la sezione teorica, potremmo:

1. Aggiungere un’introduzione contestuale per collegarlo subito alla discussione ELPD–LOO–KL.
2. Rendere più esplicito il parallelismo con la teoria (ELPD come somma delle log-predittive fuori campione).
3. Sintetizzare il codice con commenti chiave, così che lo studente possa leggerlo senza perdersi nei dettagli secondari.
4. Chiarire il senso della tabella subito dopo l’esecuzione del codice.

Esempio: confronto ELPD-LOO tra due modelli

Questo mini-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.4055  -0.6931  0.2877
#> 2 2 1 -0.4055  -0.6931  0.2877
#> 3 3 1 -0.4055  -0.6931  0.2877
#> 4 4 0 -1.7918  -0.6931 -1.0986
#> 5 5 1 -0.4055  -0.6931  0.2877
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.

77.3.6 ELPD-LOO e il problema dell’overfitting

Valutare un modello sugli stessi dati usati per addestrarlo tende a gonfiare le stime della sua capacità predittiva (overfitting). È come se uno studente ottenesse un punteggio perfetto ripetendo esercizi già svolti: non sappiamo se saprebbe risolverne di nuovi.

La Leave-One-Out Cross-Validation (LOO-CV) aggira il problema valutando ciascuna osservazione \(y_i\) usando solo i dati rimanenti (\(y_{-i}\)). Il punteggio ottenuto (ELPD-LOO) è quindi una stima out-of-sample della bontà predittiva, meno sensibile all’overfitting.

Grazie a metodi come il Pareto-smoothed importance sampling (PSIS), oggi è possibile calcolare l’ELPD-LOO senza riadattare il modello \(n\) volte. In R, la funzione loo() del pacchetto loo (integrata in brms e rstanarm) rende questa procedura rapida e diretta anche per modelli complessi.

Esempio: ELPD atteso vs ELPD-LOO stimato

Quando conosci la distribuzione vera dei dati (\(p\)), puoi calcolare l’ELPD atteso in modo esatto. Quando invece hai solo i dati osservati, lo stimi tramite Leave-One-Out (PSIS-LOO), come mostrato di seguito.

ELPD atteso (p noto)

n <- 10
p <- 0.6
q <- 0.5

y_vals   <- 0:n
p_y      <- dbinom(y_vals, size = n, prob = p)
log_q_y  <- log(dbinom(y_vals, size = n, prob = q))
elpd     <- sum(p_y * log_q_y)

cat(sprintf("ELPD atteso (modello q = 0.5): %.4f\n", elpd))
#> ELPD atteso (modello q = 0.5): -2.0549

ELPD-LOO stimato da dati simulati

set.seed(123)
df <- data.frame(
  k = rbinom(20, size = n, prob = p),
  n = n
)

fit <- brm(k | trials(n) ~ 1, data = df,
           family = binomial(),
           prior = prior(constant(0), class = "Intercept"),
           iter = 2000, chains = 2)

loo_res <- loo(fit)
print(loo_res)

L’oggetto loo_res fornisce l’ELPD-LOO, il suo errore standard, e le statistiche Pareto k per la diagnostica. Con loo_compare() puoi confrontare due modelli sulla base della differenza di ELPD-LOO e del relativo SE.

In pratica: stimare e confrontare l’ELPD-LOO

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à.

77.4 Criteri di informazione come approssimazioni della divergenza \(D_{\text{KL}}\)

Oltre alla Leave-One-Out Cross-Validation, esistono altri strumenti per stimare la qualità predittiva di un modello senza dover conoscere la distribuzione vera dei dati. Molti di questi metodi derivano, in modo più o meno diretto, dalla divergenza di Kullback–Leibler \(D_{\text{KL}}\), che — come visto — misura la distanza tra la distribuzione reale e quella stimata dal modello.

L’idea di base è sempre la stessa:

• valutare quanto bene il modello spiega i dati (bontà di adattamento);
• penalizzare la complessità del modello, per ridurre il rischio di overfitting.

Questa logica si traduce in criteri di informazione che combinano due componenti:

1. termine di fit: misura di quanto bene il modello si adatta ai dati osservati (es. log-verosimiglianza, MSE);
2. termine di penalizzazione: aumenta con il numero di parametri o con la flessibilità del modello.

Tra i criteri più usati troviamo:

• MSE (Mean Squared Error) – semplice e intuitivo, basato sugli errori di previsione;
• AIC (Akaike Information Criterion) – approssima \(D_{\text{KL}}\) tra il modello e la verità, penalizzando il numero di parametri;
• BIC (Bayesian Information Criterion) – simile all’AIC, ma con penalizzazione più forte per modelli complessi, proporzionale al numero di osservazioni;
• WAIC (Widely Applicable Information Criterion) – versione pienamente bayesiana, basata sulle previsioni del modello integrate sull’intera distribuzione a posteriori.

Nelle sezioni seguenti vedremo come ciascun criterio si calcola, quali assunzioni richiede e in quali situazioni è preferibile rispetto agli altri.

77.4.1 Errore Quadratico Medio (MSE)

L’Errore Quadratico Medio misura la media delle differenze al quadrato tra valori osservati e previsti:

\[ MSE = \frac{1}{n} \sum_{i=1}^{n} (y_i - \hat{y}_i)^2. \tag{77.10}\]

• Valori più bassi indicano previsioni più vicine ai dati osservati.
• Non tiene conto della complessità del modello, quindi può favorire modelli eccessivamente flessibili (overfitting).

Utile per valutare l’accuratezza, ma da solo non è adatto a scegliere tra modelli con diversa complessità.

77.4.2 Akaike Information Criterion (AIC)

L’AIC è un’approssimazione della divergenza \(D_{\text{KL}}\) e stima quanta informazione si perde usando un modello per descrivere i dati:

\[ AIC = -2 \sum_{i=1}^{n} \log p(y_i \mid \hat{\theta}_{\text{MLE}}) + 2k, \tag{77.11}\]

dove:

• \(\hat{\theta}_{\text{MLE}}\): stima dei parametri ottenuta massimizzando la verosimiglianza;
• \(k\): numero di parametri del modello.

Interpretazione

• Il primo termine valuta l’adattamento del modello ai dati.
• Il secondo penalizza la complessità per evitare overfitting.
• Un AIC più basso indica un miglior equilibrio tra accuratezza e semplicità.

Limiti

• Basato su assunzioni asintotiche (funziona meglio con campioni grandi).
• Usa solo stime puntuali, ignorando l’incertezza dei parametri.
• Non è pienamente coerente con l’approccio bayesiano.

77.4.3 Bayesian Information Criterion (BIC)

Il BIC valuta il compromesso tra adattamento ai dati e complessità del modello, applicando una penalizzazione più severa rispetto all’AIC — soprattutto quando il numero di osservazioni \(n\) è grande.

\[ BIC = -2 \log p(y \mid \hat{\theta}) + \log(n) \cdot k, \tag{77.12}\]

dove:

• \(p(y \mid \hat{\theta})\): massima verosimiglianza del modello (o MAP con prior piatti);
• \(n\): numero di osservazioni indipendenti;
• \(k\): numero di parametri stimati.

Interpretazione

• Il primo termine misura l’adattamento ai dati.
• Il secondo penalizza la complessità in modo crescente con \(n\) e \(k\).
• Un BIC più basso indica un compromesso migliore tra accuratezza e parsimonia.

Vantaggi

• Tende a favorire modelli più semplici quando \(n\) è elevato.
• Ha una giustificazione teorica bayesiana: in certe condizioni, approssima il log della marginal likelihood.

Limiti

• Si basa su assunzioni forti (indipendenza, modelli regolari, prior deboli).
• Può sottoselezionare modelli utili con campioni piccoli o strutture complesse.

77.4.4 Widely Applicable Information Criterion (WAIC)

Il WAIC è una versione pienamente bayesiana dell’AIC:

• utilizza tutta la distribuzione a posteriori dei parametri;
• fornisce una stima diretta della capacità predittiva del modello.

\[ WAIC = -2 \left[ \sum_{i=1}^{n} \log \left( \frac{1}{S} \sum_{s=1}^{S} p(y_i \mid \theta^{(s)}) \right) - \sum_{i=1}^{n} \mathrm{Var}_{\theta^{(s)}} \big( \log p(y_i \mid \theta^{(s)}) \big) \right], \tag{77.13}\]

dove:

• \(S\) = numero di campioni dalla distribuzione a posteriori;
• \(\theta^{(s)}\) = \(s\)-esimo campione;
• il secondo termine stima il numero effettivo di parametri basato sulla variabilità della log-verosimiglianza.

Vantaggi

• Adatto anche a modelli complessi o non regolari.
• Usa direttamente i campioni MCMC.
• Migliore dell’AIC per modelli bayesiani, perché incorpora l’incertezza dei parametri.

Nota. Il WAIC è strettamente collegato all’ELPD: è una sua stima approssimata ottenuta dalla distribuzione a posteriori, senza bisogno di eseguire la LOO-CV.

Riepilogo comparativo dei criteri di valutazione del modello
Criterio	Tipo	Penalizza la complessità	Usa stime puntuali	Basato su campioni a posteriori (es. MCMC)
MSE	Frequentista	No	Sì	No
AIC	Frequentista	Sì (modesta)	Sì	No
BIC	Frequentista/Bayesiano	Sì (forte)	Sì	No
WAIC	Bayesiano	Sì (effettiva)	No	Sì
LOO-CV	Bayesiano	Sì (empirica)	No	Sì

Riflessioni conclusive

La selezione del modello, in ottica bayesiana, 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 (e ignota) distribuzione dei dati. In termini matematici, massimizzare l’ELPD equivale a minimizzare la divergenza di Kullback–Leibler rispetto alla vera generatrice: due facce dello stesso obiettivo, rappresentare al meglio la realtà sottostante.

Poiché \(p_{\text{vera}}(y)\) è sconosciuta, l’ELPD va stimato. Le principali approssimazioni sono:

• LOO-CV (Leave-One-Out Cross-Validation): oggi lo strumento più affidabile, valuta ogni osservazione come “nuova” e stima la capacità di generalizzazione del modello.
• WAIC: alternativa completamente bayesiana, calcolata direttamente dai campioni della posteriori.
• AIC e BIC: criteri frequenstisti più rapidi ma basati su stime puntuali; utili in contesti semplici.
• MSE: misura l’accuratezza sulle osservazioni note, ma non penalizza la complessità e quindi non è adatto alla selezione del modello.

Nel confronto tra modelli, la differenza di ELPD (stimata con LOO-CV o WAIC) andrebbe interpretata insieme al relativo errore standard: una regola pratica è considerare rilevante una differenza almeno doppia rispetto all’errore standard.

In sintesi:

• la buona statistica non si limita a spiegare il passato: sa anticipare il futuro;
• la divergenza KL fornisce la misura teorica della distanza tra modello e realtà;
• l’ELPD, stimato via LOO-CV o WAIC, traduce questa misura in una valutazione pratica della capacità predittiva;
• la scelta del modello ottimale richiede un equilibrio tra accuratezza, generalizzazione e parsimonia.

Con questi strumenti possiamo individuare modelli che colgono i veri pattern nei dati, evitando di farsi ingannare dal rumore e garantendo previsioni solide anche in contesti complessi.

Informazioni sull’ambiente di sviluppo

sessionInfo()
#> R version 4.5.1 (2025-06-13)
#> Platform: aarch64-apple-darwin20
#> Running under: macOS Sequoia 15.6.1
#> 
#> 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/C/C/C/C
#> 
#> time zone: Europe/Zagreb
#> 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.0         tinytable_0.11.0     
#>  [4] patchwork_1.3.1       ggdist_3.3.3          tidybayes_3.0.7      
#>  [7] bayesplot_1.13.0      ggplot2_3.5.2         reliabilitydiag_0.2.1
#> [10] priorsense_1.1.0      posterior_1.6.1       loo_2.8.0            
#> [13] rstan_2.32.7          StanHeaders_2.32.10   brms_2.22.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.1           
#> 
#> 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.3       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.1        pkgconfig_2.0.3     
#> [13] arrayhelpers_1.1-0   fastmap_1.2.0        backports_1.5.0     
#> [16] labeling_0.4.3       cmdstanr_0.9.0       rmarkdown_2.29      
#> [19] markdown_2.0         ps_1.9.1             ragg_1.4.0          
#> [22] purrr_1.1.0          xfun_0.52            cachem_1.1.0        
#> [25] litedown_0.7         jsonlite_2.0.0       broom_1.0.9         
#> [28] parallel_4.5.1       R6_2.6.1             stringi_1.8.7       
#> [31] RColorBrewer_1.1-3   lubridate_1.9.4      estimability_1.5.1  
#> [34] knitr_1.50           zoo_1.8-14           base64enc_0.1-3     
#> [37] Matrix_1.7-3         splines_4.5.1        timechange_0.3.0    
#> [40] tidyselect_1.2.1     abind_1.4-8          yaml_2.3.10         
#> [43] codetools_0.2-20     processx_3.8.6       curl_6.4.0          
#> [46] pkgbuild_1.4.8       lattice_0.22-7       withr_3.0.2         
#> [49] bridgesampling_1.1-2 coda_0.19-4.1        evaluate_1.0.4      
#> [52] survival_3.8-3       RcppParallel_5.1.10  xml2_1.3.8          
#> [55] tensorA_0.36.2.1     checkmate_2.3.2      stats4_4.5.1        
#> [58] distributional_0.5.0 generics_0.1.4       rprojroot_2.1.0     
#> [61] commonmark_2.0.0     rstantools_2.4.0     scales_1.4.0        
#> [64] xtable_1.8-4         glue_1.8.0           emmeans_1.11.2      
#> [67] tools_4.5.1          data.table_1.17.8    mvtnorm_1.3-3       
#> [70] grid_4.5.1           QuickJSR_1.8.0       colorspace_2.1-1    
#> [73] nlme_3.1-168         cli_3.6.5            textshaping_1.0.1   
#> [76] svUnit_1.0.6         Brobdingnag_1.2-9    V8_6.0.5            
#> [79] gtable_0.3.6         sass_0.4.10          digest_0.6.37       
#> [82] TH.data_1.1-3        htmlwidgets_1.6.4    farver_2.1.2        
#> [85] memoise_2.0.1        htmltools_0.5.8.1    lifecycle_1.0.4     
#> [88] MASS_7.3-65

Bibliografia

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