注意

转到末尾 下载完整的示例代码。或通过JupyterLite或Binder在浏览器中运行此示例

直方图梯度提升树中的特征#

基于直方图的梯度提升 (HGBT) 模型可能是scikit-learn中最有用的监督学习模型之一。它们基于与LightGBM和XGBoost相当的现代梯度提升实现。因此，HGBT模型比随机森林等替代模型功能更丰富，并且通常性能更好，尤其是在样本数超过数万个时（参见 比较随机森林和直方图梯度提升模型）。

HGBT模型的主要可用性特性包括：

  1. 用于均值和分位数回归任务的几种可用损失函数，参见 分位数损失。

  2. 分类特征支持，参见 梯度提升中的分类特征支持。

  3. 提前停止。

  4. 缺失值支持，避免了对估算器的需求。

  5. 单调约束.

  6. 交互约束.

本示例旨在在一个真实场景中展示除2和6以外的所有点。

# Authors: The scikit-learn developers
# SPDX-License-Identifier: BSD-3-Clause

准备数据#

该 电力数据集 包含从澳大利亚新南威尔士州电力市场收集的数据。在这个市场上，价格不是固定的，受供求关系的影响。它们每五分钟设定一次。与邻近的维多利亚州的电力输送是为了缓解波动。

该数据集最初名为ELEC2，包含从1996年5月7日至1998年12月5日的45,312个实例。数据集的每个样本指的是30分钟的时间段，即每天有48个实例。数据集上的每个样本有7列：

  • 日期：1996年5月7日至1998年12月5日之间。归一化到0到1之间；

  • 日期：星期几（1-7）；

  • 时段：24小时内的半小时间隔。归一化到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, lables = 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"])
Hourly energy transfer for different days of the week

请注意，能源输送在周末会系统地增加。

树木数量和提前停止的影响#

为了说明（最大）树木数量的影响，我们使用整个数据集训练了一个 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()
Predicted average energy transfer during the week

只需几次迭代，HGBT模型就可以达到收敛（参见 比较随机森林和直方图梯度提升模型），这意味着添加更多树木不会再改进模型了。在上图中，5次迭代不足以获得良好的预测。使用50次迭代，我们已经能够做得很好。

设置max_iter过高可能会降低预测质量并浪费大量不必要的计算资源。因此，scikit-learn中的HGBT实现提供了一种自动的**提前停止**策略。使用它，模型会使用一部分训练数据作为内部验证集（validation_fraction），如果验证分数在n_iter_no_change次迭代后在一定容忍度（tol）内没有提高（或降低），则停止训练。

请注意，learning_rate和max_iter之间存在权衡：通常，较小的学习率更可取，但需要更多迭代才能收敛到最小损失，而较大的学习率收敛速度更快（需要较少的迭代/树），但代价是较大的最小损失。

由于学习率与迭代次数之间存在这种高度相关性，一个好的做法是与所有（重要的）其他超参数一起调整学习率，使用足够大的max_iter值在训练集上拟合HBGT，并通过提前停止和一些显式的validation_fraction来确定最佳max_iter。

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_})",
)
Loss of hgbt with early stopping (n_iter=392)

然后，我们可以将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")
Daily energy transfer predictions on data with MCAR values

正如预期的那样，随着缺失值比例的增加，模型性能会下降。

分位数损失支持#

回归中的分位数损失可以展现目标变量的可变性或不确定性。例如，预测第 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")
Daily energy transfer predictions with quantile loss

我们观察到存在高估能量转移的趋势。这可以通过计算经验覆盖率来定量确认，如置信区间校准部分所述。请记住，这些预测的百分位数只是模型的估计值。可以通过以下方法提高此类估计的质量：

单调性约束#

如果特定领域知识要求特征与目标之间的关系必须单调递增或递减，则可以使用单调性约束来强制 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()
plot hgbt regression

观察到 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.103 +/- 0.030
RMSE with constraints    = 0.107 +/- 0.034

话虽如此，需要注意的是，比较的是两个不同的模型，它们可能由不同的超参数组合进行优化。这就是本节中不使用前面使用的 common_params 的原因。

脚本总运行时间：（0 分钟 22.516 秒）

相关示例

时间序列预测的滞后特征

时间序列预测的滞后特征

比较随机森林和直方图梯度提升模型

比较随机森林和直方图梯度提升模型

单调约束

单调约束

时间相关的特征工程

时间相关的特征工程

由 Sphinx-Gallery 生成的图库