直方图梯度提升树的特性#

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

HGBT 模型的主要可用性特性是

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

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

  3. 提前停止。

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

  5. 单调约束.

  6. 交互约束.

本示例旨在展示除 2 和 6 之外的所有要点，并在现实生活中进行设置。

# Author: Arturo Amor <[email protected]>
# License: 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 分钟 42.991 秒）

相关示例

时间序列预测的滞后特征

时间序列预测的滞后特征

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

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

单调约束

单调约束

与时间相关的特征工程

与时间相关的特征工程

由 Sphinx-Gallery 生成的图库

本页内容
Launch JupyterLite
Launch binder