梯度提升回归的预测区间#
此示例展示了如何使用分位数回归创建预测区间。有关展示 HistGradientBoostingRegressor 其他一些功能的示例，请参阅 直方图梯度提升树中的特征。
通过将函数 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)
拟合非线性分位数和最小二乘回归器#
拟合使用分位数损失和 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),
)
0.9
在训练集上，校准值非常接近 90% 置信区间的预期覆盖率值。
coverage_fraction(
y_test, all_models["q 0.05"].predict(X_test), all_models["q 0.95"].predict(X_test)
)
0.868
在测试集上，估计的置信区间略窄。但是请注意，我们需要将这些指标包装在交叉验证循环中，以评估它们在数据重采样下的可变性。
调整分位数回归器的超参数#
在上图中，我们观察到第 5 个百分位数回归器似乎欠拟合，并且无法适应信号的正弦形状。
模型的超参数针对中值回归器进行了近似的手动调整，并且没有理由认为相同的超参数适用于第 5 个百分位数回归器。
为了证实这一假设，我们通过交叉验证选择 alpha=0.05 的 pinball 损失上的最佳模型参数，来调整新的第 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_)
/home/circleci/miniforge3/envs/testenv/lib/python3.9/site-packages/numpy/ma/core.py:2820: RuntimeWarning:
invalid value encountered in cast
{'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))
0.9026666666666666
coverage_fraction(y_test, search_05p.predict(X_test), search_95p.predict(X_test))
0.796
遗憾的是，调整后的对在测试集上的校准并没有更好：估计的置信区间的宽度仍然太窄。
同样，我们需要将这项研究包装在一个交叉验证循环中，以便更好地评估这些估计的可变性。
脚本总运行时间：（0 分 11.444 秒）
相关示例
