注意
转到末尾 下载完整的示例代码。或通过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的观测值生成成分布计算的确定性项和服从中心对数正态分布的随机噪声项的和。为了使问题更具挑战性，我们考虑噪声幅度取决于输入变量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)
拟合非线性分位数和最小二乘回归器#
拟合使用分位数损失和alpha=0.05、0.5、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]范围。
xx = np.atleast_2d(np.linspace(0, 10, 1000)).T
绘制真实的条件均值函数f、条件均值的预测（损失等于平方误差）、条件中位数和条件90%区间（从第5个到第95个条件百分位数）。
import matplotlib.pyplot as plt
y_pred = all_models["mse"].predict(xx)
y_lower = all_models["q 0.05"].predict(xx)
y_upper = all_models["q 0.95"].predict(xx)
y_med = all_models["q 0.50"].predict(xx)
fig = plt.figure(figsize=(10, 10))
plt.plot(xx, f(xx), "g:", linewidth=3, label=r"$f(x) = x\,\sin(x)$")
plt.plot(X_test, y_test, "b.", markersize=10, label="Test observations")
plt.plot(xx, y_med, "r-", label="Predicted median")
plt.plot(xx, y_pred, "r-", label="Predicted mean")
plt.plot(xx, y_upper, "k-")
plt.plot(xx, y_lower, "k-")
plt.fill_between(
xx.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 而言，条件中位数估计器与平方误差估计器具有竞争力：这可以解释为平方误差估计器对导致显著过拟合的大异常值非常敏感。这可以在上一个图的右侧看到。条件中位数估计器是有偏差的（对于这种非对称噪声来说是低估），但它对异常值也具有天然的鲁棒性，并且过拟合程度较低。
置信区间的校准#
我们还可以评估两个极端分位数估计器在产生良好校准的条件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,
all_models["q 0.05"].predict(X_train),
all_models["q 0.95"].predict(X_train),
)
np.float64(0.9)
在训练集上，校准非常接近90%置信区间的预期覆盖值。
coverage_fraction(
y_test, all_models["q 0.05"].predict(X_test), all_models["q 0.95"].predict(X_test)
)
np.float64(0.868)
在测试集上，估计的置信区间略窄。但是，请注意，我们需要将这些指标包装在交叉验证循环中，以评估它们在数据重采样下的变异性。
调整分位数回归器的超参数#
在上图中，我们观察到第5个百分位数回归器似乎拟合不足，无法适应信号的正弦波形状。
模型的超参数大约是为中位数回归器手动调整的，没有理由相同的超参数也适用于第5个百分位数回归器。
为了证实这一假设，我们通过在pinball损失（alpha=0.05）上进行交叉验证来选择最佳模型参数，从而调整新的第5个百分位数回归器的超参数。
from sklearn.experimental import enable_halving_search_cv # noqa
from sklearn.model_selection import HalvingRandomSearchCV
from sklearn.metrics import make_scorer
from pprint import pprint
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],
)
alpha = 0.05
neg_mean_pinball_loss_05p_scorer = make_scorer(
mean_pinball_loss,
alpha=alpha,
greater_is_better=False, # maximize the negative loss
)
gbr = GradientBoostingRegressor(loss="quantile", alpha=alpha, 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
alpha = 0.95
neg_mean_pinball_loss_95p_scorer = make_scorer(
mean_pinball_loss,
alpha=alpha,
greater_is_better=False, # maximize the negative loss
)
search_95p = clone(search_05p).set_params(
estimator__alpha=alpha,
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 = search_05p.predict(xx)
y_upper = search_95p.predict(xx)
fig = plt.figure(figsize=(10, 10))
plt.plot(xx, f(xx), "g:", linewidth=3, label=r"$f(x) = x\,\sin(x)$")
plt.plot(X_test, y_test, "b.", markersize=10, label="Test observations")
plt.plot(xx, y_upper, "k-")
plt.plot(xx, y_lower, "k-")
plt.fill_between(
xx.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.title("Prediction with tuned hyper-parameters")
plt.show()
与未调整的模型相比，该图在视觉上看起来更好，尤其是在较低分位数的形状方面。
现在，我们定量评估估计器对的联合校准。
coverage_fraction(y_train, search_05p.predict(X_train), search_95p.predict(X_train))
np.float64(0.9026666666666666)
coverage_fraction(y_test, search_05p.predict(X_test), search_95p.predict(X_test))
np.float64(0.796)
遗憾的是，调整后的对在测试集上的校准并没有更好：估计的置信区间的宽度仍然太窄。
再次强调，我们需要将此研究包装在交叉验证循环中，以更好地评估这些估计值的变异性。
脚本总运行时间：（0 分钟 12.002 秒）
相关示例
