Сравнение линейных байесовских регрессоров

Создание синтетического набора данных

Мы создаём набор данных, в котором X и y линейно связаны: 10 признаков из X будут использоваться для генерации y. Остальные признаки бесполезны для предсказания y. Кроме того, мы создаём набор данных, в котором n_samples == n_features. Такая конфигурация представляет сложность для модели МНК и потенциально приводит к произвольно большим весам. Наличие априорного распределения весов и штрафа позволяет решить эту проблему. В завершение добавляется гауссовский шум.

from sklearn.datasets import make_regression

X, y, true_weights = make_regression(
    n_samples=100,
    n_features=100,
    n_informative=10,
    noise=8,
    coef=True,
    random_state=42,
)

Обучение регрессоров

Теперь обучим обе байесовские модели и модель МНК, чтобы затем сравнить коэффициенты этих моделей.

import pandas as pd

from sklearn.linear_model import ARDRegression, BayesianRidge, LinearRegression

olr = LinearRegression().fit(X, y)
brr = BayesianRidge(compute_score=True, max_iter=30).fit(X, y)
ard = ARDRegression(compute_score=True, max_iter=30).fit(X, y)
df = pd.DataFrame(
    {
        "Weights of true generative process": true_weights,
        "ARDRegression": ard.coef_,
        "BayesianRidge": brr.coef_,
        "LinearRegression": olr.coef_,
    }
)

Построение графика истинных и оценённых коэффициентов

Теперь сравним коэффициенты каждой модели с весами истинной порождающей модели.

import matplotlib.pyplot as plt
import seaborn as sns
from matplotlib.colors import SymLogNorm

plt.figure(figsize=(10, 6))
ax = sns.heatmap(
    df.T,
    norm=SymLogNorm(linthresh=10e-4, vmin=-80, vmax=80),
    cbar_kws={"label": "coefficients' values"},
    cmap="seismic_r",
)
plt.ylabel("linear model")
plt.xlabel("coefficients")
plt.tight_layout(rect=(0, 0, 1, 0.95))
_ = plt.title("Models' coefficients")

Из-за добавленного шума ни одной модели не удаётся восстановить истинные веса. Действительно, во всех моделях всегда имеется более 10 ненулевых коэффициентов. По сравнению с оценкой МНК коэффициенты байесовской гребневой регрессии немного смещены к нулю, что делает их более стабильными. Регрессия ARD даёт более разреженное решение: некоторые неинформативные коэффициенты принимают точное значение ноль, а другие смещаются ближе к нулю. Некоторые неинформативные коэффициенты всё же сохраняются и имеют большие значения.

Построение графика предельного логарифмического правдоподобия

import numpy as np

ard_scores = -np.array(ard.scores_)
brr_scores = -np.array(brr.scores_)
plt.plot(ard_scores, color="navy", label="ARD")
plt.plot(brr_scores, color="red", label="BayesianRidge")
plt.ylabel("Log-likelihood")
plt.xlabel("Iterations")
plt.xlim(1, 30)
plt.legend()
_ = plt.title("Models log-likelihood")

Действительно, обе модели минимизируют логарифмическое правдоподобие до произвольного порога, задаваемого параметром max_iter.

Создание синтетического набора данных

Мы создаём целевую переменную, которая является нелинейной функцией входного признака. Добавляется шум со стандартным равномерным распределением.

from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import PolynomialFeatures, StandardScaler

rng = np.random.RandomState(0)
n_samples = 110

# sort the data to make plotting easier later
X = np.sort(-10 * rng.rand(n_samples) + 10)
noise = rng.normal(0, 1, n_samples) * 1.35
y = np.sqrt(X) * np.sin(X) + noise
full_data = pd.DataFrame({"input_feature": X, "target": y})
X = X.reshape((-1, 1))

# extrapolation
X_plot = np.linspace(10, 10.4, 10)
y_plot = np.sqrt(X_plot) * np.sin(X_plot)
X_plot = np.concatenate((X, X_plot.reshape((-1, 1))))
y_plot = np.concatenate((y - noise, y_plot))

Обучение регрессоров

Здесь мы пробуем полином степени 10, чтобы потенциально переобучить модель, хотя байесовские линейные модели регуляризуют величину полиномиальных коэффициентов. Поскольку fit_intercept=True по умолчанию для ARDRegression и BayesianRidge, то PolynomialFeatures не должна добавлять дополнительный признак смещения. Установив return_std=True, байесовские регрессоры возвращают стандартное отклонение апостериорного распределения параметров модели.

ard_poly = make_pipeline(
    PolynomialFeatures(degree=10, include_bias=False),
    StandardScaler(),
    ARDRegression(),
).fit(X, y)
brr_poly = make_pipeline(
    PolynomialFeatures(degree=10, include_bias=False),
    StandardScaler(),
    BayesianRidge(),
).fit(X, y)

y_ard, y_ard_std = ard_poly.predict(X_plot, return_std=True)
y_brr, y_brr_std = brr_poly.predict(X_plot, return_std=True)

Построение графиков полиномиальных регрессий со стандартными ошибками оценок

ax = sns.scatterplot(
    data=full_data, x="input_feature", y="target", color="black", alpha=0.75
)
ax.plot(X_plot, y_plot, color="black", label="Ground Truth")
ax.plot(X_plot, y_brr, color="red", label="BayesianRidge with polynomial features")
ax.plot(X_plot, y_ard, color="navy", label="ARD with polynomial features")
ax.fill_between(
    X_plot.ravel(),
    y_ard - y_ard_std,
    y_ard + y_ard_std,
    color="navy",
    alpha=0.3,
)
ax.fill_between(
    X_plot.ravel(),
    y_brr - y_brr_std,
    y_brr + y_brr_std,
    color="red",
    alpha=0.3,
)
ax.legend()
_ = ax.set_title("Polynomial fit of a non-linear feature")

Планки погрешностей соответствуют одному стандартному отклонению предсказанного гауссовского распределения для точек запроса. Обратите внимание, что при использовании параметров по умолчанию для обеих моделей регрессия ARD лучше всего воспроизводит истинную зависимость, однако дальнейшее уменьшение гиперпараметра lambda_init байесовской гребневой регрессии может снизить её смещение (см. пример Аппроксимация кривой с помощью байесовской гребневой регрессии). Наконец, из-за присущих полиномиальной регрессии ограничений обе модели не справляются с экстраполяцией.

Общее время выполнения скрипта: (0 минут 0.639 секунды)

Download Jupyter notebook: plot_ard.ipynb

Download Python source code: plot_ard.py

Download zipped: plot_ard.zip

Похожие примеры

Модели на основе L1 для разреженных сигналов

Полиномиальная и сплайновая интерполяция

SVM: взвешенные образцы

Коэффициенты гребневой регрессии в зависимости от регуляризации L2