Предиктивные интервалы для регрессии с градиентным бустингом

Обучение нелинейных квантильных регрессоров и регрессоров метода наименьших квадратов

Обучим модели градиентного бустинга с квантильной функцией потерь и alpha=0.05, alpha=0.5, alpha=0.95.

Модели, полученные для alpha=0.05 и alpha=0.95, дают интервал с покрытием 90% (95% - 5% = 90%).

Модель, обученная с alpha=0.5, оценивает медиану: в среднем одинаковое число наблюдений целевой переменной должно находиться выше и ниже предсказанных значений.

from sklearn.ensemble import GradientBoostingRegressor
from sklearn.metrics import mean_pinball_loss, mean_squared_error

all_models = {}
common_params = dict(
    learning_rate=0.05,
    n_estimators=200,
    max_depth=2,
    min_samples_leaf=9,
    min_samples_split=9,
)
for alpha in [0.05, 0.5, 0.95]:
    gbr = GradientBoostingRegressor(loss="quantile", alpha=alpha, **common_params)
    all_models["q %1.2f" % alpha] = gbr.fit(X_train, y_train)

Обратите внимание, что HistGradientBoostingRegressor работает намного быстрее, чем GradientBoostingRegressor, начиная с наборов данных среднего размера (n_samples >= 10_000), чего нельзя сказать о данном примере.

Для сравнения мы также обучим базовую модель с обычной функцией потерь — среднеквадратичной ошибкой (MSE).

gbr_ls = GradientBoostingRegressor(loss="squared_error", **common_params)
all_models["mse"] = gbr_ls.fit(X_train, y_train)

Создадим равномерно расположенный набор входных значений для оценки в диапазоне [0, 10].

x_plot = np.atleast_2d(np.linspace(0, 10, 1000)).T

Построим график истинной функции условного среднего f, предсказаний условного среднего (функция потерь равна квадрату ошибки), условной медианы и условного интервала 90% (от 5-го до 95-го условного процентиля).

import matplotlib.pyplot as plt

y_pred = all_models["mse"].predict(x_plot)
y_lower = all_models["q 0.05"].predict(x_plot)
y_upper = all_models["q 0.95"].predict(x_plot)
y_med = all_models["q 0.50"].predict(x_plot)

fig = plt.figure(figsize=(10, 10))
plt.plot(x_plot, f(x_plot), "black", linewidth=3, label=r"$f(x) = x\,\sin(x)$")
plt.plot(X_test, y_test, "b.", markersize=10, label="Test observations")
plt.plot(x_plot, y_med, "tab:orange", linewidth=3, label="Predicted median")
plt.plot(x_plot, y_pred, "tab:green", linewidth=3, label="Predicted mean")
plt.fill_between(
    x_plot.ravel(), y_lower, y_upper, alpha=0.4, label="Predicted 90% interval"
)
plt.xlabel("$x$")
plt.ylabel("$f(x)$")
plt.ylim(-10, 25)
plt.legend(loc="upper left")
plt.show()

Сравнив предсказанную медиану с предсказанным средним, мы видим, что медиана в среднем ниже среднего, поскольку распределение шума смещено в сторону больших значений (крупных выбросов). Оценка медианы также выглядит более гладкой благодаря её естественной устойчивости к выбросам.

Также обратите внимание, что индуктивное смещение деревьев градиентного бустинга, к сожалению, не позволяет нашему квантилю 0.05 полностью передать синусоидальную форму сигнала, особенно вблизи x=8. Настройка гиперпараметров может уменьшить этот эффект, как показано в последней части этой записной книжки.

Анализ метрик ошибки

Оценим модели с помощью метрик mean_squared_error и mean_pinball_loss на обучающем наборе данных.

import pandas as pd


def highlight_min(x):
    x_min = x.min()
    return ["font-weight: bold" if v == x_min else "" for v in x]


results = []
for name, gbr in sorted(all_models.items()):
    metrics = {"model": name}
    y_pred = gbr.predict(X_train)
    for alpha in [0.05, 0.5, 0.95]:
        metrics["pbl=%1.2f" % alpha] = mean_pinball_loss(y_train, y_pred, alpha=alpha)
    metrics["MSE"] = mean_squared_error(y_train, y_pred)
    results.append(metrics)

pd.DataFrame(results).set_index("model").style.apply(highlight_min)

В каждом столбце все модели оцениваются по одной и той же метрике. Минимальное значение в столбце должно быть у модели, обученной и оценённой с помощью одной и той же метрики. Если обучение сошлось, это всегда должно выполняться на обучающей выборке.

Обратите внимание: поскольку распределение целевой переменной асимметрично, условное математическое ожидание и условная медиана существенно различаются, поэтому модель с квадратичной ошибкой не даст хорошую оценку условной медианы, как и наоборот.

Если бы распределение целевой переменной было симметричным и не содержало выбросов (например, при гауссовском шуме), оценки медианы и метода наименьших квадратов давали бы похожие предсказания.

Затем проделаем то же самое на тестовой выборке.

results = []
for name, gbr in sorted(all_models.items()):
    metrics = {"model": name}
    y_pred = gbr.predict(X_test)
    for alpha in [0.05, 0.5, 0.95]:
        metrics["pbl=%1.2f" % alpha] = mean_pinball_loss(y_test, y_pred, alpha=alpha)
    metrics["MSE"] = mean_squared_error(y_test, y_pred)
    results.append(metrics)

pd.DataFrame(results).set_index("model").style.apply(highlight_min)

Ошибки выше, то есть модели немного переобучились на данных. Тем не менее, это показывает, что наилучшая тестовая метрика достигается, когда модель обучается путём минимизации этой же метрики.

Обратите внимание, что оценка условной медианы сопоставима с оценкой методом квадратичной ошибки по MSE на тестовой выборке: это можно объяснить тем, что оценка с квадратичной ошибкой очень чувствительна к крупным выбросам, которые могут привести к значительному переобучению. Это видно в правой части предыдущего графика. Оценка условной медианы смещена (занижена при данном асимметричном шуме), но при этом естественным образом устойчива к выбросам и меньше переобучается.

Настройка гиперпараметров квантильных регрессоров

На графике выше мы увидели, что регрессор 5-го процентиля, по-видимому, недообучен и не может подстроиться под синусоидальную форму сигнала.

Гиперпараметры модели для регрессора медианы были приблизительно настроены вручную, и нет оснований считать, что те же гиперпараметры подойдут для регрессора 5-го процентиля.

Чтобы проверить эту гипотезу, настроим гиперпараметры нового регрессора 5-го процентиля, выбрав наилучшие параметры модели с помощью кросс-валидации по функции потерь Пинболла с alpha=0.05:

from pprint import pprint

from sklearn.experimental import enable_halving_search_cv  # noqa: F401
from sklearn.metrics import make_scorer
from sklearn.model_selection import HalvingRandomSearchCV

param_grid = dict(
    learning_rate=[0.05, 0.1, 0.2],
    max_depth=[2, 5, 10],
    min_samples_leaf=[1, 5, 10, 20],
    min_samples_split=[5, 10, 20, 30, 50],
)
neg_mean_pinball_loss_05p_scorer = make_scorer(
    mean_pinball_loss,
    alpha=0.05,
    greater_is_better=False,  # maximize the negative loss
)
gbr = GradientBoostingRegressor(loss="quantile", alpha=0.05, random_state=0)
search_05p = HalvingRandomSearchCV(
    gbr,
    param_grid,
    resource="n_estimators",
    max_resources=250,
    min_resources=50,
    scoring=neg_mean_pinball_loss_05p_scorer,
    n_jobs=2,
    random_state=0,
).fit(X_train, y_train)
pprint(search_05p.best_params_)

{'learning_rate': 0.2,
 'max_depth': 2,
 'min_samples_leaf': 20,
 'min_samples_split': 10,
 'n_estimators': 150}

Мы видим, что гиперпараметры, настроенные вручную для регрессора медианы, находятся в том же диапазоне, что и гиперпараметры, подходящие для регрессора 5-го процентиля.

Теперь настроим гиперпараметры регрессора 95-го процентиля. Нужно заново определить метрику scoring, используемую для выбора наилучшей модели, и одновременно изменить параметр alpha самого внутреннего оценивателя градиентного бустинга:

from sklearn.base import clone

neg_mean_pinball_loss_95p_scorer = make_scorer(
    mean_pinball_loss,
    alpha=0.95,
    greater_is_better=False,  # maximize the negative loss
)
search_95p = clone(search_05p).set_params(
    estimator__alpha=0.95,
    scoring=neg_mean_pinball_loss_95p_scorer,
)
search_95p.fit(X_train, y_train)
pprint(search_95p.best_params_)

{'learning_rate': 0.05,
 'max_depth': 2,
 'min_samples_leaf': 5,
 'min_samples_split': 20,
 'n_estimators': 150}

Результат показывает, что гиперпараметры регрессора 95-го процентиля, найденные процедурой поиска, находятся примерно в том же диапазоне, что и гиперпараметры, настроенные вручную для регрессора медианы, а также гиперпараметры, найденные процедурой поиска для регрессора 5-го процентиля. Тем не менее, поиск гиперпараметров позволил улучшить интервал покрытия 90%, ограниченный предсказаниями двух настроенных квантильных регрессоров. Обратите внимание, что предсказание верхнего 95-го процентиля имеет гораздо более грубую форму, чем предсказание нижнего 5-го процентиля, из-за выбросов:

y_lower_plot = search_05p.predict(x_plot)
y_upper_plot = search_95p.predict(x_plot)

fig = plt.figure(figsize=(10, 10))
plt.plot(x_plot, f(x_plot), "black", linewidth=3, label=r"$f(x) = x\,\sin(x)$")
plt.plot(X_test, y_test, "b.", markersize=10, label="Test observations")
plt.fill_between(
    x_plot.ravel(),
    y_lower_plot,
    y_upper_plot,
    alpha=0.4,
    label="Predicted 90% interval",
)
plt.xlabel("$x$")
plt.ylabel("$f(x)$")
plt.ylim(-10, 25)
plt.legend(loc="upper left")
plt.title("Prediction with tuned hyper-parameters")
plt.show()

Качественно график выглядит лучше, чем для моделей без настройки, особенно форма нижнего квантиля.

Калибровка доверительного интервала

Мы также можем оценить способность двух крайних оценивателей квантилей формировать хорошо откалиброванные предсказания интервала с покрытием 90% (условно относительно X), то есть в среднем 90% наблюдений должны попадать в этот интервал.

Для этого можно вычислить долю покрытия, то есть долю наблюдений, попадающих в предиктивные интервалы:

def coverage_fraction(y, y_low, y_high):
    return np.mean(np.logical_and(y >= y_low, y <= y_high))


coverage_fraction(y_train, search_05p.predict(X_train), search_95p.predict(X_train))

np.float64(0.9026666666666666)

На обучающей выборке калибровка очень близка к ожидаемому значению покрытия 90%.

coverage_fraction(y_test, search_05p.predict(X_test), search_95p.predict(X_test))

np.float64(0.796)

На тестовой выборке оценённый интервал слишком узок, чтобы охватить 90% тестовых точек, но с учётом разумной статистической неопределённости покрытие всё ещё может соответствовать ожидаемому. Для измерения изменчивости доли покрытия во время предсказания без повторного обучения моделей можно использовать scipy.stats.bootstrap. Для оценённого (полученного бутстрэпом) интервала покрытия используем доверительный уровень 95%; его не следует путать с покрытием 90%, полученным из предсказаний квантилей 5% и 95%:

from scipy.stats import bootstrap

train_coverage_bs = bootstrap(
    (
        y_train,
        search_05p.predict(X_train),
        search_95p.predict(X_train),
    ),
    coverage_fraction,
    paired=True,
    confidence_level=0.95,
    n_resamples=1000,
)
ci = train_coverage_bs.confidence_interval
print(
    f"Training-set coverage lies between {ci.low:.1%} and {ci.high:.1%}, "
    f"based on a 95% bootstrap confidence interval."
)

Training-set coverage lies between 88.1% and 92.3%, based on a 95% bootstrap confidence interval.

Обратите внимание, что интервал включает целевое значение покрытия 90%.

test_coverage_bs = bootstrap(
    (
        y_test,
        search_05p.predict(X_test),
        search_95p.predict(X_test),
    ),
    coverage_fraction,
    paired=True,
    confidence_level=0.95,
    n_resamples=1000,
)
ci = test_coverage_bs.confidence_interval
print(
    f"Test-set coverage lies between {ci.low:.1%} and {ci.high:.1%}, "
    f"based on a 95% bootstrap confidence interval."
)

Test-set coverage lies between 73.8% and 84.4%, based on a 95% bootstrap confidence interval.

К сожалению, оценки квантилей настроенных моделей на тестовой выборке откалиброваны недостаточно хорошо: ширина оценённого доверительного интервала слишком мала даже с учётом её вариаций.

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

Download Jupyter notebook: plot_gradient_boosting_quantile.ipynb

Download Python source code: plot_gradient_boosting_quantile.py

Download zipped: plot_gradient_boosting_quantile.zip

Связанные примеры

Квантильная регрессия

Лаговые признаки для прогнозирования временных рядов

Регрессия с гауссовскими процессами: базовый вводный пример

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