Nei precedenti articoli abbiamo introdotto la statistica bayesiana, ricavato analiticamente una proporzione binomiale con priori coniugati e abbiamo descritto le basi della Markov Chain Monte Carlo tramite l’algoritmo Metropolis. In questo articolo introduciamo la modellazione di regressione nel framework bayesiano ed eseguiamo l’inferenza utilizzando la libreria MCMC di PyMC3 .
Iniziamo riassumendo l’approccio classico, o frequentista, alla regressione lineare multipla. Quindi descriviamo l’approccio di un bayesiano alla regressione lineare. Infine descriviamo brevemente il concetto di modello lineare generalizzato (GLM), necessario per comprendere la sintassi delle descrizioni dei modelli in PyMC3.
Successivamente alla descrizione di questi modelli simuliamo alcuni dati lineari con rumore e quindi usiamo PyMC3 per produrre distribuzioni a posteriori per i parametri del modello. Questa è la stessa procedura che abbiamo eseguito quando abbiamo discusso di modelli di serie temporali come ARMA e GARCH. Questo processo di “simula e adatta” non solo ci aiuta a capire il modello, ma controlla anche che lo stiamo adattando correttamente quando conosciamo i “veri” valori dei parametri.
Rivolgiamo ora la nostra attenzione all’approccio frequentista alla regressione lineare.
Regressione lineare frequentista
L’approccio frequentista, o classico, alla regressione lineare multipla assume 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}\)
Dove \(\beta^T\) è la trasposta del vettore coefficiente \(\beta\) e \(\epsilon \sim \mathcal{N}(0,\sigma^2)\) è l’errore della misurazione, distribuito normalmente con media zero e deviazione standard \(\sigma\).
Cioè, il nostro modello \(f(\mathbf{X})\) è lineare per i predittori,\(\mathbf{X}\), associato a qualche errore di misurazione.
Se abbiamo una serie di dati di allenamento \((x_1, y_1), \ldots, (x_N, y_N)\) l’obiettivo è quindi stimare i coefficienti \(\beta\), che forniscono il miglior adattamento lineare ai dati. Geometricamente, questo significa che dobbiamo trovare l’orientamento dell’iperpiano che meglio caratterizza linearmente i dati.
“Migliore” in questo caso significa ridurre al minimo una qualche forma di funzione di errore. Il metodo più popolare per farlo è tramite i minimi quadrati ordinari (OLS). Se consideriamo la somma residua dei quadrati (RSS), che è la somma delle differenze al quadrato tra gli output e le stime di regressione lineare, come segue:
\(\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}\)
Quindi l’obiettivo di OLS è ridurre al minimo l’RSS, tramite l’adeguamento dei coefficienti \(\beta\). La stima di massima verosimiglianza di \(\beta\), che minimizza l’RSS, è data da (vedi Hastie et al [2] per i dettagli):
\(\begin{eqnarray}\hat{\beta} = (\mathbf{X}^T \mathbf{X})^{-1} \mathbf{X}^T \mathbf{y}\end{eqnarray}\)
Per ricavare la previsione successiva \(y_{N+1}\), considerando i nuovi dati \(x_{N+1}\), è sufficiente moltiplicare i componenti di \(x_{N+1}\) con i relativi coefficienti \(\beta\) e ottenere \(y_{N+1}\).
E’ importante evidenziare che \(\hat{\beta}\) è una stima puntuale, cioè è un valore singolo di \(\mathbb{R}^{p+1}\). Nella formulazione bayesiana questa interpretazione differisce sostanzialmente.
Regressione lineare bayesiana
In un framework bayesiano, la regressione lineare è espressa in modo probabilistico. Cioè, riformuliamo il modello di regressione lineare di cui sopra per utilizzare le distribuzioni di probabilità. La sintassi di una regressione lineare in un framework bayesiano è simile alla seguente:
\(\begin{eqnarray}\mathbf{y} \sim \mathcal{N} \left(\beta^T \mathbf{X}, \sigma^2 \mathbf{I} \right)\end{eqnarray}\)
In parole, i nostri datapoint di risposta \(\mathbf{y}\) sono campionati da una distribuzione normale multivariata che ha una media uguale al prodotto tra i coefficienti \(\beta\) e predittori \(\mathbf{X}\), e una varianza di \(\sigma^2\). Inoltre \(\mathbf{I}\) si riferisce alla matrice identità, necessaria affinché la distribuzione sia multivariata.
Questa è una formulazione molto diversa dall’approccio frequentista. Nell’impostazione frequentista non si fa menzione delle distribuzioni di probabilità per qualcosa di diverso dall’errore di misurazione. Nella formulazione bayesiana l’intero problema è riformulato in modo tale che i valori \(y_i\) sono campioni di una distribuzione normale.
Una domanda comune in questa fase è “Qual è il vantaggio di farlo?”. Cosa otteniamo da questa riformulazione? Ci sono due ragioni principali per farlo ( Wiecki [6]):
- Distribuzioni a priori: se abbiamo una conoscenza preliminare dei parametri, possiamo scegliere distribuzioni a priori che riflettano questi parametri. In caso contrario, possiamo comunque scegliere priori non informativi.
- Distribuzioni a posteriori: come menzionato in precedenza, il valore MLE frequentista per i nostri coefficienti di regressione, \(\hat{\beta}\), è solo una stima puntuale. Nella formulazione bayesiana riceviamo un’intera distribuzione di probabilità che caratterizza la nostra incertezza per i differenti coefficienti \(\beta\). Il vantaggio immediato è, dopo aver preso in considerazione tutti i dati, la possibilità di quantificare l’incertezza nei parametri \(\hat{\beta}\) tramite la varianza della distribuzione a posteriori. Una varianza maggiore indica una maggiore incertezza.
Sebbene la formula dell’approccio bayesiano possa sembrare concisa, in realtà non ci dà molti indizi su come specificare un modello e campionarlo usando Markov Chain Monte Carlo. Nei prossimi paragrafi vediamo come usare PyMC3 per formulare e utilizzare un modello di regressione lineare bayesiana.
Regressione lineare bayesiana con PyMC3
In questo paragrafo vediamo un approccio consolidato con esempi statistici, cioè simuliamo alcuni dati con proprietà che già conosciamo, e quindi adattiamo un modello per verificare la presenza di queste proprietà. Abbiamo usato questa tecnica molte negli articoli precedenti, principalmente negli articoli sull’analisi delle serie temporali.
Nonostante seguire una tale procedura possa sembrare complicato, ci sono in realtà due principali vantaggi. Il primo vantaggio è la possibilità di capire esattamente come funziona il fitting del modello. Per farlo, dobbiamo prima capirlo, quindi ci aiuta a intuire come funziona il modello. Il secondo vantaggio è la possibilità di vedere come si comporta il modello (cioè i valori e l’incertezza che restituisce) in una situazione in cui conosciamo effettivamente i veri valori che stiamo cercano di stimare.
In questo approccio con Python usiamo le librerie numpy e pandas per simulare i dati, usiamo seaborn per visualizzare i grafici e infine il modulo Generalized Linear Models (GLM) di PyMC3 per effettuare una regressione lineare bayesiana e campionarla, sul nostro set di dati simulati.
La seguente 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 la regressione lineare bayesiana, è opportuno introdurre brevemente il concetto di modello lineare generalizzato (GLM), poiché è utile per formulare il nostro modello in PyMC3.
Un modello lineare generalizzato è un meccanismo flessibile per estendere la regressione lineare ordinaria a forme di regressione più generali, inclusa la regressione logistica (classificazione) e la regressione di Poisson (usata per i dati discreti), nonché la regressione lineare stessa.
I GLM permetto alle variabili endogene (le variabili di risposta) di avere distribuzioni di errore diverse dalla distribuzione normale (vedi \(\epsilon\) sopra, nel paragrafo della regressione frequentista). Nel modello lineare generalizzato ciascun valore dalla variabile dipendente \(\mathbf{y}\) si assume venga generato da una particolare variabile casuale della famiglia di distribuzioni esponenziali. Questa famiglia di distribuzioni comprende molte distribuzioni comuni tra cui normale, gamma, beta, chi quadrato, Bernoulli, Poisson e altri.
La media \(\mathbf{\mu}\) di questa distribuzione, dipende da \(\mathbf{X}\) tramite la seguente relazione:
\(\begin{eqnarray}\mathbb{E}(\mathbf{y}) = \mu = g^{-1}(\mathbf{X}\beta)\end{eqnarray}\)
Dove \(g\) è la funzione di collegamento. In questo ambito, la varianza è tipicamente 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’impostazione frequentista, come con la normale regressione lineare sopra, i coefficienti ignoti \(\)\beta\(\) sono stimati attraverso un approccio di massima verosimiglianza.
Non descriviamo in modo approfondito i GLM perchè non sono oggetto di questo articolo. Li abbiamo introdotti perché dobbiamo usare il modulo glm di PyMC3, che è stato scritto da Thomas Wiecki e altri[6], 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, dobbiamo simulare alcuni dati lineari rumorosi. Il seguente codice esegue queste operazioni (è una modifica ed estensione del codice presente nel post di Jonathan Sedar [5]):
Per prima cosa creiamo i valori x, che sono valori in virgola mobile equispaziati tra 0 e 1. Per fare questo usiamo np.linspace(). Quindi usiamo un modello lineare per simulare i valori y con un’intercetta di 1 e una pendenza di 2. Per creare il set di dati rumoroso aggiungiamo un valore, epsilon, ai valori di y. Questi valori vengono generati utilizzando il metodo RandomState.normal() di numpy, che estrae campioni casuali da una distribuzione normale. Impostiamo la media o “loc” uguale a 0 e la deviazione standard o “scala” uguale a 0,5.
import bambi as bmb
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc as pm
import seaborn as sns
sns.set(style="darkgrid", palette="muted")
def simulate_linear_data(
start, stop, N, beta_0, beta_1, eps_mean, eps_sigma_sq
):
"""
Simula un dataset random dataset usando il noisy
linear process.
Parameters
----------
N: `int`
Numero dei punti data da simulare
beta_0: `float`
Intercetta
beta_1: `float`
Pendeza del preditore univariato, X
Returns
-------
df: `pd.DataFrame`
Un DataFrame che contiene i valori x e y
"""
# Crea un DataFrame pandas con la colonna 'x' che contiene
# N valori campionati uniformemente tra 0.0 e 1.0
df = pd.DataFrame(
{"x":
np.linspace(start, stop, num=N)
}
)
# Utilizza un modello lineare (y ~ beta_0 + beta_1*x + epsilon)
# per generare una colonna "y" di risposte basate su "x"
df["y"] = beta_0 + beta_1*df["x"] + np.random.RandomState(s).normal(
eps_mean, eps_sigma_sq, N
)
return df
def plot_simulated_data(df):
"""
Tracciare i dati simulati con sns.lmplot().
Parameters
----------
df: `pd.DataFrame`
Un DataFrame che contiene i valori x e y
"""
# Visualizza i dati e la regressione lineare
# frequentista utilizzando il pacchetto seaborn
sns.lmplot(x="x", y="y", data=df, height=10)
plt.xlim(0.0, 1.0)
plt.show()
if __name__ == "__main__":
# Questi sono i nostri parametri "Reali"
beta_0 = 1.0 # Intercept
beta_1 = 2.0 # Slope
# Simula 100 punti dati tra 0 e 1, con una varianza di 0,5
start = 0
stop = 1
N = 100
eps_mean = 0.0
eps_sigma_sq = 0.5
# Random Seed
s = 42
# Simula i dati "lineari" utilizzando i parametri di cui sopra
df = simulate_linear_data(
start, stop, N, beta_0, beta_1, eps_mean, eps_sigma_sq
)
# Visualizza il grafico dei dati simulati
plot_simulated_data(df)
Abbiamo simulato 100 punti dati, con un’intercetta \(\beta_0=1\) e una pendenza \(\beta_1=2\). I valori epsilon sono normalmente distribuiti con media zero e varianza \(\sigma^2=\frac{1}{2}\). I dati sono stati tracciati utilizzando la funzione sns.lmplot. Inoltre, il metodo utilizza un approccio MLE frequentista per adattare una linea di regressione lineare ai dati.
Dopo aver eseguito la simulazione, possiamo adattare una regressione lineare bayesiana ai dati. È qui che il modulo glm entra in gioco. Insieme alla libreria Bambi, come descritto nel tutorial PyMC, il modulo utilizza una sintassi di specifica del modello. La libreria bambi considera un identificatore della formula del modello lineare da cui crea una matrice di progettazione, quindi aggiunge variabili casuali per ciascuno dei coefficienti e un’opportuna verosimiglianza al modello. Se bambi non è installato, si può installarlo con il comando <code>pip install bambi</code>.
Infine, usiamo il No-U-Turn Sampler (NUTS) per eseguire l’inferenza effettiva e quindi tracciare il grafico del modello. Consideriamo 5000 campioni dalla distribuzione a posteriori e usiamo i primi 500 campioni per mettere a punto il modello, quindi scartarli come “burn in”. PyMC consente di campionare più catene in base al numero di CPU presenti nella macchina. Per impostazione predefinita è impostato sul numero di core nel sistema, in questo caso lo abbiamo impostato esplicitamente ad 1.
import bambi as bmb
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc as pm
import seaborn as sns
...
def glm_mcmc_inference(df, iterations=5000):
"""
Calcola la traccia della catena di Markov Monte Carlo della
regressione lineare bayesiana del modello lineare
generalizzato sui dati forniti.
Parameters
----------
df: `pd.DataFrame`
DataFrame che contiene i dati
iterations: `int`
Numero di iteazioni per calcolare il MCMC
"""
# Crea il glm utilizzando la sintassi del modello Bambi
model = bmb.Model("y ~ x", df)
# Addestrare il modello utilizzando un NUTS (No-U-Turn Sampler)
trace = model.fit(
draws=5000,
tune=500,
discard_tuned_samples=True,
chains=1,
progressbar=True)
return trace
def plot_glm_model(trace):
"""
Visualizza il grafico generato dall'addestramento del modello.
Parameters
----------
trace: `tracepymc.backends.base.MultiTrace`
Un oggetto MultiTrace o ArviZ InferenceData che contiene i campioni.
"""
pm.plot_trace(trace)
plt.tight_layout()
plt.show()
if __name__ == "__main__":
# Questi sono i nostri parametri "Reali"
beta_0 = 1.0 # Intercept
beta_1 = 2.0 # Slope
# Simula 100 punti dati tra 0 e 1, con una varianza di 0,5
start = 0
stop = 1
N = 100
eps_mean = 0.0
eps_sigma_sq = 0.5
# Random Seed
s = 42
# Simula i dati "lineari" utilizzando i parametri di cui sopra
df = simulate_linear_data(
start, stop, N, beta_0, beta_1, eps_mean, eps_sigma_sq
)
# Visualizza il grafico dei dati simulati
plot_simulated_data(df)
# Addestra il GLM
trace = glm_mcmc_inference(df, iterations=5000)
# Visualizza il GLM
plot_glm_model(trace)
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [Intercept, x, y_sigma]
[-----------------]Sampling 1 chains for 500 tune and 5_000 draw iterations (500 + 5_000 draws total) took 12 seconds.
Nel precedente articolo sull’algoritmo Metropolis MCMC abbiamo descritto le basi dei traceplot. Ricordiamo che i modelli bayesiani forniscono una completa distribuzione di probabilità a posteriori per ciascuno dei parametri del modello, al contrario di una stima frequentista puntiforme.
Sul lato sinistro del pannello possiamo vedere le distribuzioni marginali per ogni parametro di interesse. Si noti che la distribuzione dell’intercetta \(\beta_0\) ha la sua stima modale/massima a posteriori quasi a 1, vicino al valore reale del parametro \(\beta_0=1\). La stima per il parametro \(\beta_1\) della pendenza(x) ha una moda a circa 2.1, vicino al valore reale del parametro \(\beta_1=2\). Il parametro di errore y_sigma associato al rumore di misurazione del modello ha una moda di circa 0.46, che è leggermente fuori rispetto al valore reale \(\epsilon=0.5\).
In tutti i casi c’è una ragionevole varianza associata a ciascun marginale a posteriori, che significa un certo grado di incertezza in ciascuno dei valori. Se possiamo simulare più dati ed eseguire più campioni, probabilmente questa varianza diminuirebbe.
Il punto chiave da sottolineare è l’assenza di una stima puntuale per una retta di regressione, cioè “una retta di migliore adattamento”, come avviene nel caso frequentista. In questo caso otteniamo una distribuzione di probabili linee di regressione.
Possiamo tracciare queste linee insieme alla vera linea di regressione e ai dati simulati che abbiamo creato in precedenza. Iniziamo creando la nostra figura in cui il nostro asse x è un numero in virgola mobile spaziato linearmente tra 0 e 1 con 100 intervalli. Quindi tracciamo i nostri dati simulati usando la funzione plt.scatter di Matplotlib e il DataFrame originale, df.
Possiamo accedere ai dati numerici per le distribuzioni posteriori dall’oggetto trace creato dal modello. I dati di inferenza della traccia hanno i seguenti gruppi: posterior, log_likelihood, sample_stats, observed_data. Il gruppo posterior contiene i valori dei 5000 campioni delle distribuzioni a posteriori dell’intercetta e della pendenza(x). Questi sono array nella forma(1, 5000). Il primo valore rappresenta il numero di catene, il secondo è il numero di estrazioni/campioni. Questo set di dati presenta anche un metodo to_numpy() che converte i dati in un array numpy. Quindi, per inserire i valori delle intercette in un array da cui possiamo campionare e tracciare, dobbiamo usare il comando trace.posterior.Intercept.to_numpy()[0]. Per ottenere i dati sulle pendenze eseguiamo un comando simile, sostituendo ‘Intercepts’ con ‘x’, cioè trace.posterior.x.to_numpy()[0]. Ora abbiamo due matrici con i valori numerici delle distribuzioni a posteriori per le nostre pendenze e intercette.
Per tracciare un campione di 100 linee, dobbiamo prelevare 100 di questi valori dai nostri array. Per prima cosa creiamo un elenco di indici di esempio usando un elenco di numeri interi casuali con il metodo random.randint() di numpy. Quindi eseguiamo il ciclo attraverso l’array con il valore dell’indice e creiamo una linea utilizzando l’intercetta e la pendenza, quindi tracciando la linea sul nostro grafico. Infine tracciamo la nostra vera linea di regressione usando le variabili beta_0 e beta_1 dei dati simulati. Il frammento di codice seguente produce questo grafico:
...
def plot_regression_lines(trace, df, N):
"""
Visualizza i dati simulati con le linee di regressione reali e stimate.
Parameters
----------
trace: `tracepymc.backends.base.MultiTrace`
Un oggetto MultiTrace o ArviZ InferenceData che contiene i campioni.
df: `pd.DataFrame`
DataFrame che contiene i dati
N: `int`
Numero di punti data da simulare
"""
fig, ax = plt.subplots(figsize=(7, 7))
# definisce i tick dell'asse x
x = np.linspace(0, 1, N)
# visualizza le osservazioni di dati simulati
ax.scatter(df['x'], df['y'])
# estrae la pendeza e l'intercetta dalla traccia PyMC
intercepts = trace.posterior.Intercept.to_numpy()[0]
slopes = trace.posterior.x.to_numpy()[0]
# visualizza il grafico di 100 campioni random dalla pendenza e l'intercetta
sample_indexes = np.random.randint(len(intercepts), size=100)
for i in sample_indexes:
y_line = intercepts[i] + slopes[i] * x
ax.plot(x, y_line, c='black', alpha=0.07)
# visualizza la regressione lineare
y = beta_0 + beta_1*x
ax.plot(x, y, label="True Regression Line", lw=3., c="green")
ax.legend(loc=0)
ax.set_xlim(0.0, 1.0)
ax.set_ylim(0.0, 4.0)
plt.show()
if __name__ == "__main__":
# Questi sono i nostri parametri "Reali"
beta_0 = 1.0 # Intercept
beta_1 = 2.0 # Slope
# Simula 100 punti dati tra 0 e 1, con una varianza di 0,5
start = 0
stop = 1
N = 100
eps_mean = 0.0
eps_sigma_sq = 0.5
# Random Seed
s = 42
# Simula i dati "lineari" utilizzando i parametri di cui sopra
df = simulate_linear_data(
start, stop, N, beta_0, beta_1, eps_mean, eps_sigma_sq
)
# Visualizza il grafico dei dati simulati
plot_simulated_data(df)
# Addestra il GLM
trace = glm_mcmc_inference(df, iterations=5000)
# Visualizza il GLM
plot_glm_model(trace)
plot_regression_lines(trace, df, N)
Prossimi passi
Nell’articolo precedente abbiamo esaminato un semplice metodo MCMC chiamato algoritmo Metropolis. In quell’articolo abbiamo espresso la volontà di descrivere come funziona i vari “algoritmi” di MCMC “sotto il cofano”. Nei prossimi articoli prendiamo in considerazione gli algoritmi di Gibbs Sampler, Hamiltonian Sampler e No-U-Turn Sampler, tutti utilizzati nei principali pacchetti software bayesiani.
Alla fine vediamo l’approccio di una robusta regressione e modelli lineari gerarchici, una potente tecnica di modellazione resa possibile dalle rapide implementazioni MCMC. Da un punto di vista della finanza quantitativa descriviamo anche un modello di volatilità stocastica che utilizza PyMC e vediamo come possiamo utilizzare questo modello per formare algoritmi di trading.
Nota bibliografica
Un’introduzione alla regressione lineare frequentista può essere trovata in James et al (2013)[4] . Una panoramica più tecnica, compresi i metodi di selezione dei sottoinsiemi, può essere trovata in Hastie et al (2009)[2]. Gelman et al (2013)[1] discutono in modo approfondito i modelli lineari bayesiani a un livello ragionevolmente tecnico.
Questo articolo è fortemente influenzato dai precedenti post sul blog di Thomas Wiecki[6], inclusa la sua discussione sui GLM bayesiani, nonché dai post di Jonathan Sedar[5].
Riferimenti
- Gelman, A. et al (2013) Bayesian Data Analysis, 3rd Edition, Chapman and Hall/CRC
- Hastie, T., Tibshirani, R., Friedman, J. (2009) The Elements of Statistical Learning, Springer
- Hoffman, M.D., and Gelman, A. (2011) “The No-U-Turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carlo, arXiv:1111.4246 [stat.CO]
- James, G., Witten, D., Hastie, T., Tibshirani, R. (2013) An Introduction to Statistical Learning, Springer
- Sedar, J. (2016) Bayesian Inference with PyMC3 – Part 1
- Wiecki, T. (2015) The Inference Button: Bayesian GLMs made easy with PyMC3, https://twiecki.github.io/blog/2013/08/12/bayesian-glms-1/
