# Metodi di regressione con `scikit-learn`¶

In questo quaderno dimostriamo l'uso della libreria Python SciKit-Learn (`sklearn`) per l'applicazione di modelli di regressione di vario tipo.

In [1]:

    # Importiamo le varie librerie
    import numpy as np
    import pandas as pd
    import matplotlib.pyplot as plt
    %matplotlib inline
    import sklearn
    import sklearn.linear_model # per la regressione lineare
    import sklearn.neighbors    # per la regressione K-Nearest Neighbor
    from sklearn.metrics import mean_squared_error, r2_score  # misure di qualità (MSE, coefficiente R²)

## Regressione lineare per il dataset `Advertising`¶

In [2]:

    # Carico il dataset
    df = pd.read_csv('data/Advertising.csv')

In [3]:

    df

Out[3]:

| Unnamed: 0 | TV | radio | newspaper | sales
--- | --- | --- | --- | --- | ---
0 | 1 | 230.1 | 37.8 | 69.2 | 22.1
1 | 2 | 44.5 | 39.3 | 45.1 | 10.4
2 | 3 | 17.2 | 45.9 | 69.3 | 9.3
3 | 4 | 151.5 | 41.3 | 58.5 | 18.5
4 | 5 | 180.8 | 10.8 | 58.4 | 12.9
... | ... | ... | ... | ... | ...
195 | 196 | 38.2 | 3.7 | 13.8 | 7.6
196 | 197 | 94.2 | 4.9 | 8.1 | 9.7
197 | 198 | 177.0 | 9.3 | 6.4 | 12.8
198 | 199 | 283.6 | 42.0 | 66.2 | 25.5
199 | 200 | 232.1 | 8.6 | 8.7 | 13.4

200 rows × 5 columns

## Regressione lineare semplice¶

In [4]:

    # Estraiamo le colonne di interesse come dataframe individuali
    X = df[['TV']]
    y = df[['sales']]
    print(X.shape)
    print(y.shape)

    (200, 1)
    (200, 1)

In [5]:

    # Costruiamo e fittiamo un modello di regressione lineare
    regr = sklearn.linear_model.LinearRegression()
    regr.fit(X, y);

In [6]:

    # Calcoliamo le sue predizioni sui dati
    y_pred = regr.predict(X)

In [7]:

    # Grafichiamo i dati originali e l'ipotesi lineare
    plt.scatter(X, y, alpha=0.5)
    plt.plot(X, y_pred, color='red', linewidth=3)
    plt.show()

In [8]:

    # Stampiamo l'intercetta e il coefficiente della regressione
    print('Intercetta: \n', regr.intercept_, '\nCoefficienti: \n', regr.coef_)

    Intercetta:
     [7.03259355]
    Coefficienti:
     [[0.04753664]]

Ogni 1000$ in più di pubblicità TV sono associati mediamente con 47.5 unità di prodotto in più vendute.

In [9]:

    # Stampiamo l'errore quadratico medio
    print('Mean squared error (MSE):\n %.3f' % mean_squared_error(y, y_pred))
    # Stampiamo la statistica R^2
    print('Coefficiente R^2:\n %.3f' % r2_score(y, y_pred))

    Mean squared error (MSE):
     10.513
    Coefficiente R^2:
     0.612

In [10]:

    # Prediciamo le vendite con un budget per la pubblicità televisiva di 100 mila dollari
    regr.predict([[100]])

Out[10]:

    array([[11.78625759]])

Per un investimento TV di 100000$ si prevede una vendita media di 11786 unità di prodotto.

## Separazione di training set e test set¶

In [11]:

    from sklearn.model_selection import train_test_split
    train, test = train_test_split(df, test_size=0.2, random_state=0) # Split 80%-20%
    X_train = train[['TV']]
    y_train = train[['sales']]
    X_test = test[['TV']]
    y_test = test[['sales']]
    regr.fit(X_train, y_train)
    y_pred_train = regr.predict(X_train)
    print('Training R^2:\n %.3f' % r2_score(y_train, y_pred_train))
    print('Training MSE:\n %.3f' % mean_squared_error(y_train, y_pred_train))
    y_pred_test = regr.predict(X_test)
    print('Test MSE:\n %.3f' % mean_squared_error(y_test, y_pred_test))

    Training R^2:
     0.588
    Training MSE:
     10.617
    Test MSE:
     10.186

## Modello basato su `radio`¶

In [12]:

    # Regressione semplice dell'output 'sales' sulla variabile di input 'radio'
    X = df[['radio']]
    regr.fit(X, y)
    y_pred = regr.predict(X)
    plt.scatter(X, y, alpha=0.5)
    plt.plot(X, y_pred, color='red', linewidth=3)
    plt.show()
    print('Intercetta: \n', regr.intercept_, '\nCoefficienti: \n', regr.coef_)
    print('Mean squared error (MSE):\n %.3f' % mean_squared_error(y, y_pred))
    print('Coefficiente R^2:\n %.3f' % r2_score(y, y_pred))

    Intercetta:
     [9.3116381]
    Coefficienti:
     [[0.20249578]]
    Mean squared error (MSE):
     18.092
    Coefficiente R^2:
     0.332

## Modello basato su `newspaper`¶

In [13]:

    # Regressione semplice dell'output 'sales' sulla variabile di input 'newspaper'
    X = df[['newspaper']]
    regr.fit(X, y)
    y_pred = regr.predict(X)
    plt.scatter(X, y, alpha=0.5)
    plt.plot(X, y_pred, color='red', linewidth=3)
    plt.show()
    print('Intercetta: \n', regr.intercept_, '\nCoefficienti: \n', regr.coef_)
    print('Mean squared error (MSE):\n %.3f' % mean_squared_error(y, y_pred))
    print('Coefficiente R^2:\n %.3f' % r2_score(y, y_pred))

    Intercetta:
     [12.35140707]
    Coefficienti:
     [[0.0546931]]
    Mean squared error (MSE):
     25.674
    Coefficiente R^2:
     0.052

## Modello con `TV` e `radio`¶

In [14]:

    # Regressione multipla dell'output 'sales' sulle variabili di input 'TV' e 'radio'
    X = df[['TV', 'radio']]
    regr.fit(X, y)
    y_pred = regr.predict(X)
    print('Intercetta: \n', regr.intercept_, '\nCoefficienti: \n', regr.coef_)
    print('Mean squared error (MSE):\n %.3f' % mean_squared_error(y, y_pred))
    print('Coefficiente R^2:\n %.3f' % r2_score(y, y_pred))

    Intercetta:
     [2.92109991]
    Coefficienti:
     [[0.04575482 0.18799423]]
    Mean squared error (MSE):
     2.785
    Coefficiente R^2:
     0.897

In [15]:

    # Grafico tridimensionale (i dettagli del codice possono essere ignorati)
    from mpl_toolkits.mplot3d import Axes3D
    ax = Axes3D(plt.figure(figsize=(15,10)))
    ax.set_xlabel('TV'); ax.set_ylabel('radio'); ax.set_zlabel('sales');
    ax.scatter(X['TV'], X['radio'], y)
    TICS = 10
    X1 = np.linspace(0.0, 300.0, TICS)
    X2 = np.linspace(0.0, 50.0, TICS)
    X1, X2 = np.meshgrid(X1, X2)
    X_all = np.c_[X1.reshape(TICS**2, 1), X2.reshape(TICS**2, 1)]
    ax.plot_surface(X1, X2, regr.predict(X_all).reshape(TICS, TICS), color='red', alpha=0.3)
    plt.show()

In [16]:

    # Distribuzione (istogramma) dei residui y - y_pred
    (y - y_pred).hist(bins=20);

In [17]:

    (y - y_pred).describe().transpose()

Out[17]:

| count | mean | std | min | 25% | 50% | 75% | max
--- | --- | --- | --- | --- | --- | --- | --- | ---
sales | 200.0 | -1.274536e-15 | 1.672891 | -8.7977 | -0.875158 | 0.242194 | 1.17077 | 2.832837

## Modello completo (3 variabili)¶

In [18]:

    # Regressione multipla dell'output 'sales' su tutte e tre le variabili di input
    X = df[['TV', 'radio', 'newspaper']]
    regr.fit(X, y)
    y_pred = regr.predict(X)
    print('Intercetta: \n', regr.intercept_, '\nCoefficienti: \n', regr.coef_)
    print('Mean squared error (MSE):\n %.3f' % mean_squared_error(y, y_pred))
    print('Coefficiente R^2:\n %.3f' % r2_score(y, y_pred))

    Intercetta:
     [2.93888937]
    Coefficienti:
     [[ 0.04576465  0.18853002 -0.00103749]]
    Mean squared error (MSE):
     2.784
    Coefficiente R^2:
     0.897

## Aggiunta di una interazione nonlineare (feature crossing)¶

In [19]:

    df2 = df.copy()
    df2['TV x radio'] = df['TV'] * df['radio']
    train, test = train_test_split(df2, test_size=0.2, random_state=0) # Split 80%-20%
    X_train = train[['TV', 'radio', 'TV x radio']]
    y_train = train[['sales']]
    X_test = test[['TV', 'radio', 'TV x radio']]
    y_test = test[['sales']]
    regr.fit(X_train, y_train)
    y_pred_train = regr.predict(X_train)
    print('Training R^2:\n %.3f' % r2_score(y_train, y_pred_train))
    print('Training MSE:\n %.3f' % mean_squared_error(y_train, y_pred_train))
    y_pred_test = regr.predict(X_test)
    print('Test MSE:\n %.3f' % mean_squared_error(y_test, y_pred_test))

    Training R^2:
     0.973
    Training MSE:
     0.690
    Test MSE:
     1.644

## Regressione polinomiale¶

In [20]:

    # Regressione polinomiale
    X = df[['TV']]
    y = df[['sales']]
    poly_features = sklearn.preprocessing.PolynomialFeatures(degree=3, include_bias=False)
    X_poly = poly_features.fit_transform(X)
    print(X.iloc[0][0])
    print(X_poly[0])

    230.1
    [2.30100000e+02 5.29460100e+04 1.21828769e+07]

In [21]:

    regr = sklearn.linear_model.LinearRegression()
    regr.fit(X_poly, y);
    # Calcoliamo le predizioni sull'intero range da 0.0 a 40.0
    y_pred = regr.predict(X_poly)
    X_poly_all = poly_features.fit_transform(np.c_[np.arange(0.0, 300.0, 0.5)]) # np.c_ rende il range un array colonna 2D
    y_pred_all = regr.predict(X_poly_all)
    # Grafichiamo i dati originali e l'ipotesi
    plt.scatter(X_poly[:,0], y, alpha=0.5)
    plt.plot(X_poly_all[:,0], y_pred_all, color='red', linewidth=3)
    plt.show()
    print('Intercetta: \n', regr.intercept_, '\nCoefficienti: \n', regr.coef_)

    Intercetta:
     [5.42010655]
    Coefficienti:
     [[ 9.64341770e-02 -3.15222433e-04  5.57199796e-07]]

In [22]:

    # Stampiamo l'errore quadratico medio
    print('Mean squared error (MSE):\n %.3f' % mean_squared_error(y, y_pred))

    Mean squared error (MSE):
     10.238

## Regressione non parametrica (K-Nearest Neighbors)¶

In [23]:

    X_train = train[['TV']]
    y_train = train[['sales']]
    knn = sklearn.neighbors.KNeighborsRegressor(n_neighbors=21) # il parametro K è controllato da n_neighbors
    knn.fit(X_train, y_train)
    y_pred_train = knn.predict(X_train)
    X = np.c_[np.arange(0.0, 300.0, 1.0)]
    y_pred = knn.predict(X)
    # Grafichiamo i dati originali...
    plt.scatter(X_train, y_train, alpha=0.5)
    # ...e le predizioni Nearest Neighbor
    plt.plot(X, y_pred, color='red', linewidth=3)
    plt.show()

In [24]:

    print('Training MSE:\n %.3f' % mean_squared_error(y_train, y_pred_train))
    X_test = test[['TV']]
    y_test = test[['sales']]
    y_pred_test = knn.predict(X_test)
    print('Test MSE:\n %.3f' % mean_squared_error(y_test, y_pred_test))

    Training MSE:
     10.410
    Test MSE:
     11.147

## Regressione lineare tramite Stochastic Gradient Descent¶

In [25]:

    X_train = train[['TV', 'radio', 'TV x radio']]
    y_train = train['sales']    # notare la parentesi singola:
                                # SGDRegressor richiede un array 1D di forma (n_samples,), non una matrice
    X_test = test[['TV', 'radio', 'TV x radio']]
    y_test = test['sales']
    from sklearn.pipeline import make_pipeline
    sgd = make_pipeline(sklearn.preprocessing.StandardScaler(), # preprocessiamo i dati scalandoli
                        sklearn.linear_model.SGDRegressor(loss='squared_loss', penalty='none',
                                                          learning_rate='constant', eta0=0.1, max_iter=2000))
    sgd.fit(X_train, y_train);
    y_pred_train = sgd.predict(X_train)
    print('Training R^2:\n %.3f' % r2_score(y_train, y_pred_train))
    print('Training MSE:\n %.3f' % mean_squared_error(y_train, y_pred_train))
    y_pred_test = sgd.predict(X_test)
    print('Test MSE:\n %.3f' % mean_squared_error(y_test, y_pred_test))

    Training R^2:
     0.970
    Training MSE:
     0.764
    Test MSE:
     1.671