Modelli Bayesiani Regressione Lineare PyMC3

Modelli di regressione lineare bayesiana con PyMC3

Sommario

Nei precedenti articoli abbiamo introdotto la statistica bayesianaricavato 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]):
				
					import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
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
):
    """
    Simulate a random dataset using a noisy
    linear process.

    Parameters
    ----------
    N: `int`
        Number of data points to simulate
    beta_0: `float`
        Intercept
    beta_1: `float`
        Slope of univariate predictor, X

    Returns
    -------
    df: `pd.DataFrame`
        A DataFrame containing the x and y values.
    """
    # Create a pandas DataFrame with column 'x' containing
    # N uniformly sampled values between 0.0 and 1.0
    df = pd.DataFrame(
        {"x": 
            np.linspace(start, stop, num=N)
        }
    )
    # Use a linear model (y ~ beta_0 + beta_1*x + epsilon) to 
    # generate a column 'y' of responses based on '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):
    """
    Plot the simulated data with sns.lmplot()

    Parameters
    ----------
    df: `pd.DataFrame`
        A DataFrame containing the x and y values.
    """
    # Plot the data, and a frequentist linear regression fit
    # using the seaborn package
    sns.lmplot(x="x", y="y", data=df, height=10)
    plt.xlim(0.0, 1.0)
    plt.show()


if __name__ == "__main__":
    # These are our "true" parameters
    beta_0 = 1.0  # Intercept
    beta_1 = 2.0  # Slope

    # Simulate 100 data points between 0 and 1, with a variance of 0.5
    start = 0
    stop = 1
    N = 100
    eps_mean = 0.0
    eps_sigma_sq = 0.5

    # Fix Random Seed
    s = 42

    # Simulate the "linear" data using the above parameters
    df = simulate_linear_data(
        start, stop, N, beta_0, beta_1, eps_mean, eps_sigma_sq
    )
    # Plot the simulated data
    plot_simulated_data(df)
				
			
L’output è dato nella figura seguente:
trading-quantitativo-bayes-linear-model-sim

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.

				
					# NEW
import bambi as bmb
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
# NEW
import pymc as pm
import seaborn as sns

...


def glm_mcmc_inference(df, iterations=5000):
    """
    Calculates the Markov Chain Monte Carlo trace of
    a Generalised Linear Model Bayesian linear regression 
    model on supplied data.

    Parameters
    ----------
    df: `pd.DataFrame`
        DataFrame containing the data
    iterations: `int`
        Number of iterations to carry out MCMC for
    """
    # Create the glm using the Bambi model syntax
    model = bmb.Model("y ~ x", df)

    # Fit the model using a 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):
    """
    Plot the trace generated from fitting the model. 

    Parameters
    ----------
    trace: `tracepymc.backends.base.MultiTrace`
        A MultiTrace or ArviZ InferenceData object that contains the samples.
    """
    pm.plot_trace(trace)
    plt.tight_layout()
    plt.show()


if __name__ == "__main__":
    # These are our "true" parameters
    beta_0 = 1.0  # Intercept
    beta_1 = 2.0  # Slope

    # Simulate 100 data points between 0 and 1, with a variance of 0.5
    start = 0
    stop = 1
    N = 100
    eps_mean = 0.0
    eps_sigma_sq = 0.5

    # Fix Random Seed
    s = 42

    # Simulate the "linear" data using the above parameters
    df = simulate_linear_data(
        start, stop, N, beta_0, beta_1, eps_mean, eps_sigma_sq
    )
    # Plot the simulated data
    plot_simulated_data(df)

    # NEW
    # Fit the GLM 
    trace = glm_mcmc_inference(df, iterations=5000)
    # Plot the GLM
    plot_glm_model(trace))
				
			
L’output dello script è il seguente:
				
					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.
				
			
Il traceplot è riportato nella figura seguente:
trading-quantitativo-bayes-linear-model-traceplot

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):
    """
    Plot the simulated data with true and estimated regression lines.

    Parameters
    ----------
    trace: `tracepymc.backends.base.MultiTrace`
        A MultiTrace or ArviZ InferenceData object that contains the samples.
    df: `pd.DataFrame`
        DataFrame containing the data
    N: `int`
        Number of data points to simulate
    """
    fig, ax = plt.subplots(figsize=(7, 7))
    # define x axis ticks
    x = np.linspace(0, 1, N)
    # plot simulated data observations
    ax.scatter(df['x'], df['y'])
    # extract slope and intercept draws from PyMC trace
    intercepts = trace.posterior.Intercept.to_numpy()[0]
    slopes = trace.posterior.x.to_numpy()[0]
    # plot 100 random samples from the slope and intercept draws
    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)
    # plot true regression line
    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__":
    # These are our "true" parameters
    beta_0 = 1.0  # Intercept
    beta_1 = 2.0  # Slope

    # Simulate 100 data points between 0 and 1, with a variance of 0.5
    start = 0
    stop = 1
    N = 100
    eps_mean = 0.0
    eps_sigma_sq = 0.5

    # Fix Random Seed
    s = 42

    # Simulate the "linear" data using the above parameters
    df = simulate_linear_data(
        start, stop, N, beta_0, beta_1, eps_mean, eps_sigma_sq
    )
    # Plot the simulated data
    plot_simulated_data(df)

    # Fit the GLM 
    trace = glm_mcmc_inference(df, iterations=5000)
    # Plot the GLM
    plot_glm_model(trace)
    # NEW
    plot_regression_lines(trace, df, N)
				
			
Nella seguente figura possiamo vedere l’intervallo campionato delle linee di regressione a posteriori:
trading-quantitativo-bayes-linear-model-post-samples
Da quanto sopra possiamo notare la presenza di incertezza nella posizione della retta di regressione campionata dal modello bayesiano. Tuttavia, si può vedere che l’intervallo è relativamente ristretto e che l’insieme dei campioni non è troppo dissimile dalla “vera” linea di regressione stessa.

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

SCRIVIMI SU TELEGRAM

Per informazioni, suggerimenti, collaborazioni...

Se è la prima volta che atterri su DataTrading, BENVENUTO!

Lascia che mi presenti. Sono Gianluca, ingegnere, trekker e appassionato di coding, finanza e trading. Leggi la mia storia.

Ho creato DataTrading per aiutare le altre persone ad utilizzare nuovi approcci e nuovi strumenti, ed applicarli correttamente al mondo del trading.

DataTrading vuole essere un punto di ritrovo per scambiare esperienze, opinioni ed idee.

TUTORIAL

TOOLBOX

Gli altri articoli di questa serie

Torna su