Comparison of various regression algorithms

_Regressor_comparison.md

Regressor comparison

I wrote a blog post about this here.

I tweeted about it here and again

Inspired, of course, by the various wonderful comparisons in the sklearn docs, like this one for classifiers.

compare_regressors.py

# Author: Matt Hall
# Email: matt@agilescientific.com
# License: BSD 3 clause
from sklearn.linear_model import Ridge, HuberRegressor
from sklearn.preprocessing import PolynomialFeatures
from sklearn.neighbors import KNeighborsRegressor
from sklearn.svm import SVR
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.tree import DecisionTreeRegressor
from sklearn.ensemble import RandomForestRegressor
from sklearn.neural_network import MLPRegressor
from sklearn.pipeline import make_pipeline
from sklearn.metrics import mean_squared_error
from sklearn.model_selection import train_test_split
import numpy as np
import matplotlib.pyplot as plt


def make_linear(N=50, noise=3, random_state=None, w=10, b=3):
    """Make x and y for a 1D linear regression problem."""
    def f(x, w, b):
        return w*x + b
    rng = np.random.default_rng(random_state)
    x = np.linspace(-2.5, 2.5, num=N) + rng.normal(0, noise/10, N)
    y = f(x, w, b) + rng.normal(0, noise, N)
    return x.reshape(-1, 1), y

def make_noisy(N=50, noise=3, random_state=None):
    xa, ya = make_linear(N=N, noise=noise, random_state=random_state)

    rng = np.random.default_rng(random_state)
    
    if noise:
        xn = np.linspace(-2.5, 2.5, num=N//2)
        yn = 70 * (0.5 - rng.random(size=N//2))
        xy = [[x, y] for x, y in zip(xn, yn) if y < 10*x -3]
        xb, yb = np.array(xy).T
    else:
        xb, yb = np.array([]), np.array([])

    return np.vstack([xa, xb.reshape(-1, 1)]), np.hstack([ya, yb])

def make_poly(N=50, noise=3, random_state=None):
    def f(x):
        return 3*x**2 + 9*x - 10
    rng = np.random.default_rng(random_state)
    x = np.linspace(-2.25, 2.25, num=N) + rng.normal(0, noise/10, N)
    y = f(x) + rng.normal(0, noise, N)
    return x.reshape(-1, 1), y

def make_periodic(N=50, noise=3, random_state=None):
    def f(x):
        return 10*np.sin(5*x) + 3*np.cos(3*x) + 5*np.sin(7*x)
    rng = np.random.default_rng(42)
    x = np.linspace(-2.25, 2.25, num=N) + rng.normal(0, noise/10, N)
    y = f(x) + rng.normal(0, noise, N)
    return x.reshape(-1, 1), y

def create_regression_datasets(N=50, noise=3, random_state=None):
    funcs = {
        'Linear': make_linear,
        'Noisy': make_noisy,
        'Polynomial': make_poly,
        'Periodic': make_periodic,
    }
    return {k: f(N, noise, random_state) for k, f in funcs.items()}

N = 60
random_state = 13

models = {
    '': dict(),
    'Linear': dict(model=Ridge(), pen='alpha', mi=0, ma=10),
    'Polynomial': dict(model=make_pipeline(PolynomialFeatures(2), Ridge()), pen='ridge__alpha', mi=0, ma=10),
    'Huber': dict(model=HuberRegressor(), pen='alpha', mi=0, ma=10),
    'Nearest Neighbours': dict(model=KNeighborsRegressor(), pen='n_neighbors', mi=3, ma=9),
    'Linear SVM': dict(model=SVR(kernel='linear'), pen='C', mi=1e6, ma=1),
    'RBF SVM': dict(model=SVR(kernel='rbf'), pen='C', mi=1e6, ma=1),
    'Gaussian Process': dict(model=GaussianProcessRegressor(random_state=random_state), pen='alpha', mi=1e-12, ma=1),
    'Decision Tree': dict(model=DecisionTreeRegressor(random_state=random_state), pen='max_depth', mi=20, ma=3),
    'Random Forest': dict(model=RandomForestRegressor(random_state=random_state), pen='max_depth', mi=20, ma=4),
    'Neural Net': dict(model=MLPRegressor(hidden_layer_sizes=(50, 50), max_iter=1000, tol=0.01, random_state=random_state), pen='alpha', mi=0, ma=10),    
}

datasets = create_regression_datasets(N=N, noise=3, random_state=random_state)
noiseless = create_regression_datasets(N=N, noise=0, random_state=0)

fig, axs = plt.subplots(nrows=len(datasets),
                        ncols=len(models),
                        figsize=(3*len(models), 3*len(datasets)),
                        sharex=True, facecolor='white'
                        )

label_rmse, label_train = True, True

for ax_row, (dataname, (x, y)), (_, (x_, y_)) in zip(axs, datasets.items(), noiseless.items()):
    for ax, (modelname, model) in zip(ax_row, models.items()):
        
        if dataname != 'Noisy':
            x_train, x_val, y_train, y_val = train_test_split(x, y, test_size=0.4, random_state=random_state)
        else:
            x_sig, x_noise = x[:N], x[N:]
            y_sig, y_noise = y[:N], y[N:]
            x_train, x_val, y_train, y_val = train_test_split(x_sig, y_sig, test_size=0.4, random_state=random_state)            
            x_train = np.vstack([x_train, x_noise])
            y_train = np.hstack([y_train, y_noise])
            
        # Plot the noise-free case.
        ax.plot(x_, y_, c='g', alpha=0.5, lw=3)

        # Plot the training and validation data.
        ax.scatter(x_train, y_train, s=18, c='k', alpha=0.67)
        ax.scatter(x_val, y_val, s=18, c='k', alpha=0.33)

        if label_train:
            ax.text(-2.75, 27, f'Train', c='k', alpha=0.67)
            ax.text(-1.75, 27, f'Validate', c='k', alpha=0.33)
            label_train = False

        if (m := model.get('model')) is not None:
            
            xm = np.linspace(-2.5, 2.5).reshape(-1, 1)
            if (pen := model.get('pen')) is not None:
                m.set_params(**{pen: model['mi']})  # Min regularization.
            m.fit(x_train, y_train)
            ŷm = m.predict(xm)
            ax.plot(xm, ŷm, 'r', lw=2, alpha=0.6)

            ŷ = m.predict(x_val)
            mscore = np.sqrt(mean_squared_error(y_val, ŷ))
            ax.text(1.8, -30, f'{mscore:.2f}', c='r')
            
            if label_rmse:
                ax.text(0.75, -32.5, f'RMSE', c='k', fontsize='large')
                label_rmse = False       

        if (pen := model.get('pen')) is not None:
            m.set_params(**{pen: model['ma']})  # Max regularization.
            r = m.fit(x_train, y_train)
            ŷr = r.predict(xm)
            ax.plot(xm, ŷr, 'b', lw=2, alpha=0.6)

            ŷ = r.predict(x_val)
            rscore = np.sqrt(mean_squared_error(y_val, ŷ))
            ax.text(1.8, -35, f'{rscore:.2f}', c='b')
            
        ax.set_ylim(-40, 40)
        ax.set_xlim(-3, 3)

        if ax.get_subplotspec().is_first_row():
            ax.set_title(modelname)

        if ax.get_subplotspec().is_first_col():
            ax.text(-2.75, 32, f'{dataname}', c='k', fontsize='x-large')
        else:
            ax.set_yticklabels([])

        ax.grid(c='k', alpha=0.15)

plt.figtext(0.5, 1.0, 'Regressor Comparison', fontsize='xx-large', color='k', ha='center', va='bottom')
plt.figtext(0.0275, 1.0, 'Underlying model', fontsize='x-large', color='g', ha='left', va='bottom')
plt.figtext(0.4, 1.0, 'Minimal regularization', fontsize='x-large', color='red', ha='center', va='bottom')
plt.figtext(0.6, 1.0, 'Lots of regularization', fontsize='x-large', color='blue', ha='center', va='bottom')
plt.tight_layout()
plt.savefig('comparison.png', dpi=90, bbox_inches = "tight")
plt.show()