TUTORIAL DATA SCIENCE

Metodo della Massima Verosimiglianza per la Regressione Lineare

Scritto da
Gianluca

TUTORIAL

Lo scopo di questa serie di articoli è quello di descrivere una tecnica molto familiare, la regressione lineare, in un contesto matematico più rigoroso sotto un’interpretazione probabilistica del machine learning. Questo permette di comprendere il quadro delle probabilità che verrà successivamente utilizzato per modelli di .apprendimento supervisionato più complessi, in un contesto più semplice.

Iniziamo definendo la regressione lineare multipla, inserendola in un framework di apprendimento supervisionato probabilistico e ricavando una stima ottimale per i suoi parametri attraverso una tecnica nota come metodo di massima verosimiglianza.

Negli articoli successivi descriviamo i metodi per ridurre o mitigare la dimensione di determinati set di dati attraverso i concetti di selezione e restringimento di sottoinsiemi. Inoltre usiamo la libreria Scitkit-Learn di Python per dimostrare la regressione lineare, la selezione di sottoinsiemi e il restringimento.

Molte di queste tecniche si possono facilmente applicare a modelli più sofisticati e ci aiutano notevolmente nella creazione di metodi statistici efficaci e robusti per lo sviluppo di strategie di trading.

Questo articolo è molto più rigoroso dal punto di vista matematico rispetto ad altri articoli. E’ necessario per introdurre un più avanzato meccanismo probabilistico che pervade la ricerca sull’apprendimento automatico. Dopo aver visto alcuni esempi di modelli più semplici in questo framework, è più facile iniziare a consultare i documenti accademici di Machine Learning più avanzati, utili per nuove idee di trading.

Regressione lineare

La regressione lineare è una delle tecniche statistiche più familiari e dirette. Viene spesso insegnato al liceo, anche se in modo semplificato. Di solito è anche la prima tecnica presa in considerazione quando si studia l’apprendimento supervisionato poiché solleva questioni importanti che riguardano molti altri modelli supervisionati.

La regressione lineare afferma che il valore della risposta \(y\) è una funzione lineare delle feature d’ingresso \({\bf x}\), cioè:

\(\begin{eqnarray}y ({\bf x}) = \beta^T {\bf x} + \epsilon = \sum_{j=0}^p \beta_j x_j + \epsilon\end{eqnarray}\)

Dove \(\beta^T, {\bf x} \in \mathbb{R}^{p+1}\) e \(\epsilon \sim \mathcal{N}(\mu, \sigma^2)\), cioè \(\beta^T\) e \(x\) sono entrambi vettori di dimensione p+1 e \(\epsilon\), l’errore o il residuo, ha distribuzione normale con media \(\mu\) e varianza \(\sigma^2\). \(\epsilon\) rappresenta la differenza tra le previsioni fatte dalla regressione lineare e il vero valore della variabile di risposta.

Da notare che \(\beta^T\), cioè la trasposizione del vettore \(\beta\), e \(x\) hanno entrambi dimensione p+1, piuttosto che dimensione p, perché dobbiamo includere un termine di intercetta. \(\beta^T = (\beta_0, \beta_1, \ldots, \beta_p)\), mentre \({\bf x} = (1, x_1, \ldots, x_p)\). Dobbiamo includere ‘1’ in \(x\) come “trucco” notazionale.

Interpretazione probabilistica

Un modo alternativo per descrivere la regressione lineare è considerarla come un modello di probabilità congiunta [2] , [3]. In altre parole siamo interessati alla probabilità congiunta di avere il comportamento della risposta y condizionato ai valori del vettore delle feature x, nonché a altri parametri del modello (descritti dal vettore \({\bf \theta}\)). Quindi siamo interessati a un modello nella forma \(p(y \mid {\bf x}, {\bf \theta})\), cioè un modello di densità di probabilità condizionata (CPD).

La regressione lineare può essere scritta come CPD nel modo seguente:

\(\begin{eqnarray}p(y \mid {\bf x}, {\bf \theta}) = \mathcal (y \mid \mu({\bf x}), \sigma^2 ({\bf x}))\end{eqnarray}\)

Per la regressione lineare assumiamo che \(\mu({\bf x})\) sia lineare e quindi \(\mu ({\bf x}) = \beta^T {\bf x}\). Dobbiamo anche assumere che la varianza nel modello sia fissa (cioè che non dipenda da x) e quindi \(\sigma^2 ({\bf x}) = \sigma^2\) è una costante. Ciò implica che il nostro vettore di parametri corrisponda a \(\theta = (\beta, \sigma^2)\).

Ricordiamo che abbiamo usato una tale interpretazione probabilistica quando abbiamo descritto la regressione lineare bayesiana in un precdente articolo.

Generalizzando l’interpretazione del modello in questo modo abbiamo il vantaggio di poter facilmente vedere come altri modelli, specialmente quelli che gestiscono le non linearità, si inseriscono nello stesso quadro probabilistico. Questo permette di ottenere risultati attraverso i modelli utilizzando tecniche simili.

Se consideriamo \({\bf x} = (1, x)\), possiamo creare un grafico bidimensionale \({\bf x} = (1, x)\) per x e y in modo da rappresentare graficamente questa distribuzione congiunta. A tale scopo dobbiamo correggere i parametri \(\beta = (\beta_0, \beta_1)\) e \(\sigma^2\) (che costituiscono i parametri di \(\theta\)). Di seguito uno script Python che usa matplotlib per visualizzare la distribuzione:

				
					
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import norm


if __name__ == "__main__":
    # Set up the X and Y dimensions
    fig = plt.figure()
    ax = fig.gca(projection='3d')
    X = np.arange(0, 20, 0.25)
    Y = np.arange(-10, 10, 0.25)
    X, Y = np.meshgrid(X, Y)

    # Create the univarate normal coefficients
    # of intercep and slope, as well as the
    # conditional probability density
    beta0 = -5.0
    beta1 = 0.5
    Z = norm.pdf(Y, beta0 + beta1*X, 1.0)

    # Plot the surface with the "coolwarm" colormap
    surf = ax.plot_surface(
        X, Y, Z, rstride=1, cstride=1, cmap=cm.coolwarm,
        linewidth=0, antialiased=False
    )

    # Set the limits of the z axis and major line locators
    ax.set_zlim(0, 0.4)
    ax.zaxis.set_major_locator(LinearLocator(5))
    ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))

    # Label all of the axes
    ax.set_xlabel('X')
    ax.set_ylabel('Y')
    ax.set_zlabel('P(Y|X)')

    # Adjust the viewing angle and axes direction
    ax.view_init(elev=30., azim=50.0)

    # Plot the probability density
    plt.show()

È chiaro che la risposta \(y\) è linearmente dipendente da \(x\). Tuttavia siamo anche in grado di verificare l’elemento probabilistico del modello dato che la probabilità si distribuisce normalmente attorno alla risposta lineare.

Estensione della funzione di base

L’uso dell’interpretazione probabilistica ha il vantaggio di verificare facilmente come modellare relazioni non lineari, semplicemente sostituendo il vettore delle feature \(x\) con una specifica funzione di trasformazione \(\phi({\bf x})\):

\(\begin{eqnarray}p(y \mid {\bf x}, {\bf \theta}) = \mathcal(y \mid \beta^T \phi({\bf x}), \sigma^2)\end{eqnarray}\)

Per \({\bf x} = (1, x_1, x_2, x_3)\) possiamo creare una \(\phi\) che include termini di ordine superiore, inclusi termini incrociati, ad es:

\(\begin{eqnarray}\phi({\bf x}) = (1, x_1, x_1^2, x_2, x^2_2, x_1 x_2, x_3, x_3^2, x_1 x_3, \ldots)\end{eqnarray}\)

Da sottolineare che questa funzione non è lineare nelle feature, \(x\), ma è ancora lineare nei parametri, \({\bf \beta}\), quindi può essere considerata una regressione lineare.

Tale modifica, utilizzando una funzione di trasformazione \(\phi\), è nota come estensione della funzione di base e può essere utilizzata per generalizzare la regressione lineare a molte configurazioni di dati non lineari.

Esistono molte altre tecniche di apprendimento supervisionato più sofisticate per catturare le non linearità. Abbiamo descritto una di queste tecniche, Support Vector Machines con il “trucco del kernel”, in questo articolo.

Metodo della massima verosimiglianza

In questa paragrafo descriviamo come i coefficienti di regressione lineare ottimali, ovvero il parametro delle componenti \(\beta\), sono scelti per adattarsi al meglio ai dati. Nel caso univariato questo è spesso noto come “trovare la retta di migliore adattamento”. Tuttavia, siamo in un caso multivariato, come il nostro vettore di feature \({\bf x} \in \mathbb{R}^{p+1}\), quindi si tratta di “trovare l’iperpiano p-dimensionale di migliore adattamento”!

Lo strumento principale per trovare i parametri dei modelli statistici è noto come metodo della massima verosimiglianza (MLE). E’ stato introdotto brevemente nell’articolo sul Deep Learning e la regressione logistica.

Verosimiglianza e verosimiglianza logaritmica negativa

L’idea fondamentale è, se i dati sono stati generati dal modello, quali parametri sono stati usati con maggiore probabilità? Cioè, qual è la probabilità di avere i dati \(\mathcal{D}\), dato uno specifico insieme di parametri \({\bf \theta}\)? Anche questo caso è un problema di densità di probabilità condizionata. Vogliamo trovare i valori di \({\bf \theta}\) che massimizzano \(p(\mathcal{D} \mid {\bf \theta})\). Questo CPD è noto come verosimiglianza e abbiamo descritto alcuni esempi nell’articolo introduttivo sulle statistiche bayesiane.

Questo problema può essere formulato come la ricerca della moda di \(p(\mathcal{D} \mid {\bf \theta})\), che corrisponde a \(\hat{{\bf \theta}}\). Per facilitare il calcolo, cerchiamo di massimizzare il logaritmo naturale del CPD, invece dello stesso CPD:

\(\begin{eqnarray}\hat{{\bf \theta}} = \text{argmax}_{\theta} \log p(\mathcal{D} \mid {\bf \theta})\end{eqnarray}\)

Nei problemi di regressione lineare dobbiamo assumere che i vettori delle feature siano tutti indipendenti e identicamente distribuiti (iid). Questo rende molto più semplice risolvere il problema della log-verosimiglianza, usando le proprietà dei logaritmi naturali. Poiché dobbiamo differenziare questi valori, è molto più facile differenziare una somma rispetto a un prodotto, da cui il logaritmo:

\(\begin{eqnarray}\mathcal{l}({\bf \theta}) &:=& \log p(\mathcal{D} \mid {\bf \theta}) \\
&=& \log \left( \prod_{i=1}^{N} p(y_i \mid {\bf x}_i, {\bf \theta}) \right) \\
&=& \sum_{i=1}^{N} \log p(y_i \mid {\bf x}_i, {\bf \theta})\end{eqnarray}\)

Come descritto nell’articolo sul Deep Learning e la Regressione Logistica, per motivi di maggiore facilità di calcolo, è spesso consigliabile ridurre al minimo il negativo della log-verosimiglianza piuttosto che massimizzare il log-verosimiglianza. Quindi, possiamo “mettere un segno meno davanti alla log-verosimiglianza” per ottenere la log-verosimiglianza negativa (NLL):

\(\begin{eqnarray}\text{NLL} ({\bf \theta}) = – \sum_{i=1}^{N} \log p(y_i \mid {\bf x}_i, {\bf \theta})\end{eqnarray}\)

Questa è la funzione che dobbiamo minimizzare. In questo modo deriviamo la stima dei minimi quadrati ordinari per i coefficienti \(\beta\).

Minimi quadrati ordinari

Il nostro obiettivo è ricavare l’insieme ottimale di coefficienti \(\beta\) che “molto probabilmente” hanno generato i dati per il nostro problema di addestramento. Questi coefficienti devono formare un iperpiano di “best fit” attraverso i dati di allenamento. Il processo da seguire è dato da:

Utilizzare la definizione della distribuzione normale per estendere la funzione di verosimiglianza logaritmica negativa.
Utilizzare le proprietà dei logaritmi per riformularlo in termini di somma residua dei quadrati (RSS), che è equivalente alla somma di ogni residuo attraverso tutte le osservazioni.
Riscrivere i residui in forma matriciale, creando la matrice X dei dati, che ha \(N \times (p+1)\) dimensioni, e formulare l’RSS come un’equazione di matrici.
Differenziare questa equazione matriciale rispetto a al vettore \(\beta\) dei parametri e impostare l’equazione a zero (con alcune ipotesi su X)
Risolvere l’equazione per \(\beta\) per ricavare \(\hat{\beta}_\text{OLS}\), la stima dei minimi quadrati ordinari (OLS).

Di seguito seguiamo da vicino gli approcci descritti in [2] e [3] . Il primo passo è espandere l’NLL utilizzando la formula per una distribuzione normale:

\(\begin{eqnarray}\text{NLL} ({\bf \theta}) &=& – \sum_{i=1}^{N} \log p(y_i \mid {\bf x}_i, {\bf \theta}) \\
&=& – \sum_{i=1}^{N} \log \left[ \left(\frac{1}{2 \pi \sigma^2}\right)^{\frac{1}{2}} \exp \left( – \frac{1}{2 \sigma^2} (y_i – {\bf \beta}^{T} {\bf x}_i)^2 \right)\right] \\
&=& – \sum_{i=1}^{N} \frac{1}{2} \log \left( \frac{1}{2 \pi \sigma^2} \right) – \frac{1}{2 \sigma^2} (y_i – {\bf \beta}^T {\bf x}_i)^2 \\
&=& – \frac{N}{2} \log \left( \frac{1}{2 \pi \sigma^2} \right) – \frac{1}{2 \sigma^2} \sum_{i=1}^N (y_i – {\bf \beta}^T {\bf x}_i)^2 \\
&=& – \frac{N}{2} \log \left( \frac{1}{2 \pi \sigma^2} \right) – \frac{1}{2 \sigma^2} \text{RSS}({\bf \beta})\end{eqnarray}\)

Dove \(\text{RSS}({\bf \beta}) := \sum_{i=1}^N (y_i – {\bf \beta}^T {\bf x}_i)^2\)
è la somma residua dei quadrati, nota anche come somma degli errori quadrati (SSE).

Poiché il primo termine nell’equazione è una costante, dobbiamo semplicemente preoccuparci di minimizzare l’RSS, che sarà sufficiente per produrre la stima del parametro ottimale.

Per semplificare la formula, possiamo scrivere quest’ultimo termine in forma matriciale. Considerando \(N \times (p+1)\) della matrice X possiamo scrivere il termine RSS come:

\(\begin{eqnarray}\text{RSS}({\bf \beta}) = ({\bf y} – {\bf X}{\bf \beta})^T ({\bf y} – {\bf X}{\bf \beta})\end{eqnarray}\)

A questo punto vogliamo ora differenziare questo termine rispetto alla variabile dei parametr \({\bf \beta}\):

\(\begin{eqnarray}\frac{\partial RSS}{\partial \beta} = -2 {\bf X}^T ({\bf y} – {\bf X} \beta)\end{eqnarray}\)

C’è un presupposto fondamentale da fare. Abbiamo bisogno che \({\bf X}^T {\bf X}\) sia definito positivamente, che si verifica solo nel caso in cui ci sono più osservazioni che dimensioni. Se questo non è il caso (cosa estremamente comune nelle configurazioni ad alta dimensione), non è possibile trovare un insieme univoco di coefficienti \(\beta\) e quindi la conseguente equazione matriciale non sarà valida.

Nell’ipotesi di una \({\bf X}^T {\bf X}\) definita positivamente possiamo impostare l’equazione differenziata a zero e ricavare \(\beta\):

\(\begin{eqnarray}{\bf X}^T ({\bf y} – {\bf X} \beta) = 0\end{eqnarray}\)

La soluzione a questa equazione matriciale fornisce \(\hat{\beta}_\text{OLS}\):

\(\begin{eqnarray}\hat{\beta}_\text{OLS} = ({\bf X}^{T} {\bf X})^{-1} {\bf X}^{T} {\bf y}\end{eqnarray}\)

Prossimi passi

Ora che abbiamo descritto la procedura MLE per produrre le stime OLS, possiamo discutere cosa succede quando ci troviamo in un ambiente ad alta dimensione (come spesso accade con i dati del mondo reale) e quindi la nostra matrice \({\bf X}^T {\bf X}\) non ha l’inverso. In questo caso è necessario utilizzare tecniche di selezione di sottoinsiemi e restringimento per ridurre la dimensionalità del problema. Questo è l’argomento del prossimo articolo.

Nota bibliografica

Un’introduzione elementare alla regressione lineare, così come al restringimento, alla regolarizzazione e alla riduzione della dimensione, nel quadro dell’apprendimento supervisionato, può essere trovata in [1] . Una spiegazione molto più rigorosa delle tecniche, compresi i recenti sviluppi, si può trovare in [2] . Un approccio probabilistico (principalmente bayesiano) alla regressione lineare, insieme a una derivazione completa della stima di massima verosimiglianza tramite i minimi quadrati ordinari e un’ampia discussione sulla contrazione e sulla regolarizzazione, può essere trovata in [3] . Una panoramica basata su esempi del “mondo reale” della regressione lineare in un regime di alta linearità, con un’ampia discussione sulla riduzione della dimensione e sui minimi quadrati parziali può essere trovata in [4].