Partiamo da quanto visto nell’ultimo articolo: la regressione lineare multipla. Con questa vengono utilizzate diverse caratteristiche di input, ognuna con il proprio peso, per predire un valore di output.

Per semplificarci la vita, consideriamo per ora di avere a disposizione un’unica caratteristica in input, ma che la sua nuvola di valori in ingresso non abbia una tendenza lineare ma segua piuttosto quella di una curva, più o meno complessa.

Prendiamo come esempio il valore LSTAT del nostro solito dataset, che come avevamo visto la volta scorsa, aveva un andamento leggermente curvo. Nonostante questo l’avevamo comunque inclusa nell’esempio di regressione lineare. Vediamo se possiamo ottenere risultati migliori.

Realisticamente se dovessimo usare questa caratteristica in una regressione lineare semplice, avremmo mediamente un errore piuttosto alto, proprio perchè una retta non può mediare bene dei punti distribuiti in quella maniera. Dovremmo fare in modo che la retta di regressione curvi (e quindi non sia più una retta…) per seguire l’andamento dei punti in modo simile a quanto riportato in figura con la linea rossa. Ma come possiamo fare?

Per fare una cosa del genere possiamo usare un polinomio della nostra caratteristica, cioè un’equazione simile alla seguente.

Come potete vedere, il polinomio non è poi molto diverso dall’equazione della regressione lineare multipla: in sostanza posso prendere proprio la regressione lineare aggiungendo nuovi termini come fossero nuove caratteristiche. In verità la caratteristica è sempre la stessa (in questo esempio sarà LSTAT) che semplicemente viene elevata a potenze mano a mano crescenti, fino ad n, detto grado del polinomio.

NOTA: il grado di una regressione polinomiale è un iperparametro del modello di regressione.

Non vi spaventate: in verità non ci stiamo complicando la vita più di tanto, perchè se vi ricordate dovremo risolvere l’equazione secondo i pesi w e non secondo i valori x. Quindi, dal punto di vista dei pesi w, questa equazione rimane lineare! In definitiva, quindi, una regressione polinomiale si risolve in modo del tutto analogo ad una regressione lineare multipla.

Purtroppo, come spesso accade, non è tutt’oro quello che luccica! Se aumentiamo il grado del polinomio, in generale l’approssimazione della curva di regressione migliorerà, ma si può incorrere facilmente nel problema dell’overfitting, dove la curva si adatta fin troppo bene al dataset di training, al punto che perderà di generalizzazione e quindi darà buoni, anzi ottimi risultati solo sui dati di allenamento e risultati pessimi su altri.

Il comportamento è riportato nella figura a lato, dove il polinomio ad alto grado in rosso, passa da tutti i punti del training set, quindi ottenendo un errore piccolissimo su questi esempi, ma si allontana dall’andamento ideale indicato dalla curva in verde, ottenendo quindi risultati peggiori su esempi che non rientrano nel training set.

Codice di esempio

NOTA: come al solito il codice di esempio lo trovate qui

Il codice di esempio non cambia molto nella sostanza, semplicemente andremo a verificare quale grado di polinomio consente di avere un risultato migliore in questo caso.

Cominciamo quindi con il prendere la sola colonna LSTAT del solito dataset come input del modello. Poi andiamo a creare i due dataset di train e di test, ma questa volta, per avere un confronto diretto, andiamo ad impostare il parametro random_state della funzione di splitting. Questo consente di mischiare i dati, ma sempre allo stesso modo, fornendo un seed fisso al generatore casuale (può essere un numero qualsiasi, io ho scelto 1234). Così facendo gli esperimenti saranno ripetibili e direttamente confrontabili tra loro.

# importiamo il solito dataset al completo
dataset = pd.read_csv("https://archive.ics.uci.edu/ml/machine-learning-databases/housing/housing.data",
                      sep="\s+",
                      names=["CRIM", "ZN", "INDUS", "CHAS", "NOX", "RM", "AGE", "DIS", "RAD", "TAX", 
                             "PRATIO", "B", "LSTAT", "MEDV"])
 
print("Utilizzo solo LSTAT")
print()
 
# associamo ad X i valori di input di LSTAT
# associamo ad Y i valori di output
X = dataset[["LSTAT"]].values
Y = dataset["MEDV"].values
 
# suddividiamo il dataset in due dataset, uno di training ed uno di test
# questa volta uso un random state fisso, in modo che i dati vengano mischiati, 
# ma sempre allo stesso modo
# così da poter confrontare direttamente i vari metodi di regressione
# in pratica random_state è il seed del generatore random
X_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size=0.3, random_state=1234)

Per vedere come si comporta la regressione al variare del grado del polinomio, eseguiamo un ciclo prendendo il range da 1 a 10:

# ciclo il procedimento per vedere i risultati dei vari gradi del polinomio
# in questo caso vado da grado 1 (lineare) a grado 10
for grado in range(1, 11):

Creiamo artificialmente delle ulteriori feature di grado superiore: in SciKitLearn sfruttiamo la classe PolynomialFeatures che si occupa appunto di creare queste feature per noi.

# istanzo il generatore di feature polinomiali per il grado corrente
# successivamente lo utilizzo per generare le feature per i set di training e di test
polyfeatures = PolynomialFeatures(degree=grado)
X_train_poly = polyfeatures.fit_transform(X_train)
X_test_poly = polyfeatures.transform(X_test)

Se andassimo a guardare il contenuto di x_train_poly, vedremmo i valori delle x replicati ed elevati a potenza per ogni grado, dal grado 0 (che sarà sempre uguale ad 1) fino al grado del ciclo corrente. Ad esempio, se fossimo al ciclo con grado 3, andando a guardare il primo esempio di train vedremmo:

> X_train_poly[0]
Out[2]: array([ 1. , 6.19 , 38.3161 , 237.176659])

dove per la potenza 0 abbiamo 1., per la potenza 1 abbiamo 6.19 (che è il normale valore di x), per la potenza 2 abbiamo 38.3161 (che è x al quadrato) e per la potenza 3 abbiamo 237.176659 (che è appunto x elevato alla terza).

Visto che elevare a potenza porta solitamente a numeri molto grandi o molto piccoli, conviene sempre eseguire una standardizzazione dei set così ottenuti, in modo da evitare che il modello consideri più del dovuto le caratteristiche con valori alti ed ignori quelle con valori bassi. In questo caso, avendo un’unica caratteristica, la standardizzazione non modifica sostanzialmente i risultati (provare per credere), ma è bene fare l’abitudine a eseguire queste trasformazioni dei set di dati. Una volta create e standardizzate le feature polinomiali, si esegue l’algoritmo di risoluzione già visto per la regressione lineare multipla, con il calcolo di MSE e R2. Eseguendo il ciclo, possiamo vedere i risultati per un confronto ed un’analisi abbastanza immediata.

Grado: 1 - MSE: 0.561293110009808 - R2: 0.561293110009808
Grado: 2 - MSE: 0.6561389171118859 - R2: 0.6561389171118859
Grado: 3 - MSE: 0.6770922322673485 - R2: 0.6770922322673485
Grado: 4 - MSE: 0.6958607823010035 - R2: 0.6958607823010035
Grado: 5 - MSE: 0.6826978684220166 - R2: 0.6826978684220166
Grado: 6 - MSE: 0.6749390736563432 - R2: 0.6749390736563432
Grado: 7 - MSE: 0.6747390224812395 - R2: 0.6747390224812395
Grado: 8 - MSE: 0.6754703185274282 - R2: 0.6754703185274282
Grado: 9 - MSE: 0.6773603201353989 - R2: 0.6773603201353989
Grado: 10 - MSE: 0.6792253949049618 - R2: 0.6792253949049618

Come si può notare, il modello migliora con l’aumentare del polinomio, fino ad grado 4. Dal grado 5 in poi, i risultati hanno un peggioramento, dovuto appunto all’overfitting di cui parlavamo in precedenza. In pratica dal grado 5 in poi, il modello di regressione è più complesso del problema e quindi predilige adattarsi ad ogni singolo esempio del training set piuttosto che ricercare una regola generalizzata. E’ un po’ come se uno studente preferisse imparare a memoria i compiti piuttosto che imparare le regole per risolverli.

Se ora replichiamo il codice per elaborare tutte le colonne del dataset (fermandoci al grado 4 che abbiamo visto essere un limite per questo caso), il risultato sarà il seguente:

Grado: 1 - MSE: 0.738410984790933 - R2: 0.738410984790933
Grado: 2 - MSE: 0.8492301327888548 - R2: 0.8492301327888548
Grado: 3 - MSE: -125.96148346549076 - R2: -125.96148346549076
Grado: 4 - MSE: -10.802897842515257 - R2: -10.802897842515257

Come si vede la regressione polinomiale questa volta regge fino al grado 2, per poi subire un netto calo. Questo è dovuto alla grande complessità di una regressione di questo tipo: per ogni caratteristica vengono create le feature polinomiali, ottenendo quindi un polinomio veramente complesso.

Si noti che però già fin dal grado 1 la regressione dà risultati migliori della regressione polinomiale di grado 4 con la sola caratteristica LSTAT. Questo perchè, in generale, più caratteristiche consentono alla regressione di comportarsi meglio.

Soluzione all’overfitting: la regolarizzazione

Visti gli esempi di poco fa, se ne può dedurre che per contrastare l’overfitting una soluzione sia quella di diminuire la complessità del modello, ad esempio rimuovendo alcune caratteristiche da elaborare o riducendo il grado del polinomio. In verità questo non è sempre possibile, perchè magari siamo nel caso iniziale dove avevamo solo LSTAT o perchè magari rimuovere alcune caratteristiche dal set in input può incidere troppo sui risultati del modello o perchè riducendo anche di solo un’unità il grado del polinomio, i risultati comunque peggiorano in modo non accettabile.

In questo caso si può agire in modo diverso mantenendo tutte le caratteristiche dei dati ed il grado del polinomio: andremo a penalizzare i pesi con valori più elevati e quindi più destabilizzanti per il modello (perchè piccole variazioni in ingresso si ripercuotono in modo più pesante sull’uscita). Questa soluzione è chiamata regolarizzazione e si attua aggiungendo alla funzione costo un nuovo termine, nel modo seguente:

dove λ è un nuovo iperparametro del modello ed indica quanto influisce la funzione di regolarizzazione R(W) sulla nuova funzione costo. Esistono divese funzioni di regolarizzazione, noi ne vedremo due in particolare, chiamate L1 e L2.

La regolarizzazione L1 è la somma dei valori assoluti dei pesi. Questo obbliga l’algoritmo di ottimizzazione della funzione costo a trovare il minimo riducendo a 0 i pesi delle caratteristiche meno influenti, quindi in sostanza a selezionare solo le proprietà del dataset che risultano essere più importanti nel raggiungimento dello scopo dell’elaborazione.

La regolarizzazione L2 è la somma dei pesi al quadrato. Viene chiamata anche weight decay in quanto costringe l’algoritmo di ottimizzazione della funzione costo a trovarne il minimo penalizzando i pesi più grandi, stabilizzando e uniformando maggiormente il comportamento della funzione stessa.

I valori dell’iperparametro λ influiscono sul risultato finale:

• con un valore 0, annulleremo l’efficacia della regolarizzazione e quindi non ci sarà alcun effetto (se il modello è in overfitting, ci resterà)
• con un valore troppo alto, l’influenza delle regolarizzazione sarà elevata e la maggior parte dei pesi tenderà a 0, quindi il modello andrà in underfitting, perchè non riuscirà ad adattarsi al set di dati.
• un metodo utile a trovare il valore ottimale di λ è quello di ricercarlo in uno spazio di potenze di 10 che va da 0.0001 a 10 (valori tipici). Una volta trovato un valore buono, si può eventualmente procedere ad affinarlo ulteriormente.

In generale la regolarizzazione L2 consente di ottenere risultati migliori della L1, ma come sempre dipenderà dalla natura del problema e dei dati. Esistono diversi algoritmi che utilizzano queste regolarizzazioni. Noi vedremo il Lasso che utilizza L1, il Bridgeche utilizza L2 ed ElasticNet che li utilizza entrambi.

Codice di esempio

NOTA: il codice di questo esempio lo trovi qui

Riprendiamo il solito esempio con il dataset del valore delle abitazioni. Come abbiamo visto nel codice precedente, utilizzando tutte le caratteristiche il modello andava decisamente in overfitting già con il grado 3. Partiamo quindi da questa situazione, evidenziando i valori di MSE ed R2 sia per il training set che per il test set. Si vede che mentre il training è perfetto (R2 = 1.0), il test è decisamente pessimo (R2 circa -126). Ciò conferma con forza l’overfitting: il modello si comporta perfettamente solo con i dati di allenamento.

Modello in overfitting già al grado 3
TRAIN --- MSE: 5.51557399807701e-23 - R2: 1.0
TEST ---- MSE: 11703.304270441036 - R2: -125.96148346549076

Ora vogliamo analizzare come si evolve il modello utilizzando diversi valori di λ. Da notare che in SciKitLearn l’iperparametro è chiamato alpha.

# preparo un array con i vari valori di lambda che vogliamo provare
# NOTA: "alpha is the new lambda"... almeno in SciKitLearn...
alphas = [0.0001, 0.001, 0.01, 0.1, 1., 10.]

Partiamo con analizzare la regolarizzazione L1, utilizzando un modello di tipo Lasso che la utilizza. Cicliamo quindi su questi valori, costruendo il modello con il valore alpha corrente e verificando l’andamento di MSE ed R2 sia per training che per test.

for alpha in alphas:
    model = Lasso(alpha=alpha)
    model.fit(X_train_poly_std, Y_train)
  
    Y_pred_train = model.predict(X_train_poly_std)
    Y_pred_test = model.predict(X_test_poly_std)
 
    # visualizzo i risultati
    mse_train = mean_squared_error(Y_train, Y_pred_train)
    r2_train = r2_score(Y_train, Y_pred_train)
    mse_test = mean_squared_error(Y_test, Y_pred_test)
    r2_test = r2_score(Y_test, Y_pred_test)
    print("Alpha:", alpha)
    print("TRAIN --- MSE:", mse_train, " - R2:", r2_train)
    print("TEST ---- MSE:", mse_test, " - R2:", r2_test)
    print()

I risultati appaiono molto interessanti, già con un alpha molto piccolo.

Alpha: 0.0001
TRAIN --- MSE: 2.4562994734900747 - R2: 0.9697055210529493
TEST ---- MSE: 18.257750879280177 - R2: 0.8019336178218133

Alpha: 0.001
TRAIN --- MSE: 2.571157433193907 - R2: 0.9682889339959952
TEST ---- MSE: 15.07945015215896 - R2: 0.8364129209220437

Alpha: 0.01
TRAIN --- MSE: 4.259481193002142 - R2: 0.947466192653042
TEST ---- MSE: 10.233757276194645 - R2: 0.8889806694599022

Alpha: 0.1
TRAIN --- MSE: 10.943459778507348 - R2: 0.8650301335623247
TEST ---- MSE: 11.953011433548612 - R2: 0.8703296070566817

Alpha: 1.0
TRAIN --- MSE: 22.15011357124248 - R2: 0.7268141948891317
TEST ---- MSE: 19.100207088500945 - R2: 0.7927943621376276

Alpha: 10.0
TRAIN --- MSE: 81.08076319065403 - R2: 0.0
TEST ---- MSE: 92.2020720754774 - R2: -0.00023989625398068704

Il modello migliora tanto fino ad un alpha = 0.01, per poi degradare nuovamente, perchè mano a mano che alpha cresce, con L1diverse proprietà tenderanno ad essere ignorate dal modello.

Ora sperimentiamo la regolarizzazione L2. Il codice è esattamente lo stesso, ma invece della classe Lasso useremo la classe Ridgeche appunto usa L2.

Alpha: 0.0001
TRAIN --- MSE: 0.25039809765814675 - R2: 0.9969117446875857
TEST ---- MSE: 184.51331542152369 - R2: -1.0016641201254397

Alpha: 0.001
TRAIN --- MSE: 0.5515880625040435 - R2: 0.993197053890489
TEST ---- MSE: 49.99729451544472 - R2: 0.457612095331544

Alpha: 0.01
TRAIN --- MSE: 1.0761277840796133 - R2: 0.9867277052936811
TEST ---- MSE: 25.462605572016866 - R2: 0.7237728677630881

Alpha: 0.1
TRAIN --- MSE: 1.9022251024735395 - R2: 0.9765391317543395
TEST ---- MSE: 15.681778479011296 - R2: 0.8298786553724753

Alpha: 1.0
TRAIN --- MSE: 3.261859967499238 - R2: 0.9597702360074082
TEST ---- MSE: 11.206560698455538 - R2: 0.8784273622266199

Alpha: 10.0
TRAIN --- MSE: 6.069436692029711 - R2: 0.9251433206448985
TEST ---- MSE: 10.240505693312462 - R2: 0.88890746030217

In questo caso vediamo che il modello migliora anche con alpha = 10.0.

Ora proviamo ad utilizzare ElasticNet che è un modello che utilizza entrambe le regolarizzazioni. Attraverso il parametro l1_ratioè possibile indicare a quale regolarizzazione, tra L1 ed L2, dare più peso nel calcolo della regolarizzazione complessiva:

• l1_ratio = 0.5: le due regolarizzazioni sono utilizzate con lo stesso peso
• l1_ratio > 0.5: la regolarizzazione L1 ha peso man mano maggiore
• l1_ratio < 0.5: la regolarizzazione L1 ha peso man mano minore

Il codice di calcolo rimane lo stesso visto in precedenza, semplicemente si utilizzerà la classe ElasticNet per istanziare il modello.

Alpha: 0.0001
TRAIN --- MSE: 2.4594625050931245 - R2: 0.9696665101768971
TEST ---- MSE: 17.94089266646625 - R2: 0.8053710050602776

Alpha: 0.001
TRAIN --- MSE: 2.6175632431282003 - R2: 0.9677165934295261
TEST ---- MSE: 13.841321280174073 - R2: 0.849844570196145

Alpha: 0.01
TRAIN --- MSE: 4.762397978930566 - R2: 0.9412635279747895
TEST ---- MSE: 10.095068335777444 - R2: 0.8904852149462651

Alpha: 0.1
TRAIN --- MSE: 10.094327399066263 - R2: 0.8755028073017718
TEST ---- MSE: 11.939901080809046 - R2: 0.8704718326875043

Alpha: 1.0
TRAIN --- MSE: 20.722913241578453 - R2: 0.7444164013003872
TEST ---- MSE: 19.50835365378118 - R2: 0.7883666473485486

Alpha: 10.0
TRAIN --- MSE: 66.83301748568552 - R2: 0.1757228859756319
TEST ---- MSE: 75.31698600782207 - R2: 0.18293534434934888

Come si vede dai risultati, il modello migliora fino ad alpha=0.01, ma il risultato è decisamente migliore di entrambe le prove precedenti con le sole L1 ed L2, in quanto in questo caso entrambe concorrono a migliorare il modello, che sfiora un risultato ottimo.

See you soon!

Sperimentate ed esercitatevi con le varie regressioni. La prossima volta inizieremo a trattare le classificazioni, ma lo faremo cominciando con una… regressione…

…suspense….

Bye!