Modelli di regressione lineare bayesiana con PyMC3

Nelle precedenti lezioni di questo corso abbiamo introdotto la statistica bayesiana, ricavato analiticamente una proporzione binomiale con priori coniugati e descritto le basi della Markov Chain Monte Carlo tramite l’algoritmo Metropolis. In questa lezione introduciamo i modelli di regressione lineare bayesiana ed eseguiamo l’inferenza con la libreria MCMC di PyMC3.

Cominciamo riassumendo l’approccio classico, o frequentista, alla regressione lineare multipla. Proseguiamo con una descrizione dell’approccio bayesiano alla regressione lineare. Infine, introduciamo brevemente il concetto di modello lineare generalizzato (GLM), utile per comprendere la sintassi dei modelli in PyMC3.

Dopo questa panoramica, simuliamo alcuni dati lineari con rumore e utilizziamo PyMC3 per ottenere le distribuzioni a posteriori dei parametri del modello. Usiamo lo stesso procedimento applicato nei modelli di serie temporali come ARMA e GARCH. Questo approccio “simula e adatta” ci consente sia di comprendere meglio il modello, sia di verificare l’accuratezza dell’adattamento confrontandolo con i valori reali dei parametri.

Ora concentriamo l’attenzione sull’approccio frequentista alla regressione lineare.

Regressione lineare frequentista

L’approccio frequentista, o classico, alla regressione lineare multipla considera un modello della forma (Hastie et al [2]):

\(\begin{eqnarray}f \left( \mathbf{X} \right) = \beta_0 + \sum_{j=1}^p \mathbf{X}_j \beta_j + \epsilon = \beta^T \mathbf{X} + \epsilon\end{eqnarray}\)

Qui \(\beta^T\) rappresenta la trasposta del vettore coefficiente \(\beta\) e \(\epsilon \sim \mathcal{N}(0,\sigma^2)\) indica l’errore di misurazione, distribuito normalmente con media zero e deviazione standard \(\sigma\).

In altre parole, il modello \(f(\mathbf{X})\) risulta lineare rispetto ai predittori \(\mathbf{X}\), con l’aggiunta di un errore di misurazione.

Se disponiamo di una serie di dati di allenamento \((x_1, y_1), \ldots, (x_N, y_N)\), puntiamo a stimare i coefficienti \(\beta\) per ottenere il miglior adattamento lineare ai dati. Dal punto di vista geometrico, cerchiamo l’orientamento dell’iperpiano che caratterizza meglio i dati in modo lineare.

In questo contesto, “migliore” significa minimizzare una funzione di errore. Il metodo più diffuso per raggiungere questo obiettivo è quello dei minimi quadrati ordinari (OLS). Consideriamo la somma residua dei quadrati (RSS), cioè la somma dei quadrati delle differenze tra output osservati e stime di regressione lineare:

\(\begin{eqnarray}\text{RSS}(\beta) &=& \sum_{i=1}^{N} (y_i – f(x_i))^2 \\
&=& \sum_{i=1}^{N} (y_i – \beta^T x_i)^2\end{eqnarray}\)

Il metodo OLS mira quindi a minimizzare l’RSS, regolando i coefficienti \(\beta\). La stima di massima verosimiglianza di \(\beta\), che riduce l’RSS al minimo, risulta:

\(\begin{eqnarray}\hat{\beta} = (\mathbf{X}^T \mathbf{X})^{-1} \mathbf{X}^T \mathbf{y}\end{eqnarray}\)

Per prevedere il valore successivo \(y_{N+1}\) in base ai nuovi dati \(x_{N+1}\), basta moltiplicare i componenti di \(x_{N+1}\) per i relativi coefficienti \(\beta\) e ottenere \(y_{N+1}\).

Vale la pena sottolineare che \(\hat{\beta}\) rappresenta una stima puntuale, ovvero un singolo valore di \(\mathbb{R}^{p+1}\). Nella formulazione bayesiana, questa interpretazione cambia radicalmente.

Modelli di regressione lineare bayesiana

I modelli di regressione lineare bayesiana esprimiamo la regressione lineare in modo probabilistico. Riformuliamo quindi il modello di regressione lineare introducendo le distribuzioni di probabilità. La sintassi assume la seguente forma:

\(\begin{eqnarray}\mathbf{y} \sim \mathcal{N} \left(\beta^T \mathbf{X}, \sigma^2 \mathbf{I} \right)\end{eqnarray}\)

In parole semplici, i valori di risposta \(\mathbf{y}\) derivano da una distribuzione normale multivariata con media uguale al prodotto tra coefficienti \(\beta\) e predittori \(\mathbf{X}\), e varianza pari a \(\sigma^2\). La matrice identità \(\mathbf{I}\) garantisce che la distribuzione sia multivariata.

Questa formulazione si discosta nettamente dall’approccio frequentista. Quest’ultimo non considera distribuzioni di probabilità se non per l’errore di misurazione, mentre il modello bayesiano descrive ogni valore \(y_i\) come campione da una distribuzione normale.

Una domanda naturale che sorge a questo punto è: “Qual è il vantaggio di questa riformulazione?” Possiamo individuare due principali motivi per adottarla (Wiecki [6]):

• Distribuzioni a priori: Se disponiamo di conoscenze preliminari sui parametri, possiamo impostare distribuzioni a priori che le riflettano. In alternativa, possiamo scegliere priori non informativi.
• Distribuzioni a posteriori: Come già osservato, il valore MLE frequentista \(\hat{\beta}\) rappresenta solo una stima puntuale. Invece, il framework bayesiano fornisce un’intera distribuzione di probabilità che caratterizza l’incertezza sui coefficienti \(\beta\). Dopo aver considerato tutti i dati, possiamo quantificare l’incertezza sui parametri \(\hat{\beta}\) attraverso la varianza della distribuzione a posteriori: una varianza più alta implica maggiore incertezza.

Nonostante la concisione della formula bayesiana, essa non offre indicazioni pratiche su come specificare e campionare un modello usando Markov Chain Monte Carlo. Nei prossimi paragrafi vediamo come utilizzare PyMC3 per formulare e sfruttare un modello di regressione lineare bayesiana.

Regressione lineare bayesiana con PyMC3

In questo paragrafo seguiamo un approccio consolidato con esempi statistici: simuliamo alcuni dati con proprietà già note e adattiamo i modelli di regressione lineare bayesiana per verificare la presenza di queste proprietà. Abbiamo già utilizzato questa tecnica in precedenti lezioni, soprattutto nel corso sull’analisi delle serie temporali.

Seguire una tale procedura potrebbe sembrare complicato, ma offre in realtà due principali vantaggi. Il primo consiste nella possibilità di comprendere esattamente come funziona il fitting del modello. Per farlo, dobbiamo prima capirlo, il che ci aiuta a intuire il funzionamento del modello. Il secondo vantaggio consiste nella possibilità di osservare come si comporta il modello (cioè i valori e l’incertezza che restituisce) in una situazione in cui conosciamo i veri valori che stiamo stimando.

In questo approccio con Python usiamo le librerie numpy e pandas per simulare i dati, usiamo seaborn per visualizzare i grafici e infine utilizziamo il modulo Generalized Linear Models (GLM) di PyMC3 per effettuare una regressione lineare bayesiana e campionarla sul nostro set di dati simulati.

Questa analisi si basa principalmente su una raccolta di post di blog scritti da Thomas Wiecki [6] e Jonathan Sedar [5], insieme a basi bayesiane più teoriche di Gelman et al [1].

Cosa sono i modelli lineari generalizzati?

Prima di iniziare a descrivere i modelli di regressione lineare bayesiana, introduciamo brevemente il concetto di modello lineare generalizzato (GLM), utile per formulare il nostro modello in PyMC3.

Un modello lineare generalizzato estende in modo flessibile la regressione lineare ordinaria a forme più generali, come la regressione logistica (classificazione), la regressione di Poisson (per dati discreti) e la stessa regressione lineare.

I GLM permettono alle variabili endogene (le variabili di risposta) di seguire distribuzioni di errore diverse dalla normale (vedi \(\epsilon\) sopra, nel paragrafo della regressione frequentista). Nel modello lineare generalizzato ogni valore della variabile dipendente \(\mathbf{y}\) deriva da una variabile casuale appartenente alla famiglia delle distribuzioni esponenziali. Questa famiglia comprende molte distribuzioni comuni tra cui normale, gamma, beta, chi quadrato, Bernoulli, Poisson e altre.

La media \(\mathbf{\mu}\) di questa distribuzione dipende da \(\mathbf{X}\) secondo la relazione:

\(\begin{eqnarray}\mathbb{E}(\mathbf{y}) = \mu = g^{-1}(\mathbf{X}\beta)\end{eqnarray}\)

dove \(g\) rappresenta la funzione di collegamento. In questo ambito, la varianza si presenta tipicamente come una funzione \(V\) della media:

\(\begin{eqnarray}\text{Var}(\mathbf{y}) = V(\mathbb{E}(\mathbf{y})) = V(g^{-1}(\mathbf{X}\beta))\end{eqnarray}\)

Nell’approccio frequentista, come nella regressione lineare sopra, stimiamo i coefficienti ignoti \(\beta\) attraverso l’approccio di massima verosimiglianza.

Non approfondiamo i GLM perché non costituiscono l’oggetto di questa lezione. Li abbiamo introdotti perché dobbiamo utilizzare il modulo glm di PyMC3, che Thomas Wiecki e altri [6] hanno sviluppato, per specificare facilmente la nostra regressione lineare bayesiana.

Simulazione dei dati e adattamento del modello con PyMC3

Prima di utilizzare PyMC3 per specificare e campionare un modello bayesiano, simuliamo alcuni dati lineari rumorosi. Il seguente codice svolge queste operazioni (modifica ed estende il codice presente nel post di Jonathan Sedar [5]):

Per prima cosa creiamo i valori x, equispaziati in virgola mobile tra 0 e 1, utilizzando np.linspace(). Usiamo poi un modello lineare per simulare i valori y, impostando un’intercetta di 1 e una pendenza di 2. Per creare il set di dati rumoroso aggiungiamo ai valori di y un errore casuale, epsilon, generato con il metodo RandomState.normal() di numpy, che estrae campioni da una distribuzione normale. Impostiamo la media (“loc”) uguale a 0 e la deviazione standard (“scale”) uguale a 0,5.