Особенности деревьев градиентного бустинга на гистограммах

Подготовка данных

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

Набор данных, первоначально называвшийся ELEC2, содержит 45 312 наблюдений за период с 7 мая 1996 года по 5 декабря 1998 года. Каждое наблюдение относится к 30-минутному периоду, то есть на каждый день приходится 48 наблюдений. Каждое наблюдение в наборе данных содержит 7 столбцов:

• date: период с 7 мая 1996 года по 5 декабря 1998 года. Нормализовано в диапазоне от 0 до 1;
• day: день недели (1–7);
• period: получасовые интервалы в течение суток. Нормализовано в диапазоне от 0 до 1;
• nswprice/nswdemand: цена/спрос на электроэнергию в Новом Южном Уэльсе;
• vicprice/vicdemand: цена/спрос на электроэнергию в Виктории.

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

from sklearn.datasets import fetch_openml

electricity = fetch_openml(
    name="electricity", version=1, as_frame=True, parser="pandas"
)
df = electricity.frame

Для первых 17 760 наблюдений в этом наборе данных целевая переменная имеет ступенчато-постоянное значение:

df["transfer"][:17_760].unique()

array([0.414912, 0.500526])

Удалим эти записи и рассмотрим почасовую передачу электроэнергии в разные дни недели:

import matplotlib.pyplot as plt
import seaborn as sns

df = electricity.frame.iloc[17_760:]
X = df.drop(columns=["transfer", "class"])
y = df["transfer"]

fig, ax = plt.subplots(figsize=(15, 10))
pointplot = sns.lineplot(x=df["period"], y=df["transfer"], hue=df["day"], ax=ax)
handles, labels = ax.get_legend_handles_labels()
ax.set(
    title="Hourly energy transfer for different days of the week",
    xlabel="Normalized time of the day",
    ylabel="Normalized energy transfer",
)
_ = ax.legend(handles, ["Sun", "Mon", "Tue", "Wed", "Thu", "Fri", "Sat"])

Обратите внимание, что передача энергии в выходные дни систематически увеличивается.

Влияние числа деревьев и ранней остановки

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

from sklearn.ensemble import HistGradientBoostingRegressor
from sklearn.model_selection import train_test_split

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.4, shuffle=False)

print(f"Training sample size: {X_train.shape[0]}")
print(f"Test sample size: {X_test.shape[0]}")
print(f"Number of features: {X_train.shape[1]}")

Training sample size: 16531
Test sample size: 11021
Number of features: 7

max_iter_list = [5, 50]
average_week_demand = (
    df.loc[X_test.index].groupby(["day", "period"], observed=False)["transfer"].mean()
)
colors = sns.color_palette("colorblind")
fig, ax = plt.subplots(figsize=(10, 5))
average_week_demand.plot(color=colors[0], label="recorded average", linewidth=2, ax=ax)

for idx, max_iter in enumerate(max_iter_list):
    hgbt = HistGradientBoostingRegressor(
        max_iter=max_iter, categorical_features=None, random_state=42
    )
    hgbt.fit(X_train, y_train)

    y_pred = hgbt.predict(X_test)
    prediction_df = df.loc[X_test.index].copy()
    prediction_df["y_pred"] = y_pred
    average_pred = prediction_df.groupby(["day", "period"], observed=False)[
        "y_pred"
    ].mean()
    average_pred.plot(
        color=colors[idx + 1], label=f"max_iter={max_iter}", linewidth=2, ax=ax
    )

ax.set(
    title="Predicted average energy transfer during the week",
    xticks=[(i + 0.2) * 48 for i in range(7)],
    xticklabels=["Sun", "Mon", "Tue", "Wed", "Thu", "Fri", "Sat"],
    xlabel="Time of the week",
    ylabel="Normalized energy transfer",
)
_ = ax.legend()

Модели HGBT могут сходиться всего за несколько итераций (см. Сравнение моделей случайного леса и градиентного бустинга на гистограммах), то есть дальнейшее добавление деревьев уже не улучшает модель. На рисунке выше 5 итераций недостаточно для получения хороших предсказаний. Уже 50 итераций позволяют добиться хорошего результата.

Слишком высокое значение max_iter может ухудшить качество предсказаний и привести к большим, но предотвратимым затратам вычислительных ресурсов. Поэтому в реализации HGBT в scikit-learn предусмотрена автоматическая стратегия ранней остановки. При её использовании модель выделяет часть тренировочных данных для внутренней проверки (validation_fraction) и прекращает обучение, если оценка на проверочной выборке не улучшается (или ухудшается) в течение n_iter_no_change итераций при заданном допуске (tol).

Обратите внимание на компромисс между learning_rate и max_iter: как правило, предпочтительны меньшие скорости обучения, однако для сходимости к минимуму функции потерь требуется больше итераций; большие скорости обучения обеспечивают более быструю сходимость (требуется меньше итераций/деревьев), но минимум функции потерь при этом выше.

Из-за высокой корреляции между скоростью обучения и числом итераций рекомендуется настраивать скорость обучения вместе со всеми остальными важными гиперпараметрами, обучать HBGT на тренировочной выборке с достаточно большим значением max_iter и определять оптимальное значение max_iter с помощью ранней остановки и явной оценки validation_fraction.

common_params = {
    "max_iter": 1_000,
    "learning_rate": 0.3,
    "validation_fraction": 0.2,
    "random_state": 42,
    "categorical_features": None,
    "scoring": "neg_root_mean_squared_error",
}

hgbt = HistGradientBoostingRegressor(early_stopping=True, **common_params)
hgbt.fit(X_train, y_train)

_, ax = plt.subplots()
plt.plot(-hgbt.validation_score_)
_ = ax.set(
    xlabel="number of iterations",
    ylabel="root mean squared error",
    title=f"Loss of hgbt with early stopping (n_iter={hgbt.n_iter_})",
)

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

import math

common_params["max_iter"] = math.ceil(hgbt.n_iter_ / 100) * 100
common_params["early_stopping"] = False
hgbt = HistGradientBoostingRegressor(**common_params)

Поддержка пропущенных значений

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

В этом примере показано, как модели регрессии HGBT обрабатывают значения, пропущенные полностью случайным образом (MCAR), то есть пропуски не зависят ни от наблюдаемых, ни от ненаблюдаемых данных. Смоделируем такую ситуацию, случайным образом заменив значения в случайно выбранных признаках на значения nan.

import numpy as np

from sklearn.metrics import root_mean_squared_error

rng = np.random.RandomState(42)
first_week = slice(0, 336)  # first week in the test set as 7 * 48 = 336
missing_fraction_list = [0, 0.01, 0.03]


def generate_missing_values(X, missing_fraction):
    total_cells = X.shape[0] * X.shape[1]
    num_missing_cells = int(total_cells * missing_fraction)
    row_indices = rng.choice(X.shape[0], num_missing_cells, replace=True)
    col_indices = rng.choice(X.shape[1], num_missing_cells, replace=True)
    X_missing = X.copy()
    X_missing.iloc[row_indices, col_indices] = np.nan
    return X_missing


fig, ax = plt.subplots(figsize=(12, 6))
ax.plot(y_test.values[first_week], label="Actual transfer")

for missing_fraction in missing_fraction_list:
    X_train_missing = generate_missing_values(X_train, missing_fraction)
    X_test_missing = generate_missing_values(X_test, missing_fraction)
    hgbt.fit(X_train_missing, y_train)
    y_pred = hgbt.predict(X_test_missing[first_week])
    rmse = root_mean_squared_error(y_test[first_week], y_pred)
    ax.plot(
        y_pred[first_week],
        label=f"missing_fraction={missing_fraction}, RMSE={rmse:.3f}",
        alpha=0.5,
    )
ax.set(
    title="Daily energy transfer predictions on data with MCAR values",
    xticks=[(i + 0.2) * 48 for i in range(7)],
    xticklabels=["Mon", "Tue", "Wed", "Thu", "Fri", "Sat", "Sun"],
    xlabel="Time of the week",
    ylabel="Normalized energy transfer",
)
_ = ax.legend(loc="lower right")

Как и ожидалось, качество модели снижается по мере увеличения доли пропущенных значений.

Поддержка квантильной функции потерь

Квантильная функция потерь в регрессии позволяет оценить изменчивость или неопределённость целевой переменной. Например, предсказание 5-го и 95-го процентилей может дать 90%-й предиктивный интервал, то есть диапазон, в который новое наблюдаемое значение, как ожидается, попадёт с вероятностью 90%.

from sklearn.metrics import mean_pinball_loss

quantiles = [0.95, 0.05]
predictions = []

fig, ax = plt.subplots(figsize=(12, 6))
ax.plot(y_test.values[first_week], label="Actual transfer")

for quantile in quantiles:
    hgbt_quantile = HistGradientBoostingRegressor(
        loss="quantile", quantile=quantile, **common_params
    )
    hgbt_quantile.fit(X_train, y_train)
    y_pred = hgbt_quantile.predict(X_test[first_week])

    predictions.append(y_pred)
    score = mean_pinball_loss(y_test[first_week], y_pred)
    ax.plot(
        y_pred[first_week],
        label=f"quantile={quantile}, pinball loss={score:.2f}",
        alpha=0.5,
    )

ax.fill_between(
    range(len(predictions[0][first_week])),
    predictions[0][first_week],
    predictions[1][first_week],
    color=colors[0],
    alpha=0.1,
)
ax.set(
    title="Daily energy transfer predictions with quantile loss",
    xticks=[(i + 0.2) * 48 for i in range(7)],
    xticklabels=["Mon", "Tue", "Wed", "Thu", "Fri", "Sat", "Sun"],
    xlabel="Time of the week",
    ylabel="Normalized energy transfer",
)
_ = ax.legend(loc="lower right")

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

• собирать больше данных;
• лучше настраивать гиперпараметры модели, см. Предиктивные интервалы для регрессии с градиентным бустингом;
• создавать более информативные признаки на основе тех же данных, см. Создание признаков, связанных со временем.

Монотонные ограничения

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

В этом примере передача энергии из Виктории в Новый Южный Уэльс должна сглаживать колебания цен. Это означает, что предсказания модели должны обеспечивать достижение этой цели: объём передачи должен увеличиваться вместе с ценой и спросом в Новом Южном Уэльсе, но уменьшаться вместе с ценой и спросом в Виктории, чтобы это было выгодно населению обоих штатов.

Если тренировочные данные содержат названия признаков, монотонные ограничения можно задать, передав словарь со следующими значениями:

• 1: монотонное возрастание
• 0: без ограничений
• -1: монотонное убывание

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

from sklearn.inspection import PartialDependenceDisplay

monotonic_cst = {
    "date": 0,
    "day": 0,
    "period": 0,
    "nswdemand": 1,
    "nswprice": 1,
    "vicdemand": -1,
    "vicprice": -1,
}
hgbt_no_cst = HistGradientBoostingRegressor(
    categorical_features=None, random_state=42
).fit(X, y)
hgbt_cst = HistGradientBoostingRegressor(
    monotonic_cst=monotonic_cst, categorical_features=None, random_state=42
).fit(X, y)

fig, ax = plt.subplots(nrows=2, figsize=(15, 10))
disp = PartialDependenceDisplay.from_estimator(
    hgbt_no_cst,
    X,
    features=["nswdemand", "nswprice"],
    line_kw={"linewidth": 2, "label": "unconstrained", "color": "tab:blue"},
    ax=ax[0],
)
PartialDependenceDisplay.from_estimator(
    hgbt_cst,
    X,
    features=["nswdemand", "nswprice"],
    line_kw={"linewidth": 2, "label": "constrained", "color": "tab:orange"},
    ax=disp.axes_,
)
disp = PartialDependenceDisplay.from_estimator(
    hgbt_no_cst,
    X,
    features=["vicdemand", "vicprice"],
    line_kw={"linewidth": 2, "label": "unconstrained", "color": "tab:blue"},
    ax=ax[1],
)
PartialDependenceDisplay.from_estimator(
    hgbt_cst,
    X,
    features=["vicdemand", "vicprice"],
    line_kw={"linewidth": 2, "label": "constrained", "color": "tab:orange"},
    ax=disp.axes_,
)
_ = plt.legend()

Обратите внимание, что nswdemand и vicdemand и без ограничений, по-видимому, уже ведут себя монотонно. Это наглядно показывает, что модель с ограничениями монотонности оказывается «излишне ограниченной».

Кроме того, можно убедиться, что введение монотонных ограничений не приводит к существенному ухудшению качества предсказаний модели. Для этого воспользуемся кросс-валидацией TimeSeriesSplit, чтобы оценить дисперсию оценки на тестовой выборке. Это позволяет гарантировать, что тренировочные данные не идут после тестовых, что крайне важно при работе с данными, связанными во времени.

from sklearn.metrics import make_scorer, root_mean_squared_error
from sklearn.model_selection import TimeSeriesSplit, cross_validate

ts_cv = TimeSeriesSplit(n_splits=5, gap=48, test_size=336)  # a week has 336 samples
scorer = make_scorer(root_mean_squared_error)

cv_results = cross_validate(hgbt_no_cst, X, y, cv=ts_cv, scoring=scorer)
rmse = cv_results["test_score"]
print(f"RMSE without constraints = {rmse.mean():.3f} +/- {rmse.std():.3f}")

cv_results = cross_validate(hgbt_cst, X, y, cv=ts_cv, scoring=scorer)
rmse = cv_results["test_score"]
print(f"RMSE with constraints    = {rmse.mean():.3f} +/- {rmse.std():.3f}")

RMSE without constraints = 0.109 +/- 0.037
RMSE with constraints    = 0.108 +/- 0.036

При этом следует учитывать, что сравниваются две разные модели, оптимальные комбинации гиперпараметров которых могут различаться. Поэтому в этом разделе мы не используем common_params, как делали ранее.

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

Download Jupyter notebook: plot_hgbt_regression.ipynb

Download Python source code: plot_hgbt_regression.py

Download zipped: plot_hgbt_regression.zip

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

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

Объединение предикторов с помощью стекинга

Монотонные ограничения

Сравнение моделей случайного леса и градиентного бустинга на гистограммах