Note
Go to the end to download the full example code. or to run this example in your browser via Binder
梯度提升回归的预测区间#
本示例展示了如何使用分位数回归来创建预测区间。有关展示
HistGradientBoostingRegressor 其他功能的示例，请参见
直方图梯度提升树的特性 。
生成一些用于合成回归问题的数据，方法是将函数 f 应用于均匀采样的随机输入。
import numpy as np
from sklearn.model_selection import train_test_split
def f(x):
"""要预测的函数。"""
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)
请注意，HistGradientBoostingRegressor 比 GradientBoostingRegressor 快得多，尤其是在处理中等规模的数据集（ n_samples >= 10_000 ）时，而这并不是当前示例的情况。
为了进行比较，我们还拟合了一个使用常规（均值）平方误差（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百分位数回归器。
为了验证这一假设，我们通过在 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_)
{'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)
调优后的配对在测试集上的校准效果仍然不佳：估计的置信区间宽度仍然过窄。
我们需要将这项研究包装在交叉验证循环中，以更好地评估这些估计的变异性。
Total running time of the script: (0 minutes 5.393 seconds)
Related examples
