Bayesian Linear Models: A compact general set of results

1. Problema e Contesto

Il documento affronta la necessità di una presentazione compatta, generale e accessibile delle derivazioni algebriche fondamentali per i Modelli Lineari Bayesiani (BLM) con prior coniugati Normal-Gamma.
Sebbene i risultati teorici siano disponibili in numerosi testi accademici, spesso sono dispersi o presentati in modo non uniforme. Il problema specifico risolto è la mancanza di un riferimento unificato che includa:

• La derivazione completa della distribuzione a posteriori.
• La gestione di osservazioni correlate (covarianza non diagonale), cruciale per serie temporali e statistica spaziale.
• La predizione, il confronto tra modelli (selezione delle variabili) e l'implementazione numerica efficiente.
• L'estensione a casi particolari come regressione, processi Gaussiani (GP) e Modelli Lineari Dinamici Bayesiani (DLM).

2. Metodologia

L'autore utilizza un approccio algebrico rigoroso basato sulla distribuzione Normal-Gamma coniugata.

Il Modello

Il modello lineare generale è definito come:
$$Y = X\theta + \epsilon; \quad \epsilon \sim N_n(0, \lambda^{-1}\Sigma)$$
Dove:

• $Y$ è il vettore delle risposte ($n \times 1$).
• $X$ è la matrice di design ($n \times p$).
• $\theta$ è il vettore dei parametri sconosciuti ($p \times 1$).
• $\Sigma$ è la matrice di covarianza nota (con $A = \Sigma^{-1}$ come matrice di precisione).
• $\lambda$ è il parametro di precisione (incognito), che scala la varianza.

Prior Coniugata

La prior è definita come una distribuzione Normal-Gamma:

• Condizionata a $\lambda$, $\theta \sim N_p(\theta_0, \lambda^{-1}A_0^{-1})$.
• $\lambda \sim Ga(\alpha_0, \beta_0)$ (Gamma con parametri di forma e tasso).

Derivazione della Posterior

L'autore espande la forma quadratica nell'esponente della funzione di verosimiglianza combinata con la prior. Utilizzando identità matriciali chiave (Lemma di Sherman-Woodbury-Morrison e Lemma del determinante matriciale), dimostra che la distribuzione a posteriori congiunta $\pi(\theta, \lambda | y)$ rimane nella famiglia Normal-Gamma con parametri aggiornati:

• Precisione dei parametri: $A_n = A_0 + X'AX$
• Media a posteriori: $\theta_n = A_n^{-1}(A_0\theta_0 + X'Ay)$
• Parametri Gamma: $\alpha_n = \alpha_0 + n/2$ e $\beta_n = \beta_0 + \frac{1}{2}q(y|X, \theta_n, A_0, \theta_0)$.

Il documento offre tre espressioni equivalenti per $\beta_n$, evidenziando quella che separa la somma degli errori quadratici pesati dalla penalità sulla discrepanza tra media a priori e a posteriori:
$$\beta_n = \beta_0 + \frac{1}{2} \left[ (y - X\theta_n)'A(y - X\theta_n) + (\theta_0 - \theta_n)'A_0(\theta_0 - \theta_n) \right]$$

Distribuzioni Predittive e Marginali

• Predizione: La distribuzione predittiva per nuove osservazioni $Z$ è una distribuzione t multivariata (se $\lambda$ è incognito) o Normale (se $\lambda$ è noto).
• Marginali: La distribuzione marginale di $\lambda$ è Gamma, mentre quella dei singoli parametri $\theta_j$ è una distribuzione t di Student.

3. Contributi Chiave

1. Generalità delle Correlazioni: A differenza di molte trattazioni standard che assumono dati indipendenti ($\Sigma = I$), questo lavoro tratta esplicitamente il caso in cui $\Sigma$ è una matrice di covarianza generica (simmetrica e definita positiva). Questo rende i risultati immediatamente applicabili a serie temporali e geostatistica.
2. Formule Compatte per il Confronto di Modelli: Viene derivata una forma compatta per la "evidenza del modello" (normalizing constant), essenziale per il calcolo delle probabilità a posteriori dei modelli e la selezione delle variabili, evitando calcoli ridondanti.
3. Ottimizzazione Numerica: Il documento fornisce consigli pratici per l'implementazione efficiente in Python, suggerendo l'uso della decomposizione di Cholesky e identità matriciali per evitare il calcolo diretto di determinanti di matrici grandi ($n \times n$), riducendo il carico computazionale da $O(n^3)$ a operazioni su matrici più piccole ($p \times p$).
4. Unificazione di Casi Particolari: Dimostra come Regressione Lineare, Processi Gaussiani e Modelli Dinamici (DLM) siano casi specifici di questo quadro generale, fornendo le formule di aggiornamento specifiche per ciascuno.
5. Risorsa Software: Viene fornita un'implementazione Python completa (classe BLM) disponibile su GitHub, che traduce la teoria in codice pratico.

4. Risultati ed Esempi

L'autore illustra l'applicazione dei risultati attraverso:

• Regressione Polinomiale: Simulazione di dati con rumore correlato e non, mostrando l'adattamento MAP (Maximum A Posteriori) e le distribuzioni marginali a posteriori per diversi gradi polinomiali ($p=1..6$). Viene mostrato come calcolare la probabilità a posteriori di ciascun modello.
• Processi Gaussiani (GP): Applicazione alla predizione spaziale. Viene mostrato come il modello funzioni come un interpolatore (varianza zero nei punti osservati) o come smussatore (se si include un "nugget" di rumore), utilizzando la distribuzione predittiva t.
• Confronto Visivo: Le figure includono l'adattamento del modello, le distribuzioni marginali dei parametri e le probabilità dei modelli, confermando la coerenza tra la teoria e l'implementazione numerica.

5. Significato

Questo lavoro è significativo perché funge da ponte tra la teoria bayesiana astratta e l'implementazione pratica.

• Accessibilità: Fornisce un "manuale di riferimento" compatto che elimina la necessità di consultare molteplici testi per derivare le stesse formule algebriche.
• Versatilità: La capacità di gestire osservazioni correlate lo rende uno strumento potente per statistici che lavorano su dati spaziali o temporali, dove l'assunzione di indipendenza è spesso violata.
• Efficienza: I suggerimenti sull'ottimizzazione numerica sono cruciali per l'analisi di dataset di grandi dimensioni ($n \gg p$), rendendo l'approccio bayesiano computazionalmente fattibile.
• Riproducibilità: La disponibilità del codice Python garantisce che i risultati teorici siano immediatamente verificabili e riutilizzabili dalla comunità di ricerca.

In sintesi, il documento è una risorsa tecnica essenziale per chiunque desideri implementare o comprendere a fondo l'inferenza bayesiana in modelli lineari, offrendo un equilibrio raro tra rigore matematico, generalità statistica e utilità pratica computazionale.