注意
跳转至末尾可下载完整示例代码，或通过 JupyterLite 或 Binder 在浏览器中运行此示例。
梯度提升回归的预测区间#
本示例展示了如何使用分位数回归来创建预测区间。有关 HistGradientBoostingRegressor 其他特性的展示，请参阅直方图梯度提升树的特性。
# Authors: The scikit-learn developers
# SPDX-License-Identifier: BSD-3-Clause
通过将函数 f 应用于均匀采样的随机输入，为合成回归问题生成一些数据。
import numpy as np
from sklearn.model_selection import train_test_split
def f(x):
"""The function to predict."""
return x * np.sin(x)
rng = np.random.RandomState(42)
X = np.atleast_2d(rng.uniform(0, 10.0, size=1000)).T
expected_y = f(X).ravel()
为了使问题更有趣，我们将目标值 y 生成为：由函数 f 计算的确定性项与遵循中心化 对数正态分布 的随机噪声项之和。为了更进一步，我们考虑噪声幅度依赖于输入变量 x 的情况（异方差噪声）。
对数正态分布是非对称且长尾的：观察到大的离群值是有可能的，但不可能观察到小的离群值。
sigma = 0.5 + X.ravel() / 10
noise = rng.lognormal(sigma=sigma) - np.exp(sigma**2 / 2)
y = expected_y + noise
划分为训练集和测试集
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=0)
拟合非线性分位数回归器和最小二乘回归器#
拟合使用分位数损失（quantile loss）训练的梯度提升模型，设置 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)
请注意，在中大规模数据集（n_samples >= 10_000）上，HistGradientBoostingRegressor 的速度比 GradientBoostingRegressor 快得多，但本示例的情况并非如此。
为了进行比较，我们还拟合了一个使用常规均方误差 (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 百分位数回归器。
为了证实这一假设，我们通过在 alpha=0.05 的 Pinball 损失上进行交叉验证来选择最佳模型参数，从而调整一个新的第 5 百分位数回归器的超参数。
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% 覆盖区间，该区间由这两个经过调优的分位数回归器的预测组成。请注意，由于离群值的存在，上 95th 百分位数的预测形状比下 5th 百分位数的预测形状粗糙得多。
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% 的置信水平；这不要与源自 5% 和 95% 分位数预测的 90% 覆盖率混淆。
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 秒)
相关示例
