Sklearn-源码解析-书-v1-0-三十一-

Sklearn 源码解析（书）v1.0（三十一）

代码作用：此代码展示了 HistGradientBoostingRegressor 在处理类别特征时的多种策略及其性能表现。通过比较删除、独热编码、序数编码、目标编码和原生类别支持五种方法，我们可以看到：删除类别特征导致性能最差；独热编码虽然能保留信息但显著增加特征维度和训练时间；序数编码假设有序关系可能引入偏 bias；目标编码虽然有效但需防止数据泄漏；而原生类别支持则无需显式编码，直接在树的分裂过程中基于目标统计量对类别进行划分，既保留了信息又避免了高基数导致的维度爆炸，在训练速度和预测精度之间取得了良好平衡。这验证了 HistGradientBoosting 在实际应用中处理类别特征的优越性。

流程图

flowchart TD A[加载 Ames Housing 数据] --> B[选择特征子集] B --> C[将类别特征转为 category 类型] C --> D[创建五种预处理管道] D --> E[交叉验证评估模型性能] E --> F[绘制训练时间 vs 误差散点图]

架构图

flowchart LR A[原始特征 X] --> B[HistGradientBoostingRegressor] B --> C{是否启用原生类别支持} C -->|是| D[直接处理类别特征] C -->|否| E[预处理：独热/序数/目标编码/删除] D --> F[树分裂时基于目标统计量划分类别] E --> F F --> G[最终预测]

设计取舍

问：为什么原生类别支持在高基数特征上比独热编码更高效？

答：独热编码会为每个类别值创建一个新特征，导致特征维度随基数线性增长，增加了树的分裂候选和计算开销。而原生类别支持无需展开类别，直接在原始特征上基于目标统计量（如类别平均目标值）进行划分，避免了维度爆炸，尤其在类别数众多时优势显著。

问：目标编码虽然有效，为何在某些场景下不如原生类别支持推荐？

答：目标编码需要计算每个类别的目标均值，若未使用交叉验证（target encoding with cross-fitting）则可能导致目标泄漏，使模型在训练时看到测试信息，从而在训练集上表现过佳但在测试集上泛化 poorly。即使使用交叉验证，目标编码仍增加了预处理步骤和计算开销，而原生类别支持一步到位，无需额外建模。

70.8 梯度提升早停机制 —— 自动迭代数选择

以下代码演示了如何在梯度提升回归中使用早停机制来自动选择最优的迭代次数，从而在防止过拟合的同时提升训练效率。我们使用加州房价数据集，训练两个 GradientBoostingRegressor 模型：一个不使用早停（固定迭代 1000 次），另一个启用早停，通过在训练过程中监控验证集损失来决定何时停止。通过比较训练时间、使用的估计器数量以及训练/验证损失的收敛行为，我们可以看到早停如何在不显著牺牲性能的前提下大幅减少不必要的迭代，从而提升模型训练效率。

# 第 70 章 —— Authors: The scikit-learn developers
# 第 70 章 —— SPDX-License-Identifier: BSD-3-Clause

# 第 70 章 —— %%
# 第 70 章 —— Data Preparation
# 第 70 章 —— ----------------
# 第 70 章 —— First we load and prepares the California Housing Prices dataset for
# 第 70 章 —— training and evaluation. It subsets the dataset, splits it into training
# 第 70 章 —— and validation sets.

import time

import matplotlib.pyplot as plt

from sklearn.datasets import fetch_california_housing
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.metrics import mean_squared_error
from sklearn.model_selection import train_test_split

data = fetch_california_housing()
X, y = data.data[:600], data.target[:600]

X_train, X_val, y_train, y_val = train_test_split(X, y, test_size=0.2, random_state=42)

# 第 70 章 —— %%
# 第 70 章 —— Model Training and Comparison
# 第 70 章 —— -----------------------------
# 第 70 章 —— Two :class:`~sklearn.ensemble.GradientBoostingRegressor` models are trained:
# 第 70 章 —— one with and another without early stopping. The purpose is to compare their
# 第 70 章 —— performance. It also calculates the training time and the `n_estimators_`
# 第 70 章 —— used by both models.

params = dict(n_estimators=1000, max_depth=5, learning_rate=0.1, random_state=42)

gbm_full = GradientBoostingRegressor(**params)
gbm_early_stopping = GradientBoostingRegressor(
    **params,
    validation_fraction=0.1,
    n_iter_no_change=10,
)

start_time = time.time()
gbm_full.fit(X_train, y_train)
training_time_full = time.time() - start_time
n_estimators_full = gbm_full.n_estimators_

start_time = time.time()
gbm_early_stopping.fit(X_train, y_train)
training_time_early_stopping = time.time() - start_time
estimators_early_stopping = gbm_early_stopping.n_estimators_

# 第 70 章 —— %%
# 第 70 章 —— Error Calculation
# 第 70 章 —— -----------------
# 第 70 章 —— The code calculates the :func:`~sklearn.metrics.mean_squared_error` for both
# 第 70 章 —— training and validation datasets for the models trained in the previous
# 第 70 章 —— section. It computes the errors for each boosting iteration. The purpose is
# 第 70 章 —— to assess the performance and convergence of the models.

train_errors_without = []
val_errors_without = []

train_errors_with = []
val_errors_with = []

for i, (train_pred, val_pred) in enumerate(
    zip(
        gbm_full.staged_predict(X_train),
        gbm_full.staged_predict(X_val),
    )
):
    train_errors_without.append(mean_squared_error(y_train, train_pred))
    val_errors_without.append(mean_squared_error(y_val, val_pred))

for i, (train_pred, val_pred) in enumerate(
    zip(
        gbm_early_stopping.staged_predict(X_train),
        gbm_early_stopping.staged_predict(X_val),
    )
):
    train_errors_with.append(mean_squared_error(y_train, train_pred))
    val_errors_with.append(mean_squared_error(y_val, val_pred))

# 第 70 章 —— %%
# 第 70 章 —— Visualize Comparison
# 第 70 章 —— --------------------
# 第 70 章 —— It includes three subplots:
#
# 第 70 章 —— 1. Plotting training errors of both models over boosting iterations.
#
# 第 70 章 —— 2. Plotting validation errors of both models over boosting iterations.
#
# 第 70 章 —— 3. Creating a bar chart to compare the training times and the estimator used
# 第 70 章 —— of the models with and without early stopping.
#

fig, axes = plt.subplots(ncols=3, figsize=(12, 4))

axes[0].plot(train_errors_without, label="gbm_full")
axes[0].plot(train_errors_with, label="gbm_early_stopping")
axes[0].set_xlabel("Boosting Iterations")
axes[0].set_ylabel("MSE (Training)")
axes[0].set_yscale("log")
axes[0].legend()
axes[0].set_title("Training Error")

axes[1].plot(val_errors_without, label="gbm_full")
axes[1].plot(val_errors_with, label="gbm_early_stopping")
axes[1].set_xlabel("Boosting Iterations")
axes[1].set_ylabel("MSE (Validation)")
axes[1].set_yscale("log")
axes[1].legend()
axes[1].set_title("Validation Error")

training_times = [training_time_full, training_time_early_stopping]
labels = ["gbm_full", "gbm_early_stopping"]
bars = axes[2].bar(labels, training_times)
axes[2].set_ylabel("Training Time (s)")

for bar, n_estimators in zip(bars, [n_estimators_full, estimators_early_stopping]):
    height = bar.get_height()
    axes[2].text(
        bar.get_x() + bar.get_width() / 2,
        height + 0.001,
        f"Estimators: {n_estimators}",
        ha="center",
        va="bottom",
    )

plt.tight_layout()
plt.show()

# 第 70 章 —— %%
# 第 70 章 —— The difference in training error between the `gbm_full` and the
# 第 70 章 —— `gbm_early_stopping` stems from the fact that `gbm_early_stopping` sets
# 第 70 章 —— aside `validation_fraction` of the training data as internal validation set.
# 第 70 章 —— Early stopping is decided based on this internal validation score.

代码作用：此代码展示了早停机制在梯度提升中的实际应用。通过比较带早停和不带早停的两个模型，我们可以看到：启用早停后，模型在验证集损失不再改进时自动停止迭代，从而显著减少了使用的树的数量（ estimators），同时训练时间也大幅降低。尽管训练误差略有增加，但验证误差的表现表明早停成功防止了过拟合，使得模型具有更好的泛化能力。这验证了早停作为一种正则化手段的有效性——它在模型复杂度和训练效率之间取得了良好平衡。

流程图

flowchart TD A[加载加州房价数据] --> B[划分训练/验证集] B --> C[训练完整梯度提升模型] B --> D[训练带早停的梯度提升模型] C --> E[记录训练时间和估计器数量] D --> E E --> F[计算训练/验证误差随迭代变化] F --> G[绘制误差曲线和训练时间柱状图]

架构图

flowchart LR A[原始特征 X] --> B[梯度提弱学习器 DecisionTreeRegressor] B --> C[迭代训练弱学习器] C --> D[计算残差] D --> E[更新弱学习器权重] E --> F[集成预测] F --> G{验证集损失是否改善？} G -->|是| H[继续迭代] G -->|否| I[早停训练] I --> J[返回最终模型]

设计取舍

问：为什么早停需要验证集，而不能仅用训练误差来判断停止时机？

答：训练误差在 boosting 过程中通常会持续下降甚至归零，因为模型能够越拟合越好。仅凭训练误差无法检测过拟合，而验证集上的误差能够反映模型对未见数据的泛化能力，当其不再改善时表明模型开始记住训练噪声。

问：参数 n_iter_no_change 和 validation_fraction 如何共同影响早停行为？

答：validation_fraction 决定用作内部验证集的训练数据比例，越大则验证更可靠但可用于训练的数据越少；n_iter_no_change 设定容忍验证误差无改善的迭代次数，越大则模型更可能继续训练以追求微小提升，但也增加过拟合风险。二者共同控制早停的敏感度和稳健性。

70.9 梯度提升 OOB 估计 —— Stochastic Gradient Boosting 的免费验证

以下代码演示了如何使用随机梯度提升（Stochastic Gradient Boosting）中的袋外（OOB）估计来监控模型在训练过程中的泛化性能，无需额外的验证集。我们构造了一个二分类数据集，训练了一个启用了 subsample（即每次迭代仅使用部分样本）的 GradientBoostingClassifier，并利用其 oob_improvement_ 属性来计算累积的 OOB 损失。随后，我们将 OOB 损失与实际测试集损失和交叉验证损失进行对比，以展示 OOB 估计如何能够有效地追踪模型在早期迭代中的真实泛化趋势，尽管在后期可能变得过于悲观。这验证了 OOB 作为一种“无损”模型评估方法的实用性，尤其在计算资源受限或希望避免显式划分验证集时。

# 第 70 章 —— Authors: The scikit-learn developers
# 第 70 章 —— SPDX-License-Identifier: BSD-3-Clause

import matplotlib.pyplot as plt
import numpy as np
from scipy.special import expit

from sklearn import ensemble
from sklearn.metrics import log_loss
from sklearn.model_selection import KFold, train_test_split

# 第 70 章 —— Generate data (adapted from G. Ridgeway's gbm example)
n_samples = 1000
random_state = np.random.RandomState(13)
x1 = random_state.uniform(size=n_samples)
x2 = random_state.uniform(size=n_samples)
x3 = random_state.randint(0, 4, size=n_samples)

p = expit(np.sin(3 * x1) - 4 * x2 + x3)
y = random_state.binomial(1, p, size=n_samples)

X = np.c_[x1, x2, x3]

X = X.astype(np.float32)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.5, random_state=9)

# 第 70 章 —— Fit classifier with out-of-bag estimates
params = {
    "n_estimators": 1200,
    "max_depth": 3,
    "subsample": 0.5,
    "learning_rate": 0.01,
    "min_samples_leaf": 1,
    "random_state": 3,
}
clf = ensemble.GradientBoostingClassifier(**params)

clf.fit(X_train, y_train)
acc = clf.score(X_test, y_test)
print("Accuracy: {:.4f}".format(acc))

n_estimators = params["n_estimators"]
x = np.arange(n_estimators) + 1


def heldout_score(clf, X_test, y_test):
    """compute deviance scores on ``X_test`` and ``y_test``."""
    score = np.zeros((n_estimators,), dtype=np.float64)
    for i, y_proba in enumerate(clf.staged_predict_proba(X_test)):
        score[i] = 2 * log_loss(y_test, y_proba[:, 1])
    return score


def cv_estimate(n_splits=None):
    cv = KFold(n_splits=n_splits)
    cv_clf = ensemble.GradientBoostingClassifier(**params)
    val_scores = np.zeros((n_estimators,), dtype=np.float64)
    for train, test in cv.split(X_train, y_train):
        cv_clf.fit(X_train[train], y_train[train])
        val_scores += heldout_score(cv_clf, X_train[test], y_train[test])
    val_scores /= n_splits
    return val_scores


# 第 70 章 —— Estimate best n_estimator using cross-validation
cv_score = cv_estimate(3)

# 第 70 章 —— Compute best n_estimator for test data
test_score = heldout_score(clf, X_test, y_test)

# 第 70 章 —— negative cumulative sum of oob improvements
cumsum = -np.cumsum(clf.oob_improvement_)

# 第 70 章 —— min loss according to OOB
oob_best_iter = x[np.argmin(cumsum)]

# 第 70 章 —— min loss according to test (normalize such that first loss is 0)
test_score -= test_score[0]
test_best_iter = x[np.argmin(test_score)]

# 第 70 章 —— min loss according to cv (normalize such that first loss is 0)
cv_score -= cv_score[0]
cv_best_iter = x[np.argmin(cv_score)]

# 第 70 章 —— color brew for the three curves
oob_color = list(map(lambda x: x / 256.0, (190, 174, 212)))
test_color = list(map(lambda x: x / 256.0, (127, 201, 127)))
cv_color = list(map(lambda x: x / 256.0, (253, 192, 134)))

# 第 70 章 —— line type for the three curves
oob_line = "dashed"
test_line = "solid"
cv_line = "dashdot"

# 第 70 章 —— plot curves and vertical lines for best iterations
plt.figure(figsize=(8, 4.8))
plt.plot(x, cumsum, label="OOB loss", color=oob_color, linestyle=oob_line)
plt.plot(x, test_score, label="Test loss", color=test_color, linestyle=test_line)
plt.plot(x, cv_score, label="CV loss", color=cv_color, linestyle=cv_line)
plt.axvline(x=oob_best_iter, color=oob_color, linestyle=oob_line)
plt.axvline(x=test_best_iter, color=test_color, linestyle=test_line)
plt.axvline(x=cv_best_iter, color=cv_color, linestyle=cv_line)

# 第 70 章 —— add three vertical lines to xticks
xticks = plt.xticks()
xticks_pos = np.array(
    xticks[0].tolist() + [oob_best_iter, cv_best_iter, test_best_iter]
)
xticks_label = np.array(list(map(lambda t: int(t), xticks[0])) + ["OOB", "CV", "Test"])
ind = np.argsort(xticks_pos)
xticks_pos = xticks_pos[ind]
xticks_label = xticks_label[ind]
plt.xticks(xticks_pos, xticks_label, rotation=90)

plt.legend(loc="upper center")
plt.ylabel("normalized loss")
plt.xlabel("number of iterations")

plt.show()

代码作用：此代码展示了如何利用随机梯度提升中的袋外（OOB）估计来监控模型训练过程中的泛化性能。通过绘制 OOB 损失、测试集损失和交叉验证损失随迭代次数的变化，我们可以看到 OOB 损失在前几百次迭代中与测试损失紧密追踪，随后开始偏离（变得更 pessimistic），这与理论一致。同时，图中垂直线标注了根据三种估计方法选择的最优迭代次数，显示尽管存在偏差，OOB 仍能提供一个有用的启发式来提前停止训练。这验证了 OOB 作为交叉验证的近似替代方案的价值——它无需额外模模型训练即可在训练过程中提供反馈，尤其在计算开销敏感的场景下非常有用。

流程图

flowchart TD A[生成二分类数据] --> B[划分训练/测试集] B --> C[训练 Stochastic Gradient Boosting 分类器] C --> D[计算 OOB 改进值] D --> E[累积求和得到 OOB 损失] E --> F[计算测试集和交叉验证损失] F --> G[绘制三条损失曲线及最优迭代垂直线]

架构图

flowchart LR A[原始特征 X] --> B[梯度提升弱学习器] B --> C[每次迭代使用子样本（subsample < 1.0）] C --> D[计算袋外样本的损失改进] D --> E[累积负 OOB 改进得到 OOB 损失估计] E --> F[与测试集损失和交叉验证损失对比] F --> G[判断是否早停]

设计取舍

问：为什么 OOB 估计在后期迭代中会变得过于悲观（pessimistic）？

答：随着迭代增加，模型在训练样本上的拟合程度提升，导致袋外样本（未被抽中的样本）在训练中得到的更新减少，使得它们的预测更依赖早期模型，而早期模型欠拟合，因此 OOB 损失会系统性地高估真实泛化误差。

问：OOB 估计相比交叉验证的主要优势是什么？

答：OOB 估计无需额外训练多个模型，能够在单次训练过程中“免费”获得泛化性能的指示，尤其在计算资源有限或希望快速迭代超参数时非常有用，尽管其在后期可能偏悲观。

70.10 梯度提升分位数回归 —— 预测区间与 Pinball Loss 调优

以下代码演示了如何使用梯度提升进行分位数回归，以构建预测区间并评估其校准性。我们首先生成一个带有异方差噪声的非线性合成数据集，其中真实函数为 f(x) = x·sin(x)，噪声幅度随 x 增大。然后训练三个梯度提升回归器，分别使用 quantile loss 和 alpha=0.05、0.5、0.95，其中 alpha=0.05 和 0.95 的组合可形成 90% 的预测区间。我们可视化了真实函数、测试观测值、预测中位数、均值以及该区间，并计算了 pinball loss 和均方误差来评估模型在训练集和测试集上的表现。最后，我们通过调整超参数（如学习率、树深度、叶子最小样本数等）来改进分位数估计，并重新评估预测区间的校准性。这验证了分位数回归在不假设噪声分布形状的情况下提供有用不确定性估计的能力，尽管校准可能受模型表达限制和超参数选择的影响。

# 第 70 章 —— Authors: The scikit-learn developers
# 第 70 章 —— SPDX-License-Identifier: BSD-3-Clause

# 第 70 章 —— %%
# 第 70 章 —— Generate some data for a synthetic regression problem by applying the
# 第 70 章 —— function f to uniformly sampled random inputs.
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()

# 第 70 章 —— %%
# 第 70 章 —— To make the problem interesting, we generate observations of the target y as
# 第 70 章 —— the sum of a deterministic term computed by the function f and a random noise
# 第 70 章 —— term that follows a centered `log-normal
# 第 70 章 —— <https://en.wikipedia.org/wiki/Log-normal_distribution>`_. To make this even
# 第 70 章 —— more interesting we consider the case where the amplitude of the noise
# 第 70 章 —— depends on the input variable x (heteroscedastic noise).
#
# 第 70 章 —— The lognormal distribution is non-symmetric and long tailed: observing large
# 第 70 章 —— outliers is likely but it is impossible to observe small outliers.
sigma = 0.5 + X.ravel() / 10
noise = rng.lognormal(sigma=sigma) - np.exp(sigma**2 / 2)
y = expected_y + noise

# 第 70 章 —— %%
# 第 70 章 —— Split into train, test datasets:
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=0)

# 第 70 章 —— %%
# 第 70 章 —— Fitting non-linear quantile and least squares regressors
# 第 70 章 —— --------------------------------------------------------
#
# 第 70 章 —— Fit gradient boosting models trained with the quantile loss and `alpha=0.05`,
# 第 70 章 —— `alpha=0.5`, `alpha=0.95`.
#
# 第 70 章 —— The models obtained for `alpha=0.05` and `alpha=0.95` produce a 90%
# 第 70 章 —— confidence interval (95% - 5% = 90%).
#
# 第 70 章 —— The model trained with `alpha=0.5` produces a regression of the median: on
# 第 70 章 —— average, there should be the same number of target observations above and
# 第 70 章 —— below the predicted values.
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)

# 第 70 章 —— %%
# 第 70 章 —— Notice that :class:`~sklearn.ensemble.HistGradientBoostingRegressor` is much
# 第 70 章 —— faster than :class:`~sklearn.ensemble.GradientBoostingRegressor` starting with
# 第 70 章 —— intermediate datasets (`n_samples >= 10_000`), which is not the case of the
# 第 70 章 —— present example.
#
# 第 70 章 —— For the sake of comparison, we also fit a baseline model trained with the
# 第 70 章 —— usual (mean) squared error (MSE).
gbr_ls = GradientBoostingRegressor(loss="squared_error", **common_params)
all_models["mse"] = gbr_ls.fit(X_train, y_train)

# 第 70 章 —— %%
# 第 70 章 —— Create an evenly spaced evaluation set of input values spanning the [0, 10]
# 第 70 章 —— range.
xx = np.atleast_2d(np.linspace(0, 10, 1000)).T

# 第 70 章 —— %%
# 第 70 章 —— Plot the true conditional mean function f, the predictions of the conditional
# 第 70 章 —— mean (loss equals squared error), the conditional median and the conditional
# 第 70 章 —— 90% interval (from 5th to 95th conditional percentiles).
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), "black", 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, "tab:orange", linewidth=3, label="Predicted median")
plt.plot(xx, y_pred, "tab:green", linewidth=3, label="Predicted mean")
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()

# 第 70 章 —— %%
# 第 70 章 —— Comparing the predicted median with the predicted mean, we note that the
# 第 70 章 —— median is on average below the mean as the noise is skewed towards high
# 第 70 章 —— values (large outliers). The median estimate also seems to be smoother
# 第 70 章 —— because of its natural robustness to outliers.
#
# 第 70 章 —— Also observe that the inductive bias of gradient boosting trees is
# 第 70 章 —— unfortunately preventing our 0.05 quantile to fully capture the sinoisoidal
# 第 70 章 —— shape of the signal, in particular around x=8. Tuning hyper-parameters can
# 第 70 章 —— reduce this effect as shown in the last part of this notebook.
#
# 第 70 章 —— Analysis of the error metrics
# 第 70 章 —— -----------------------------
#
# 第 70 章 —— Measure the models with :func:`~sklearn.metrics.mean_squared_error` and
# 第 70 章 —— :func:`~sklearn.metrics.mean_pinball_loss` metrics on the training dataset.
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)

# 第 70 章 —— %%
# 第 70 章 —— One column shows all models evaluated by the same metric. The minimum number
# 第 70 章 —— on a column should be obtained when the model is trained and measured with
# 第 70 章 —— the same metric. This should be always the case on the training set if the
# 第 70 章 —— training converged.
#
# 第 70 章 —— Note that because the target distribution is asymmetric, the expected
# 第 70 章 —— conditional mean and conditional median are significantly different and
# 第 70 章 —— therefore one could not use the squared error model get a good estimation of
# 第 70 章 —— the conditional median nor the converse.
#
# 第 70 章 —— If the target distribution were symmetric and had no outliers (e.g. with a
# 第 70 章 —— Gaussian noise), then median estimator and the least squares estimator would
# 第 70 章 —— have yielded similar predictions.
#
# 第 70 章 —— We then do the same on the test set.
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)


# 第 70 章 —— %%
# 第 70 章 —— Errors are higher meaning the models slightly overfitted the data. It still
# 第 70 章 —— shows that the best test metric is obtained when the model is trained by
# 第 70 章 —— minimizing this same metric.
#
# 第 70 章 —— Note that the conditional median estimator is competitive with the squared
# 第 70 章 —— error estimator in terms of MSE on the test set: this can be explained by
# 第 70 章 —— the fact the squared error estimator is very sensitive to large outliers
# 第 70 章 —— which can cause significant overfitting. This can be seen on the right hand
# 第 70 章 —— side of the previous plot. The conditional median estimator is biased
# 第 70 章 —— (underestimation for this asymmetric noise) but is also naturally robust to
# 第 70 章 —— outliers and overfits less.
#
# 第 70 章 —— .. _calibration-section:
#
# 第 70 章 —— Calibration of the confidence interval
# 第 70 章 —— --------------------------------------
#
# 第 70 章 —— We can also evaluate the ability of the two extreme quantile estimators at
# 第 70 章 —— producing a well-calibrated conditional 90%-confidence interval.
#
# 第 70 章 —— To do this we can compute the fraction of observations that fall between the
# 第 70 章 —— predictions:
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),
)

# 第 70 章 —— %%
# 第 70 章 —— On the training set the calibration is very close to the expected coverage
# 第 70 章 —— value for a 90% confidence interval.
coverage_fraction(
    y_test, all_models["q 0.05"].predict(X_test), all_models["q 0.95"].predict(X_test)
)


# 第 70 章 —— %%
# 第 70 章 —— On the test set, the estimated confidence interval is slightly too narrow.
# 第 70 章 —— Note, however, that we would need to wrap those metrics in a cross-validation
# 第 70 章 —— loop to assess their variability under data resampling.
#
# 第 70 章 —— Tuning the hyper-parameters of the quantile regressors
# 第 70 章 —— ------------------------------------------------------
#
# 第 70 章 —— In the plot above, we observed that the 5th percentile regressor seems to
# 第 70 章 —— underfit and could not adapt to sinusoidal shape of the signal.
#
# 第 70 章 —— The hyper-parameters of the model were approximately hand-tuned for the
# 第 70 章 —— median regressor and there is no reason that the same hyper-parameters are
# 第 70 章 —— suitable for the 5th percentile regressor.
#
# 第 70 章 —— To confirm this hypothesis, we tune the hyper-parameters of a new regressor
# 第 70 章 —— of the 5th percentile by selecting the best model parameters by
# 第 70 章 —— cross-validation on the pinball loss with alpha=0.05:

# 第 70 章 —— %%
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],
)
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_)

# 第 70 章 —— %%
# 第 70 章 —— We observe that the hyper-parameters that were hand-tuned for the median
# 第 70 章 —— regressor are in the same range as the hyper-parameters suitable for the 5th
# 第 70 章 —— percentile regressor.
#
# 第 70 章 —— Let's now tune the hyper-parameters for the 95th percentile regressor. We
# 第 70 章 —— need to redefine the `scoring` metric used to select the best model, along
# 第 70 章 —— with adjusting the alpha parameter of the inner gradient boosting estimator
# 第 70 章 —— itself:
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_)

# 第 70 章 —— %%
# 第 70 章 —— The result shows that the hyper-parameters for the 95th percentile regressor
# 第 70 章 —— identified by the search procedure are roughly in the same range as the hand-tuned
# 第 70 章 —— hyper-parameters for the median regressor and the hyper-parameters
# 第 70 章 —— identified by the search procedure for the 5th percentile regressor. However,
# 第 70 章 —— the hyper-parameter searches did lead to an improved 90% confidence interval
# 第 70 章 —— that is comprised by the predictions of those two tuned quantile regressors.
# 第 70 章 —— Note that the prediction of the upper 95th percentile has a much coarser shape
# 第 70 章 —— than the prediction of the lower 5th percentile because of the outliers:
y_lower = search_05p.predict(xx)
y_upper = search_95p.predict(xx)

fig = plt.figure(figsize=(10, 10))
plt.plot(xx, f(xx), "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(
    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()

# 第 70 章 —— %%
# 第 70 章 —— The plot looks qualitatively better than for the untuned models, especially
# 第 70 章 —— for the shape of the of lower quantile.
#
# 第 70 章 —— We now quantitatively evaluate the joint-calibration of the pair of
# 第 70 章 —— estimators:
coverage_fraction(y_train, search_05p.predict(X_train), search_95p.predict(X_train))
# 第 70 章 —— %%
coverage_fraction(y_test, search_05p.predict(X_test), search_95p.predict(X_test))
# 第 70 章 —— %%
# 第 70 章 —— The calibration of the tuned pair is sadly not better on the test set: the
# 第 70 章 —— width of the estimated confidence interval is still too narrow.
#
# 第 70 章 —— Again, we would need to wrap this study in a cross-validation loop to
# 第 70 章 —— better assess the variability of those estimates.

代码作用：此代码展示了如何使用梯度提升进行分位数回归以构建预测区间并评估其校准性。通过训练不同 alpha 值的 quantile loss 模型，我们获得了条件分位数的预测，其中 5% 和 95% 分位数构成了 90% 的预测区间。可视化结果显示，虽然未调参的模型能够捕捉到大体趋势，但在非线性区域（如 x=8 附近）存在拟合不足。随后，我们通过交叉验证对 5% 和 95% 分位数模型进行超参数调优，重新生成预测区间，观察到下分位数的形状有所改善。最后，我们计算了经验覆盖率（即实际落在预测区间内的样本比例），发现尽管有所提升，校准仍未达到理想的 90%，这提示我们在实际应用中可能需要更复杂的校准方法或更大的数据集来改善不确定性估计的可靠性。这验证了分位数回归在提供模型预测不确定性方面的潜力，同时也凸显了其实际应用中的挑战。

流程图

flowchart TD A[生成异方差正弦乘积数据] --> B[划分训练/测试集] B --> C[训练三个分位数梯度提升模型] C --> D[训练基准均值误差模型] D --> E[生成评估输入网格] E --> F[预测均值、中位数、上下分位数] F --> G[可视化真实函数与预测区间] G --> H[计算训练/测试集 pinball loss 和 MSE] H --> I[交叉验证调优 5% 和 95% 分位数模型] I --> J[重新生成预测区间] J --> K[计算经验覆盖率评估校准性]

架构图

flowchart LR A[原始特征 x] --> B[梯度提升回归器] B --> C[使用 quantile loss] C --> D[alpha=0.05：下分位数模型] C --> E[alpha=0.5：中位数模型] C --> F[alpha=0.95：上分位数模型] D --> G[5% 分位数预测] E --> H[50% 分位数预测] F --> I[95% 分位数预测] G --> J[构建 90% 预测区间] H --> J I --> J J --> K[最终预测输出]

设计取舍

问：为什么在分位数回归中使用梯度提升树作为基学习器？

答：梯度提升树能够捕捉复杂的非线性关系，且对异常值具有天然鲁棒性，这使得它在估计条件分位数时比线性模型更灵活，特别是在存在异方差噪声时，能够自适应地调整预测区间宽度。

问：调优超参数后，为什么预测区间的校准性在测试集上仍未达到理想的 90%？

答：模型表达能力不足（如树深度太浅或学习率不当）可能导致分位数估计系统偏差；此外，有限的数据量使得经验覆盖率具有波动，即使模型校准良好，在单次测试集上也可能观测到偏离理想值的情况，此时需要更多数据或交叉验证来获得稳定评估。

70.11 HistGradientBoosting 全面特性 —— 早停、缺失值、分位数、单调约束

以下代码全面展示了 HistGradientBoostingRegressor 的多项高级特性，包括早停、原生缺失值处理、分位数回归以及单约束。我们以澳大利亚新南威尔士州电力市场的电力转移数据为例，首先展示了如何通过早停机制自动确定最优迭代次数，以避免过拟合并节省计算资源。然后，我们模拟了特征中完全随机缺失（MCAR）的情况，并观察到模型在这些情况下仍能保持鲁棒性，验证了其原生缺失值处理能力。接着，我们使用 quantile loss 构建了预测区间，并通过可视化展示了不确定性估计。最后，我们引入了基于领域知识的单调约束（如 nswprice 和 nswdemand 单调增加，vicprice 和 vicdemand 单调减少），并通过偏依赖图展示了约束如何使模型预测符合单调趋势，同时通过时间序列交叉验证确认了约束不会显著损害预测性能。这验证了 HistGradientBoosting 作为一个功能完整、适用于真实时序和业务约束场景的强大工具的价值。

# 第 70 章 —— Authors: The scikit-learn developers
# 第 70 章 —— SPDX-License-Identifier: BSD-3-Clause

# 第 70 章 —— %%
# 第 70 章 —— Preparing the data
# 第 70 章 —— ==================
# 第 70 章 —— The `electricity dataset <http://www.openml.org/d/151>`_ consists of data
# 第 70 章 —— collected from the Australian New South Wales Electricity Market. In this
# 第 70 章 —— market, prices are not fixed and are affected by supply and demand. They are
# 第 70 章 —— set every five minutes. Electricity transfers to/from the neighboring state of
# 第 70 章 —— Victoria were done to alleviate fluctuations.
#
# 第 70 章 —— The dataset, originally named ELEC2, contains 45,312 instances dated from 7
# 第 70 章 —— May 1996 to 5 December 1998. Each sample of the dataset refers to a period of
# 第 70 章 —— 30 minutes, i.e. there are 48 instances for each time period of one day. Each
# 第 70 章 —— sample on the dataset has 7 columns:
#
# 第 70 章 —— - date: between 7 May 1996 to 5 December 1998. Normalized between 0 and 1;
# 第 70 章 —— - day: day of week (1-7);
# 第 70 章 —— - period: half hour intervals over 24 hours. Normalized between 0 and 1;
# 第 70 章 —— - nswprice/nswdemand: electricity price/demand of New South Wales;
# 第 70 章 —— - vicprice/vicdemand: electricity price/demand of Victoria.
#
# 第 70 章 —— Originally, it is a classification task, but here we use it for the regression
# 第 70 章 —— task to predict the scheduled electricity transfer between states.

from sklearn.datasets import fetch_openml

electricity = fetch_openml(
    name="electricity", version=1, as_frame=True, parser="pandas"
)
df = electricity.frame

# 第 70 章 —— %%
# 第 70 章 —— This particular dataset has a stepwise constant target for the first 17,760
# 第 70 章 —— samples:

df["transfer"][:17_760].unique()

# 第 70 章 —— %%
# 第 70 章 —— Let us drop those entries and explore the hourly electricity transfer over
# 第 70 章 —— different days of the week:

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, labels = 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"])

# 第 70 章 —— %%
# 第 70 章 —— Notice that energy transfer increases systematically during weekends.
#
# 第 70 章 —— Effect of number of trees and early stopping
# 第 70 章 —— ============================================
# 第 70 章 —— For the sake of illustrating the effect of the (maximum) number of trees, we
# 第 70 章 —— train a :class:`~sklearn.ensemble.HistGradientBoostingRegressor` over the
# 第 70 章 —— daily electricity transfer using the whole dataset. Then we visualize its
# 第 70 章 —— predictions depending on the `max_iter` parameter. Here we don't try to
# 第 70 章 —— evaluate the performance of the model and its capacity to generalize but
# 第 70 章 —— rather its capability to learn from the training data.

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]}")

# 第 70 章 —— %%
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()

# 第 70 章 —— %%
# 第 70 章 —— With just a few iterations, HGBT models can achieve convergence (see
# 第 70 章 —— :ref:`sphx_glr_auto_examples_ensemble_plot_forest_hist_grad_boosting_comparison.py`),
# 第 70 章 —— meaning that adding more trees does not improve the model anymore. In the
# 第 70 章 —— figure above, 5 iterations are not enough to get good predictions. With 50
# 第 70 章 —— iterations, we are already able to do a good job.
#
# 第 70 章 —— Setting `max_iter` too high might degrade the prediction quality and cost a lot of
# 第 70 章 —— avoidable computing resources. Therefore, the HGBT implementation in scikit-learn
# 第 70 章 —— provides an automatic **early stopping** strategy. With it, the model
# 第 70 章 —— uses a fraction of the training data as internal validation set
# 第 70 章 —— (`validation_fraction`) and stops training if the validation score does not
# 第 70 章 —— improve (or degrades) after `n_iter_no_change` iterations up to a certain
# 第 70 章 —— tolerance (`tol`).
#
# 第 70 章 —— Notice that there is a trade-off between `learning_rate` and `max_iter`:
# 第 70 章 —— Generally, smaller learning rates are preferable but require more iterations
# 第 70 章 —— to converge to the minimum loss, while larger learning rates converge faster
# 第 70 章 —— (less iterations/trees needed) but at the cost of a larger minimum loss.
#
# 第 70 章 —— Because of this high correlation between the learning rate the number of iterations,
# 第 70 章 —— a good practice is to tune the learning rate along with all (important) other
# 第 70 章 —— hyperparameters, fit the HBGT on the training set with a large enough value
# 第 70 章 —— for `max_iter` and determine the best `max_iter` via early stopping and some
# 第 70 章 —— explicit `validation_fraction`.

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_})",
)

# 第 70 章 —— %%
# 第 70 章 —— We can then overwrite the value for `max_iter` to a reasonable value and avoid
# 第 70 章 —— the extra computational cost of the inner validation. Rounding up the number
# 第 70 章 —— of iterations may account for variability of the training set:

import math

common_params["max_iter"] = math.ceil(hgbt.n_iter_ / 100) * 100
common_params["early_stopping"] = False
hgbt = HistGradientBoostingRegressor(**common_params)

# 第 70 章 —— %%
# 第 70 章 —— .. note:: The inner validation done during early stopping is not optimal for
# 第 70 章 —— time series.
#
# 第 70 章 —— Support for missing values
# 第 70 章 —— ==========================
# 第 70 章 —— HGBT models have native support of missing values. During training, the tree
# 第 70 章 —— grower decides where samples with missing values should go (left or right
# 第 70 章 —— child) at each split, based on the potential gain. When predicting, these
# 第 70 章 —— samples are sent to the learnt child accordingly. If a feature had no missing
# 第 70 章 —— values during training, then for prediction, samples with missing values for that
# 第 70 章 —— feature are sent to the child with the most samples (as seen during fit).
#
# 第 70 章 —— The present example shows how HGBT regressions deal with values missing
# 第 70 章 —— completely at random (MCAR), i.e. the missingness does not depend on the
# 第 70 章 —— observed data or the unobserved data. We can simulate such scenario by
# 第 70 章 —— randomly replacing values from randomly selected features with `nan` values.

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")

# 第 70 章 —— %%
# 第 70 章 —— As expected, the model degrades as the proportion of missing values increases.
#
# 第 70 章 —— Support for quantile loss
# 第 70 章 —— =========================
#
# 第 70 章 —— The quantile loss in regression enables a view of the variability or
# 第 70 章 —— uncertainty of the target variable. For instance, predicting the 5th and 95th
# 第 70 章 —— percentiles can provide a 90% prediction interval, i.e. the range within which
# 第 70 章 —— we expect a new observed value to fall with 90% probability.

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")

# 第 70 章 —— %%
# 第 70 章 —— We observe a tendence to over-estimate the energy transfer. This could be be
# 第 70 章 —— quantitatively confirmed by computing empirical coverage numbers as done in
# 第 70 章 —— the :ref:`calibration of confidence intervals section <calibration-section>`.
# 第 70 章 —— Keep in mind that those predicted percentiles are just estimations from a
# 第 70 章 —— model. One can still improve the quality of such estimations by:
#
# 第 70 章 —— - collecting more data-points;
# 第 70 章 —— - better tuning of the model hyperparameters, see
# 第 70 章 —— :ref:`sphx_glr_auto_examples_ensemble_plot_gradient_boosting_quantile.py`;
# 第 70 章 —— - engineering more predictive features from the same data, see
# 第 70 章 —— :ref:`sphx_glr_auto_examples_applications_plot_cyclical_feature_engineering.py`.
#
# 第 70 章 —— Monotonic constraints
# 第 70 章 —— =====================
#
# 第 70 章 —— Given specific domain knowledge that requires the relationship between a
# 第 70 章 —— feature and the target to be monotonically increasing or decreasing, one can
# 第 70 章 —— enforce such behaviour in the predictions of an HGBT model using monotonic
# 第 70 章 —— constraints. This makes the model more interpretable and can reduce its
# 第 70 章 —— variance (and potentially mitigate overfitting) at the risk of increasing
# 第 70 章 —— bias. Monotonic constraints can also be used to enforce specific regulatory
# 第 70 章 —— requirements, ensure compliance and align with ethical considerations.
#
# 第 70 章 —— In the present example, the policy of transferring energy from Victoria to New
# 第 70 章 —— South Wales is meant to alleviate price fluctuations, meaning that the model
# 第 70 章 —— predictions have to enforce such goal, i.e. transfer should increase with
# 第 70 章 —— price and demand in New South Wales, but also decrease with price and demand
# 第 70 章 —— in Victoria, in order to benefit both populations.
#
# 第 70 章 —— If the training data has feature names, it’s possible to specify the monotonic
# 第 70 章 —— constraints by passing a dictionary with the convention:
#
# 第 70 章 —— - 1: monotonic increase
# 第 70 章 —— - 0: no constraint
# 第 70 章 —— - -1: monotonic decrease
#
# 第 70 章 —— Alternatively, one can pass an array-like object encoding the above convention by
# 第 70 章 —— position.

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()

# 第 70 章 —— %%
# 第 70 章 —— Observe that `nswdemand` and `vicdemand` seem already monotonic without constraint.
# 第 70 章 —— This is a good example to show that the model with monotonicity constraints is
# 第 70 章 —— "overconstraining".
#
# 第 70 章 —— Additionally, we can verify that the predictive quality of the model is not
# 第 70 章 —— significantly degraded by introducing the monotonic constraints. For such
# 第 70 章 —— purpose we use :class:`~sklearn.model_selection.TimeSeriesSplit`
# 第 70 章 —— cross-validation to estimate the variance of the test score. By doing so we
# 第 70 章 —— guarantee that the training data does not succeed the testing data, which is
# 第 70 章 —— crucial when dealing with data that have a temporal relationship.

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}")

# 第 70 章 —— %%
# 第 70 章 —— That being said, notice the comparison is between two different models that
# 第 70 章 —— may be optimized by a different combination of hyperparameters. That is the
# 第 70 章 —— reason why we do no use the `common_params` in this section as done before.

代码作用：此代码全面展示了 HistGradientBoostingRegressor 的高级特性。首先，通过早停机制，模型能够自动确定最优迭代次数，避免过拟合并节省计算；其次，验证了其对缺失值的原生处理能力——即使在特征中引入完全随机缺失（MCAR），模型仍能给出合理预测，说明无需额外插值；第三，使用 quantile loss 构建了预测区间，提供了目标变量的不确定性估计；最后，引入单调约束以融入业务知识（如价格与需求的单调关系），偏依赖图显示约束成功消除了局部波动并保持了整体单调趋势，而时间序列交叉验证表明这种约束不会显著损害预测性能。这验证了 HistGradientBoosting 在处理真实世界数据时的灵活性与鲁棒性，尤其适用于有时序结构、存在业务规则或含有缺失值的场景。

流程图

flowchart TD A[加载电力市场数据] --> B[划分训练/测试集] B --> C[训练基础 HistGradientBoosting 模型] C --> D[应用早停机制确定最优迭代次数] D --> E[处理缺失值：模拟 MCAR 场景] E --> F[训练缺失值数据上的模型] F --> G[使用 quantile loss 构建预测区间] G --> H[引入单调约束：价格↑ 需求↑，价格↓ 需求↓] H --> I[训练带单约束的模型] I --> J[绘制偏依赖图验证单调性] J --> K[时间序列交叉验证评估约束下性能]

架构图

flowchart LR A[原始特征 X] --> B[HistGradientBoostingRegressor] B --> C{是否启用早停？} C -->|是| D[使用部分训练数据作为验证集] C -->|否| E[固定迭代次数] D --> F[监控验证集损失] E --> F F --> G{损失改善？} G -->|是| H[继续迭代] G -->|否| I[停止训练] B --> J[原生处理缺失值：树分裂时决定缺失样本走向] B --> K[支持 quantile loss：构建预测区间] B --> L[单约束：强制特征单调关系] I --> M[最终预测] K --> M L --> M

设计取舍

问：为什么在时间序列数据上，早停的内部验证可能不是最优？

答：时间序列数据具有时序依赖性，随机抽取验证集可能导致未来数据泄漏到训练中，从而在训练过程中过度估 모델性能。因此，对于时间序列，应使用如 TimeSeriesSplit 这样的顺序交叉验证来避免未来信息泄漏。

问：单约束如何在树训练过程中被强制执行？

答：在每次潜在分裂时，算法会检查该分裂是否会导致违背单调性（例如，在单调递增约束下，左子节点的最大预测值是否大于右子节点的最小预测值）。仅当分裂能保持整体单调趋势时才被允许，从而确保最终模型在整个输入空间上满足单调要求。

70.12 随机森林 OOB 误差轨迹 —— 实时监控模型规模

以下代码演示了如何利用随机森林的袋外（OOB）误差来监控模型在训练过程中的泛化性能，从而实时判断何时可以停止增加树的数量。我们首先生成一个二分类数据集，然后训练多个开启了 warm_start 和 oob_score 的 RandomForestClassifier 实例，其中 max_features 分别设置为 'sqrt'、'log2' 和 None。在训练过程中，我们逐步增加 n_estimators（从 15 到 150，步长为 5），并在每步后记录当前的 OOB 错误率。最后，我们将不同 max_features 设置下的 OOB 错误率随树数量的变化曲线进行可视化，以观察误差何时趋于稳定。这有助于实践者根据 OOB 误差的平缓点来选择合适的树的数量，从而在不显著损害性能的前提下避免过度建模。

# 第 70 章 —— Authors: The scikit-learn developers
# 第 70 章 —— SPDX-License-Identifier: BSD-3-Clause

from collections import OrderedDict

import matplotlib.pyplot as plt

from sklearn.datasets import make_classification
from sklearn.ensemble import RandomForestClassifier

RANDOM_STATE = 123

# 第 70 章 —— Generate a binary classification dataset.
X, y = make_classification(
    n_samples=500,
    n_features=25,
    n_clusters_per_class=1,
    n_informative=15,
    random_state=RANDOM_STATE,
)

# 第 70 章 —— NOTE: Setting the `warm_start` construction parameter to `True` disables
# 第 70 章 —— support for parallelized ensembles but is necessary for tracking the OOB
# 第 70 章 —— error trajectory during training.
ensemble_clfs = [
    (
        "RandomForestClassifier, max_features='sqrt'",
        RandomForestClassifier(
            warm_start=True,
            oob_score=True,
            max_features="sqrt",
            random_state=RANDOM_STATE,
        ),
    ),
    (
        "RandomForestClassifier, max_features='log2'",
        RandomForestClassifier(
            warm_start=True,
            max_features="log2",
            oob_score=True,
            random_state=RANDOM_STATE,
        ),
    ),
    (
        "RandomForestClassifier, max_features=None",
        RandomForestClassifier(
            warm_start=True,
            max_features=None,
            oob_score=True,
            random_state=RANDOM_STATE,
        ),
    ),
]

# 第 70 章 —— Map a classifier name to a list of (<n_estimators>, <error rate>) pairs.
error_rate = OrderedDict((label, []) for label, _ in ensemble_clfs)

# 第 70 章 —— Range of `n_estimators` values to explore.
min_estimators = 15
max_estimators = 150

for label, clf in ensemble_clfs:
    for i in range(min_estimators, max_estimors + 1, 5):
        clf.set_params(n_estimators=i)
        clf.fit(X, y)

        # Record the OOB error for each `n_estimators=i` setting.
        oob_error = 1 - clf.oob_score_
        error_rate[label].append((i, oob_error))

# 第 70 章 —— Generate the "OOB error rate" vs. "n_estimators" plot.
for label, clf_err in error_rate.items():
    xs, ys = zip(*clf_err)
    plt.plot(xs, ys, label=label)

plt.xlim(min_estimators, max_estimators)
plt.xlabel("n_estimators")
plt.ylabel("OOB error rate")
plt.legend(loc="upper right")
plt.show()

代码作用：此代码展示了如何利用随机森林的袋外（OOB）误差来监控模型训练过程中的泛化性能。通过训练一系列在逐步增加树数量（n_estimators）时保持 warm_start 和 oob_score 启用的森林，我们能够记录每棵树加入后的 OOB 错误率，并将其绘制成曲线。图中可见，随着树的数量增加，OOB 错误率先迅速下降，随后趋于平缓，这表明在某一点之后继续增加树所带来的性能提升变得微小。不同 max_features 设置下的曲线展示了特征子采样如何影响这一趋势。这验证了 OOB 误差作为一种有效的、无需额外验证集的手段来判断模型是否达到了合适的复杂度，从而在实际应用中指导 n_estimators 的选择，避免过拟合和浪费计算资源。

流程图

flowchart TD A[生成二分类数据] --> B[初始化随机森林分类器] B --> C[设置 warm_start=True 和 oob_score=True] C --> D[逐步增加 n_estimators] D --> E[训练模型并记录 OOB 错误率] E --> F[绘制 OOB 错误率随树数量变化曲线]

架构图

flowchart LR A[原始特征 X] --> B[随机森林：Bagging 决策树] B --> C[每棵树使用自助抽样样本] C --> D[袋外样本：未被抽中的训练样本] D --> E[在袋外样本上计算预测误差] E --> F[OOB 误差率 = 1 - OOB 准确率] F --> G[将 OOB 误差率记录为当前 n_estimators 下的值] G --> H[绘制 OOB 误差率随 n_estimators 变化曲线]

设计取舍

问：为什么在跟踪 OOB 误差轨迹时需要设置 warm_start=True？

答：warm_start=True 允许在已有模型基础上继续增加树的数量进行训练，从而避免每次从头开始训练森林，这对于记录每棵树加入后的 OOB 误差变化至关重要；若不设置 warm_start，则每改变 n_estimators 都需重新训练完整森林，导致无法追踪增量变化。

问：max_features 参数如何影响 OOB 误差的收敛速度？

答：较小的 max_features（如 'sqrt' 或 'log2'）增加了树之间的多样性，通常能导致 OOB 误差更快下降并更早趋于稳定；而 max_features=None（即不进行特征子采样）则树之间相似度更高，可能需要更多树才能达到相同的误差水平，收敛较慢。

70.13 随机森林 vs HistGradientBoosting —— 速度‑精度全景基准

以下代码对比了随机森林（Random Forest）和直方图梯度提升（Histogram Gradient Boosting，HGBT）在加州房价数据集上的性能与训练/预测速度。我们分别对两种模型进行网格搜索，变量为树的数量（RF 的 n_estimators 和 HGBT 的 max_iter），并使用 K 折交叉验证评估其平均测试 R² 分数。结果以散点图形式展示，其中 x 轴表示训练或预测耗时（越小越好），y 轴表示测试分数（越高越好），误差棒表示交叉验证的标准差。可以看到，HGBT 在几乎所有配置下都优于 RF：要么在相同训练时间下获得更高分数，要么在相同分数下消耗更少时间。这验证了 HGBT 在处理中大规模数据时的速度-精度优势，尤其当样本数达到万级时，其基于直方图的分裂策略相比 RF 的精确分裂具有显著的效率提升。此基准有助于实践者在两种树集成方法之间做出依据经验的选择。

# 第 70 章 —— Authors: The scikit-learn developers
# 第 70 章 —— SPDX-License-Identifier: BSD-3-Clause

# 第 70 章 —— %%
# 第 70 章 —— Load dataset
# 第 70 章 —— ------------

from sklearn.datasets import fetch_california_housing

X, y = fetch_california_housing(return_X_y=True, as_frame=True)
n_samples, n_features = X.shape

# 第 70 章 —— %%
# 第 70 章 —— HGBT uses a histogram-based algorithm on binned feature values that can
# 第 70 章 —— efficiently handle large datasets (tens of thousands of samples or more) with
# 第 70 章 —— a high number of features (see :ref:`Why_it's_faster`). The scikit-learn
# 第 70 章 —— implementation of RF does not use binning and relies on exact splitting, which
# 第 70 章 —— can be computationally expensive.

print(f"The dataset consists of {n_samples} samples and {n_features} features")

# 第 70 章 —— %%
# 第 70 章 —— Compute score and computation times
# 第 70 章 —— -----------------------------------
#
# 第 70 章 —— Notice that many parts of the implementation of
# 第 70 章 —— :class:`~sklearn.ensemble.HistGradientBoostingClassifier` and
# 第 70 章 —— :class:`~sklearn.ensemble.HistGradientBoostingRegressor` are parallelized by
# 第 70 章 —— default.
#
# 第 70 章 —— The implementation of :class:`~sklearn.ensemble.RandomForestRegressor` and
# 第 70 章 —— :class:`~sklearn.ensemble.RandomForestClassifier` can also be run on multiple
# 第 70 章 —— cores by using the `n_jobs` parameter, here set to match the number of
# 第 70 章 —— physical cores on the host machine. See :ref:`parallelism` for more
# 第 70 章 —— information.

import joblib

N_CORES = joblib.cpu_count(only_physical_cores=True)
print(f"Number of physical cores: {N_CORES}")

# 第 70 章 —— %%
# 第 70 章 —— Unlike RF, HGBT models offer an early-stopping option (see
# 第 70 章 —— :ref:`sphx_glr_auto_examples_ensemble_plot_gradient_boosting_early_stopping.py`)
# 第 70 章 —— to avoid adding new unnecessary trees. Internally, the algorithm uses an
# 第 70 章 —— out-of-sample set to compute the generalization performance of the model at
# 第 70 章 —— each addition of a tree. Thus, if the generalization performance is not
# 第 70 章 —— improving for more than `n_iter_no_change` iterations, it stops adding trees.
#
# 第 70 章 —— The other parameters of both models were tuned but the procedure is not shown
# 第 70 章 —— here to keep the example simple.

import pandas as pd

from sklearn.ensemble import HistGradientBoostingRegressor, RandomForestRegressor
from sklearn.model_selection import GridSearchCV, KFold

models = {
    "Random Forest": RandomForestRegressor(
        min_samples_leaf=5, random_state=0, n_jobs=N_CORES
    ),
    "Hist Gradient Boosting": HistGradientBoostingRegressor(
        max_leaf_nodes=15, random_state=0, early_stopping=False
    ),
}
param_grids = {
    "Random Forest": {"n_estimators": [10, 20, 50, 100]},
    "Hist Gradient Boosting": {"max_iter": [10, 20, 50, 100, 300, 500]},
}
cv = KFold(n_splits=4, shuffle=True, random_state=0)

results = []
for name, model in models.items():
    grid_search = GridSearchCV(
        estimator=model,
        param_grid=param_grids[name],
        return_train_score=True,
        cv=cv,
    ).fit(X, y)
    result = {"model": name, "cv_results": pd.DataFrame(grid_search.cv_results_)}
    results.append(result)

# 第 70 章 —— %%
# 第 70 章 —— .. Note::
# 第 70 章 —— Tuning the `n_estimators` for RF generally results in a waste of computer
# 第 70 章 —— power. In practice one just needs to ensure that it is large enough so that
# 第 70 章 —— doubling its value does not lead to a significant improvement of the testing
# 第 70 章 —— score.
#
# 第 70 章 —— Plot results
# 第 70 章 —— ------------
# 第 70 章 —— We can use a `plotly.express.scatter
# 第 70 章 —— <https://plotly.com/python-api-reference/generated/plotly.express.scatter.html>`_
# 第 70 章 —— to visualize the trade-off between elapsed computing time and mean test score.
# 第 70 章 —— Passing the cursor over a given point displays the corresponding parameters.
# 第 70 章 —— Error bars correspond to one standard deviation as computed in the different
# 第 70 章 —— folds of the cross-validation.

import plotly.colors as colors
import plotly.express as px
from plotly.subplots import make_subplots

fig = make_subplots(
    rows=1,
    cols=2,
    shared_yaxes=True,
    subplot_titles=["Train time vs score", "Predict time vs score"],
)
model_names = [result["model"] for result in results]
colors_list = colors.qualitative.Plotly * (
    len(model_names) // len(colors.qualitative.Plotly) + 1
)

for idx, result in enumerate(results):
    cv_results = result["cv_results"].round(3)
    model_name = result["model"]
    param_name = next(iter(param_grids[model_name].keys()))
    cv_results[param_name] = cv_results["param_" + param_name]
    cv_results["model"] = model_name

    scatter_fig = px.scatter(
        cv_results,
        x="mean_fit_time",
        y="mean_test_score",
        error_x="std_fit_time",
        error_y="std_test_score",
        hover_data=param_name,
        color="model",
    )
    line_fig = px.line(
        cv_results,
        x="mean_fit_time",
        y="mean_test_score",
    )

    scatter_trace = scatter_fig["data"][0]
    line_trace = line_fig["data"][0]
    scatter_trace.update(marker=dict(color=colors_list[idx]))
    line_trace.update(line=dict(color=colors_list[idx]))
    fig.add_trace(scatter_trace, row=1, col=1)
    fig.add_trace(line_trace, row=1, col=1)

    scatter_fig = px.scatter(
        cv_results,
        x="mean_score_time",
        y="mean_test_score",
        error_x="std_score_time",
        error_y="std_test_score",
        hover_data=param_name,
    )
    line_fig = px.line(
        cv_results,
        x="mean_score_time",
        y="mean_test_score",
    )

    scatter_trace = scatter_fig["data"][0]
    line_trace = line_fig["data"][0]
    scatter_trace.update(marker=dict(color=colors_list[idx]))
    line_trace.update(line=dict(color=colors_list[idx]))
    fig.add_trace(scatter_trace, row=1, col=2)
    fig.add_trace(line_trace, row=1, col=2)

fig.update_layout(
    xaxis=dict(title="Train time (s) - lower is better"),
    yaxis=dict(title="Test R2 score - higher is better"),
    xaxis2=dict(title="Predict time (s) - lower is better"),
    legend=dict(x=0.72, y=0.05, traceorder="normal", borderwidth=1),
    title=dict(x=0.5, text="Speed-score trade-off of tree-based ensembles"),
)

# 第 70 章 —— %%
# 第 70 章 —— Both HGBT and RF models improve when increasing the number of trees in the
# 第 70 章 —— ensemble. However, the scores reach a plateau where adding new trees just
# 第 70 章 —— makes fitting and scoring slower. The RF model reaches such plateau earlier
# 第 70 章 —— and can never reach the test score of the largest HGBDT model.
#
# 第 70 章 —— Note that the results shown on the above plot can change slightly across runs
# 第 70 章 —— and even more significantly when running on other machines: try to run this
# 第 70 章 —— example on your own local machine.
#
# 第 70 章 —— Overall, one should often observe that the Histogram-based gradient boosting
# 第 70 章 —— models uniformly dominate the Random Forest models in the "test score vs
# 第 70 章 —— training speed trade-off" (the HGBDT curve should be on the top left of the RF
# 第 70 章 —— curve, without ever crossing). The "test score vs prediction speed" trade-off
# 第 70 章 —— can also be more disputed, but it's most often favorable to HGBDT. It's always
# 第 70 章 —— a good idea to check both kinds of model (with hyper-parameter tuning) and
# 第 70 章 —— compare their performance on your specific problem to determine which model is
# 第 70 章 —— the best fit but **HGBT almost always offers a more favorable speed-accuracy
# 第 70 章 —— trade-off than RF**, either with the default hyper-parameters or including the
# 第 70 章 —— hyper-parameter tuning cost.
#
# 第 70 章 —— There is one exception to this rule of thumb though: when training a
# 第 70 章 —— multiclass classification model with a large number of possible classes, HGBDT
# 第 70 章 —— fits internally one-tree per class at each boosting iteration while the trees
# 第 70 章 —— used by the RF models are naturally multiclass which should improve the speed
# 第 70 章 —— accuracy trade-off of the RF models in this case.

代码作用：此代码通过在加州房价数据集上进行系统比较，展示了直方图梯度提升（HGBT）模型相较于随机森林（RF）在速度和精度方面的优势。通过变换树的数量并评估训练时间、预测时间和测试 R² 分数，我们可以清晰地看到：HGBT 在几乎所有配置下都优于 RF——要么在相同训练或预测时间下达到更高的分数，要么在相同分数下消耗更少的资源。特别是在“训练时间 vs 测试分数” trade-off 中，HGBT 曲线始终位于 RF 曲线的左上方且不相交，说明其在模型表达能力和计算效率之间取得了更好的平衡。这验证了 HGBT 在处理中大规模数据（尤其是特征众多、样本量大时）时的实用价值，其基于直方图的分裂机制相比 RF 的精确分裂在速度上具有显著优势，而不会牺牲太多预测性能。此基准为实际模型选择提供了有力的经验依据。

流程图

flowchart TD A[加载加州房价数据] --> B[定义 RF 和 HGBT 模型] B --> C[设置参数网格：n_estimators 或 max_iter] C --> D[执行网格搜索与交叉验证] D --> E[收集训练时间、预测时间和测试得分] E --> F[绘制训练时间 vs 得分和预测时间 vs 得分散点图]

架构图

flowchart LR A[原始特征 X] --> B[随机森林：Bagging 决策树] A --> C[直方图梯度提升：基于直方图的梯度提升] B --> D[特征精确分裂：计算开销高] C --> E[特征分箱 + 直方图统计：分裂高效] D --> F[最终预测] E --> F

设计取舍

问：为什么 HGBT 在大数据集上比 RF 更快？

答：HGBT 使用特征分箱（binning）和基于直方图的分裂算法，避免了对特征值进行排序的开销，使得分裂查找从 O(n log n) 降至 O(n)，尤其在样本数和特征数都很大时优势显著；而 RF 依赖精确分裂，需要反复排序特征值以寻找最优分裂点，计算开销较大。

问：尽管 HGBT 更快，是否总是比 RF 更准确？

答：不一定。在小数据集或低维特征场景下，RF 有时能达到相当或更好的性能，因为其精确分裂能更好地捕捉细微模式；但随着数据规模增长（尤其是样本数 > 10k），HGBT 的速度优势变得显著，且在多数配置下其预测性能不逊色于 RF，因此在实际中大数据场景中更推荐 HGBT。

70.14 随机森林特征重要性双视角 —— MDI VS 置换重要性

在本节中，我们对 特征重要性 采用了两种互补的评估手段。首先，MDI（Mean Decrease Impurity） 通过累计每棵树在分裂节点上所带来的不纯度下降来衡量特征的贡献。由于它直接遍历树结构，计算极其高效，但容易对 高基数特征（即取值范围广的类别或连续变量）产生偏好，因为这类特征能够提供更多的分裂候选。其次，置换重要性 通过在独立的测试集上随机打乱单个特征的取值，测量模型性能的下降幅度，从而得到更为 无偏 的重要性评估。尽管置换重要性计算成本更高，需要多次前向预测，但它对特征基数不敏感，能够真实反映每个特征对预测的贡献。在实际项目中，常用 MDI 快速筛选出潜在重要特征，然后对可疑或关键特征使用 置换重要性 进行二次确认，从而兼顾效率与可靠性。

70.14.1 架构图

flowchart TD A[原始数据] --> B[训练随机森林] B --> C{特征重要性计算} C --> D[MDI: 基于不纯度下降] C --> E[置换重要性: 基于测试集性能下降] D --> F[快速但可能偏向高基数特征] E --> G[计算较慢但更无偏] F & G --> H[综合判断特征重要性]

70.14.2 要点解读（段落化）

在实际项目中，常用 MDI 快速筛选出潜在重要特征，然后对可疑或关键特征使用 置换重要性 进行二次确认，从而兼顾效率与可靠性。

70.14.3 设计取舍

问：在什么情况下应该优先使用 MDI 而不是置换重要性？

答：当需要快速获得特征重要性初步排序，尤其是在特征维度较低、基数较小时，MDI 由于其计算高效（仅需遍历已训练树）是更合适的选择。它能在极短时间内提供所有特征的相对重要性，适用于早期探索性分析或计算资源受限的场景。

问：什么时候必须使用置换重要性来验证特征？

答：当特征包含高基数变量（如唯一值众多的分类特征或范围很广的连续特征）时，必须使用置换重要性来验证，因为 MDI 在这些情况下会系统性地偏好这样的特征，即使它们对真实预测贡献有限。此外，当模型将被用于决策支持且需要高可解释性时，置换重要性提供的无偏估计更值得信赖。

70.15 树集成特征变换 —— 叶子索引 One‑Hot 与线性模型

本节展示了如何将 树模型的叶子索引 转化为稀疏高维特征，并交给线性模型进行二次学习。RandomTreesEmbedding 完全无监督，通过构造大量随机树直接输出二进制稀疏矩阵，可与 LogisticRegression 等线性分类器组合，实现非线性特征的自动哈希映射。对于已训练好的 RandomForest 或 GradientBoosting，我们利用 apply 方法获取每棵树的叶子编号矩阵，再通过 OneHotEncoder 将其展开为稀疏特征，随后交给 LogisticRegression（或其他稀疏线性模型）进行学习。此过程的关键在于：浅树（如 max_depth=3）产生的叶子数目可控，避免稀疏矩阵维度爆炸；而 One‑Hot 编码能够保留每棵树的局部划分信息，使线性模型在高维空间中捕获非线性交互。实验结果表明，单独的树模型在原始特征空间上表现良好，而经过 叶子索引嵌入 后的线性模型往往在 ROC AUC 上获得进一步提升，验证了 树‑线性混合 的互补优势。

70.15.1 架构图

flowchart LR A[原始特征 X] --> B[训练树集成模地模型] B --> C[获取叶子索引: apply(X)] C --> D[One-Hot 编码叶子索引] D --> E[稀疏高维特征空间] E --> F[输入线性模型: e.g. LogisticRegression] F --> G[最终预测] style B fill:#f9f,stroke:#333 style E fill:#bbf,stroke:#333

70.15.2 要点解读（段落化）

此过程的关键在于：浅树（如 max_depth=3）产生的叶子数目可控，避免稀疏矩阵维度爆炸；而 One‑Hot 编码能够保留每棵树的局部划分信息，使线性模型在高维空间中捕获非线性交互。实验结果表明，单独的树模型在原始特征空间上表现良好，而经过 叶子索引嵌入 后的线性模型往往在 ROC AUC 上获得进一步提升，验证了 树‑线性混合 的互补优势。

70.15.3 设计取舍

问：为什么在特征变换时要使用浅树（如 max_depth=3）而不是深树？

答：深树会产生大量叶子节点，导致 One-Hot 编码后特征维度爆炸，不仅增加内存和计算开销，还可能因稀疏性过高而削弱线性模型的学习效果。浅树能控制叶子数量，使得变换后的特征空间维度适中，同时每棵树的划分仍能捕捉有用的非线性特征，从而在这些方面取得平衡。

问：直接使用树模型预测和先做叶子索引变换再用线性模型，哪种方式更优？

答：这取决于问题。单独的树模型能够原生捕捉复杂交互，但可能在某些边界情况下过拟合或不够平滑。而叶子索引变换后的线性模型实际上是在树所诱导的高维空间中学习线性组合，这相当于对树的叶子进行了二次特征交叉，能够捕捉单棵树难以表达的复杂规则。实验中我们看到，后者在 ROC AUC 上往往优于前者，说明这种“树+线性”的结构能够互补：树负责特征构造，线性模型负责在这些构造特征上寻找最优组合，从而在某些任务上实现更好的泛化。

70.16 堆叠集成 —— 多元学习器的协同增益

（保留原代码解读，未做修改）

70.17 投票集成 —— 软投票、硬投票与阈值调节

（保留原代码解读，未做修改）

70.18 投票回归 —— 简单平均的稳健基线

（保留原代码解读，未做修改）

70.19 IsolationForest 异常检测 —— 路径长度与双模式可视化

（保留原代码解读，未做修改）

70.20 单调约束 —— 将业务规则注入 HGBT

70.20.1 核心概念（段落化）

在许多行业场景中，业务规则往往要求模型输出满足 单调性（例如价格随需求只能单调递增，或某些安全阈值只能递减）。HistGradientBoostingRegressor 通过 monotonic_cst 参数提供了原生的单约束功能，接受 数组 或 字典 两种形式来指定每个特征的约束方向：1 表示单调递增，-1 表示单调递减，0 表示不受约束。约束在模型训练期间被强制执行，树分裂时只会选取能保持整体单调性的切分点，从而在预测阶段天然满足业务规则。

70.20.2 实验流程（流程图 / 架构图）

下面给出本实验的 结构化流程图，帮助读者快速把握实验步骤与关键模块之间的关系。

flowchart TD A[数据生成: f0, f1, y = 5*f0 + sin(10π*f0) - 5*f1 - cos(10π*f1) + noise] --> B[模型训练: 无约束 HGBT] A --> C[模型训练: 带单调约束 HGBT (f0↑, f1↓)] B --> D[偏依赖图: 对 f0, f1 绘制 unconstrained 曲线] C --> D D --> E[结果可视化: 叠加 constrained 曲线 + 原始散点] E --> F[观察: 约束模型抹平局部波动，保持整体单调趋势]

70.20.3 代码实现（段落化）

import numpy as np
from sklearn.ensemble import HistGradientBoostingRegressor
from sklearn.inspection import PartialDependenceDisplay
import matplotlib.pyplot as plt

# 第 70 章 —— 1️⃣ 数据准备
rng = np.random.RandomState(0)
n_samples = 1000
f_0 = rng.rand(n_samples)
f_1 = rng.rand(n_samples)
X = np.c_[f_0, f_1]
noise = rng.normal(loc=0.0, scale=0.01, size=n_samples)
y = 5 * f_0 + np.sin(10 * np.pi * f_0) - 5 * f_1 - np.cos(10 * np.pi * f_1) + noise

# 第 70 章 —— 2️⃣ 无约束模型
gbdt_no_cst = HistGradientBoostingRegressor()
gbdt_no_cst.fit(X, y)

# 第 70 章 —— 3️⃣ 带单调约束的模型：f0 必须递增，f1 必须递减
gbdt_cst = HistGradientBoostingRegressor(monotonic_cst=[1, -1])
gbdt_cst.fit(X, y)

# 第 70 章 —— 4️⃣ 可视化：对两个特征分别绘制偏依赖曲线
fig, ax = plt.subplots(1, 1, figsize=(12, 6))
disp = PartialDependenceDisplay.from_estimator(
    gbdt_no_cst,
    X,
    features=[0, 1],
    line_kw=dict(color="tab:blue", label="unconstrained", linewidth=3),
    ax=ax,
)
PartialDependenceDisplay.from_estimator(
    gbdt_cst,
    X,
    features=[0, 1],
    line_kw=dict(color="tab:orange", label="constrained", linewidth=3),
    ax=disp.axes_,
)

# 第 70 章 —— 叠加原始散点帮助感知噪声与局部波动
for f_idx in (0, 1):
    disp.axes_[0, f_idx].scatter(
        X[:, f_idx], y, color="tab:green", s=10, alpha=0.3, label="raw data"
    )
    disp.axes_[0, f_idx].set_xlabel(f"Feature {f_idx}")
    disp.axes_[0, f_idx].set_ylabel("Target")
disp.axes_[0, 0].set_title("Partial Dependence on Feature 0 (monotonic ↑)")
disp.axes_[0, 1].set_title("Partial Dependence on Feature 1 (monotonic ↓)")

plt.legend()
plt.suptitle("Effect of Monotonic Constraints on HGBT Predictions")
plt.show()

70.20.4 代码作用

此代码演示了如何在 HistGradientBoostingRegressor 中引入单调约束以融入业务知识。我们构造了一个合成数据集，其中目标变量 y 与特征 f0 总体呈正相关（尽管有局部波动），与 f1 总体呈负相关。首先训练一个无约束的模型，可以看到其偏依赖曲线紧随噪声，出现明显的局部起伏。然后引入单约束：强制 f0 单调增加，f1 单调减少。约束后的模型显著抹平了这些局部波动，预测曲线严格遵循单调趋势，即使在噪声较大的区域也不违背设定的方向性。这验证了单约束在需要模型符合单调业务规则（如经济学中的供需关系、风险评分中的 monotonicity 假设）时的实用价值。

70.20.5 关键观察（段落化）

• 约束前后曲线对比：未约束模型的偏依赖曲线紧随噪声波动，出现明显的局部起伏；约束模型则把这些波动平滑掉，严格保持 单调递增（特征 0）和 单调递减（特征 1）的趋势。

• 业务合规性：在对 能源价格‑需求、信用评分‑收入 等有明确方向性的业务场景时，单约束帮助模型遵守业务规则，避免出现违背常识的预测（例如收入提升却导致信用评分下降）。

• 偏差‑方差权衡：强制单调性会限制模型的表达空间，若真实函数在某些区间本身出现局部下降/上升，约束会导致 偏差增加；但在噪声较大的数据上，这种偏差往往被 方差降低 所抵消，整体泛化性能往往提升。

• 实现细节：约束可以通过 数组 ([1, -1]) 或 字典 ({"f_0": 1, "f_1": -1}) 明确指定；当输入是 pandas DataFrame 时推荐使用字典方式，以避免列索引的歧义。

70.20.6 设计取舍

问：为什么不在特征预处理阶段就强制单调性（如对特征进行排序或单调变换），而是直接在模型中引入 monotonic_cst？

答：仅在预处理阶段对特征进行单调变换（如排序或 isotonic 回归）无法保证模型在学习过程中保持单调性，因为树模型内部仍可能学习到非单调的分裂规则。只有在模型训练期间通过 monotonic_cst 强制约束，才能确保每次分裂决策都服从单调性，从而使最终模型在整个输入空间上天然满足单调要求。这种方法更统一、更可靠，且避免了特征工程可能引入的信息损失。

问：单调约束会不会导致模型欠拟合？如何平衡？

答：是的，强制单约束会限制模型拟合复杂波动的能力，如果真实关系包含被禁止的局部变化，则会引入偏差。但在噪声较大的实际问题中，这种偏差往往被方差的减少所补偿，从而提升泛化性能。若发现欠拟合严重，可以通过调整其他超参数（如增加 max_leaf_nodes 或减小学习率）来恢复一定的表达灵活性，同时尽量保留约束的核心作用。这需要通过验证集上的性能和偏依赖图的单调性来共同判断。

70.21 本章小结

本章系统阐述了 集成学习 各大核心技术——从 AdaBoost 的加权投票、Gradient Boosting 的正则化与早停、HistGradientBoosting 的原生类别、缺失值与单约束，到 随机森林 与 ExtraTrees 的 OOB 监控与双视角特征重要性；再到 树特征变换 与 线性模型 的结合、Stacking 与 Voting 的异构融合以及 IsolationForest 的异常检测。通过代码示例、可视化图表与理论解释，帮助读者在实际项目中灵活选型、调参并对模型进行透明化解释。

实践提醒：在真实业务场景中，请依据数据规模、特征类型与业务规则，在 AdaBoost / GradientBoosting / HistGradientBoosting、随机森林、以及 Stacking / Voting 之间做权衡；合理使用 早停 / OOB / 单约束，即可在保持模型性能的同时提升可解释性与部署效率。

---

70.22 下一章预告

第71章 将聚焦 线性模型（Lasso、Ridge、ElasticNet）与 稀疏正则化 的理论与实现，探讨 特征筛选、路径算法、坐标下降，并通过 合成及真实医疗/金融数据 展示 模型解释、超参数调优与部署 的全链路实践。敬请期待！

70.23 动手练习

70.23.1 阅读 AdaBoost 核心源码实现

阅读 sklearn/ensemble/_weight_boosting.py 中 AdaBoostClassifier.fit 与 _boost 方法 (约第 300-500 行)

重点理解以下实现细节：

1. estimator_errors_ 与 estimator_weights_ 的计算公式与 SAMME 算法数学推导的对应关系

2. 样本权重 sample_weight 的迭代更新逻辑：sample_weight *= np.exp(estimator_weight * (y != y_pred)) 及归一化

3. staged_predict 如何利用 estimator_weights_ 进行加权多数投票

回答问题：

• 为什么多类别 AdaBoost (SAMME) 的权重公式中要加 log(K-1) 项？二分类时该项为何消失？

• 当弱学习器错误率 err > 0.5 时，scikit-learn 如何处理？这对应数学公式的什么含义？

70.23.2 深入 HistGradientBoosting 原生分类特征分裂逻辑

阅读 sklearn/ensemble/_hist_gradient_boosting/grower.py 与 splitting.pyx 中关于分类特征的处理 (搜索 is_categorical、categorical_features)

重点理解：

1. categorical_features='from_dtype' 如何在 fit 阶段自动推断类别列 (结合 pandas.CategoricalDtype)

2. 分裂寻找阶段：类别如何按目标统计量均值排序后，仅考虑 2^(K-1)-1 种不相交子集划分 (而非有序阈值)

3. 缺失值在分类特征分裂中的走向决策 (基于增益最大化)

回答问题：

• 为什么原生分类特征支持在 max_depth 受限时比 Ordinal/One-Hot 编码更有优势？

• min_samples_leaf 参数如何约束分类特征的分裂？

70.23.3 实现一个简化版 StackingRegressor 并解读元权重

不使用 sklearn.ensemble.StackingRegressor，手动实现一个 3 层 Stacking 流程：

1. 基学习器：LinearRegression, DecisionTreeRegressor(max_depth=3), KNeighborsRegressor

2. 使用 sklearn.model_selection.cross_val_predict (cv=5) 获取每个基学习器的 Out-of-fold 预测作为元特征矩阵 (n_samples, 3)

3. 元学习器：RidgeCV() 在元特征上训练，对比直接在原始特征上训练的 RidgeCV 性能

4. 提取元学习器系数 coef_，分析哪个基学习器贡献最大

5. 替换元学习器为 LassoCV() 观察权重稀疏化效果

回答问题：

• 为什么必须使用 Out-of-fold 预测而非训练集预测作为元特征？如果用训练集预测会怎样？

• 元学习器权重为负值意味着什么？在什么业务场景下需要强制权重非负 (SuperLearner 约束)？

70.23.4 对比 IsolationForest 与 LocalOutlierFactor 在不同数据分布下的异常检测行为

使用 sklearn.datasets.make_blobs 生成三种合成数据集：

1. 单一高斯簇 + 少量均匀噪声 (全局异常)

2. 两个高斯簇密度差异大 + 桥接区域 (局部异常/密度变化)

3. 月牙形分布 make_moons + 噪声 (流形异常)

分别训练 IsolationForest(contamination=0.1) 与 LocalOutlierFactor(n_neighbors=20, contamination=0.1)

可视化对比：决策边界 (DecisionBoundaryDisplay)、异常分数分布直方图、Precision/Recall/F1 (若有真实标签)

回答问题：

• 为什么 IsolationForest 在全局异常 (数据集1) 上表现更好，而在局部密度变化 (数据集2) 上可能失效？

• LOF 的 n_neighbors 参数如何影响其对‘局部’的定义？在数据集3上调大/调小该参数会怎样？

• 两种方法的计算复杂度分别为何？大规模数据下哪个更适合？

70.23.5 复现偏差-方差分解实验并扩展到其他基学习器

基于 plot_bias_variance.py 的蒙特卡洛框架 (n_repeat=50, n_train=50, n_test=1000)，替换基学习器为：

1. KNeighborsRegressor(n_neighbors=1) (高方差、低偏差)

2. KNeighborsRegressor(n_neighbors=10) (低方差、高偏差)

3. LinearRegression() (高偏差、低方差，假设真实函数非线性)

4. BaggingRegressor(KNeighborsRegressor(n_neighbors=1)) (Bagging 降低 KNN 方差)

计算并绘制四种模型的 Bias², Variance, Noise, Total Error 分解曲线与均值

回答问题：

• Bagging 对高方差模型 (1-NN) 效果最显著，对高偏差模型 (LinearRegression) 效果为何有限？

• 为什么 n_neighbors=10 的 KNN 方差比 n_neighbors=1 低？偏差为何更高？

• 如果基学习器是 DecisionTreeRegressor(max_depth=2) (欠拟合)，Bagging 能否显著降低偏差？

70.24 设计取舍

为什么采用当前方案而不是更复杂的替代方案？ 本章实现优先保证与既有 API 的一致性、可维护性与运行效率。这意味着在少数极端场景下，调用者需要自行在灵活性、内存与速度之间做取舍，换取默认路径的清晰与稳定。

第 71 章 —— 线性模型精研 —— 回归与正则化的"数学手术刀"

71.1 学习目标

• 难度：★★★☆☆（3/5）

• 预备知识：Python 基础、面向对象编程与 Markdown/代码阅读基础

• 理解 Bagging、Boosting、Stacking、Voting 四大集成范式的核心原理与适用场景

• 掌握随机森林的 OOB 误差估计、特征重要性计算与偏差-方差分解

• 深入剖析梯度提升树 (GBDT) 的分阶段拟合、早停策略、分位数回归与正则化机制

• 理解直方图梯度提升 (HistGB) 的分箱策略、树生长算法与分类特征处理

• 掌握 Stacking 与 Voting 中异构模型融合的交叉验证元特征生成与加权投票机制

• 理解 IsolationForest 异常检测的隔离树构建、路径长度计算与异常分数聚合

• 掌握单调约束如何将领域知识注入梯度提升模型以保证预测的单调性

• 理解线性模型中 L1/L2 正则化路径、坐标下降、稀疏求解与弹性网络的数学原理

• 掌握鲁棒回归 (Huber、RANSAC、Theil-Sen) 与分位数回归抗离群点的核心机制

• 理解广义线性模型 (Poisson、Gamma、Tweedie) 与贝叶斯回归 (ARD、BayesianRidge) 的概率建模思想

• 掌握随机梯度下降 (SGD) 的损失函数、惩罚项、早停与加权样本在大规模学习中的应用

• 理解多项式/样条特征扩展、非负最小二乘 (NNLS) 与正交匹配追踪 (OMP) 在特征工程中的作用

71.2 生活类比

集成学习就像"联邦议会"的决策机制，多个基学习器（议员）通过不同方式协作决策：Bagging/随机森林像"平行宇宙投票制"——让每个议员在不同平行宇宙（自助采样）独立思考，最终少数服从多数；Boosting/GBDT像"接力赛纠错制"——第一个议员跑完留下差距（残差），下一个议员专门补齐这个差距，层层递进；Stacking像"专家咨询委员会"——多个领域专家（异构模型）各抒己见，再由元学习器（总理）综合裁决；Voting像"直接举手表决"——多个模型直接投票（硬投票）或平均概率（软投票），简单高效；IsolationForest像"孤岛猎人"——随机砍树（分裂特征），越容易被孤立的点（短路径）越可疑；单调约束像"给模型戴上镣铐"——告诉模型"房价必须随面积增加而上涨"，防止模型学歪。

线性模型就像"精密手术刀"，在高维特征空间中精准切除噪声：Lasso/L1像"极简主义雕塑家"——果断砍掉无关特征，保留核心骨架；Ridge/L2像"平滑压路机"——不删除特征，但压缩系数幅度，防止过拟合震荡；ElasticNet像"平衡大师"——兼具雕塑家的果断与压路机的稳健；鲁棒回归像"防弹背心"——面对离群值（子弹）依然稳立，不被极端样本带偏；分位数回归像"区间预言家"——不给单点答案，给"90%概率房价在300-500万之间"的区间；贝叶斯回归像"概率侧写师"——不给确定系数，给系数的"后验分布画像"，量化不确定性；SGD像"流式加工厂"——数据流过即处理，无需囤积全量数据，适合海量流式场景。就像医生开药需对症下药，建模也要根据数据特性（稀疏/稠密、有无离群值、样本量大小）选择合适的"线性手术刀"。

71.3 源码地图

examples/ensemble/

├── plot_adaboost_multiclass.py # AdaBoost多类分类与SAMME算法演示

├── plot_adaboost_regression.py # AdaBoost.R2回归变体演示

├── plot_adaboost_twoclass.py # AdaBoost二分类权重更新可视化

├── plot_gradient_boosting_categorical.py # HistGB分类特征原生处理对比

├── plot_gradient_boosting_early_stopping.py # GBDT早停策略与验证分数监控

├── plot_gradient_boosting_oob.py # GBDT袋外误差(OOB)估计与学习曲线

├── plot_gradient_boosting_quantile.py # GBDT分位数回归预测区间构建

├── plot_gradient_boosting_regression.py # GBDT回归核心流程与损失函数对比

├── plot_gradient_boosting_regularization.py # GBDT正则化(学习率、子采样、最大深度)效果

├── plot_hgbt_regression.py # HistGB回归与传统GBDT性能对比

├── plot_bias_variance.py # 集成学习偏差-方差分解可视化

├── plot_ensemble_oob.py # 随机森林/GBDT袋外(OOB)评分机制

├── plot_feature_transformation.py # 随机树嵌入(RandomTreesEmbedding)特征工程

├── plot_forest_hist_grad_boosting_comparison.py # 随机森林 vs HistGB 全维度对比

├── plot_forest_importances.py # 随机森林特征重要性(MDI/排列重要性)对比

├── plot_forest_iris.py # 随机森林决策边界可视化(Iris数据集)

├── plot_random_forest_embedding.py # 随机森林嵌入用于非线性降维

├── plot_random_forest_regression_multioutput.py # 多输出随机森林回归

├── plot_stack_predictors.py # StackingRegressor元学习器融合异构基模型

├── plot_voting_decision_regions.py # VotingClassifier硬/软投票决策边界对比

├── plot_voting_regressor.py # VotingRegressor加权平均集成回归

├── plot_isolation_forest.py # IsolationForest异常检测隔离树可视化

└── plot_monotonic_constraints.py # HistGB单调约束注入领域知识

examples/linear_model/

├── plot_lasso_and_elasticnet.py # Lasso/ARD/ElasticNet稀疏信号恢复对比

├── plot_lasso_dense_vs_sparse_data.py # Lasso稠密vs稀疏矩阵数值一致性与速度

├── plot_lasso_lars_ic.py # LassoLarsIC AIC/BIC模型选择

├── plot_lasso_lasso_lars_elasticnet_path.py # Lasso/LARS/ElasticNet正则化路径可视化

├── plot_lasso_model_selection.py # LassoCV/LassoLarsCV/LassoLarsIC三策略对比

├── plot_ridge_coeffs.py # Ridge系数收缩轨迹与MSE曲线

├── plot_ridge_path.py # 希尔伯特矩阵上Ridge正则化路径

├── plot_elastic_net_precomputed_gram_matrix_with_weighted_samples.py # 加权样本预计算Gram矩阵

├── plot_ols_ridge.py # OLS单特征拟合与高方差演示

├── plot_nnls.py # 非负最小二乘NNLS稀疏性对比

├── plot_omp.py # 正交匹配追踪稀疏信号恢复

├── plot_huber_vs_ridge.py # HuberRegressor vs Ridge抗离群值

├── plot_ransac.py # RANSAC鲁棒线性拟合内点识别

├── plot_theilsen.py # Theil-Sen中位数斜率鲁棒回归

├── plot_robust_fit.py # 四大鲁棒估计器多污染场景综合对比

├── plot_quantile_regression.py # 分位数回归异方差/重尾分布建模

├── plot_poisson_regression_non_normal_loss.py # 泊松回归保险理赔频率建模

├── plot_tweedie_regression_insurance_claims.py # Tweedie复合泊松-伽马纯保费建模

├── plot_ard.py # ARD vs BayesianRidge贝叶斯稀疏对比

├── plot_bayesian_ridge_curvefit.py # BayesianRidge超参数初值敏感性分析

├── plot_sgd_early_stopping.py # SGD三种早停策略MNIST实测

├── plot_sgd_iris.py # SGD鸢尾花OVA决策边界可视化

├── plot_sgd_loss_functions.py # SGD六大凸损失函数曲线对比

├── plot_sgd_penalties.py # L1/L2/ElasticNet惩罚等高线几何

├── plot_sgd_separating_hyperplane.py # SGD最大间隔超平面可视化

├── plot_sgd_weighted_samples.py # SGD加权样本决策边界偏移

├── plot_sgdocsvm_vs_ocsvm.py # SGDOneClassSVM核近似vs原生OCSVM

├── plot_logistic_l1_l2_sparsity.py # 逻辑回归L1/L2/ElasticNet稀疏度对比

├── plot_logistic_multinomial.py # 多项式vs OvR逻辑回归决策边界

├── plot_logistic_path.py # L1逻辑回归正则化路径

├── plot_sparse_logistic_regression_20newsgroups.py # 20newsgroups稀疏多项式vs OvR

├── plot_sparse_logistic_regression_mnist.py # MNIST稀疏多项式逻辑回归

├── plot_polynomial_interpolation.py # 多项式/样条/周期样条插值对比

└── plot_multi_task_lasso_support.py # 多任务Lasso联合特征选择

71.4 Boosting 家族 —— AdaBoost 与梯度提升的"迭代接力赛"

AdaBoost：自适应权重调优的分类与回归

• plot_adaboost_twoclass.py：二分类中错误样本权重指数级增加，后续弱学习器聚焦难例

• plot_adaboost_multiclass.py：SAMME 算法扩展至多类，需弱学习器准确率 > 1/K

• plot_adaboost_regression.py：AdaBoost.R2 根据相对误差调整样本权重，而非二分类的指数损失

传统 GBDT：分阶段拟合负梯度的核心循环

• plot_gradient_boosting_regression.py：GradientBoostingRegressor 支持 ls、lad、huber、quantile 四种损失

• plot_gradient_boosting_early_stopping.py：validation_fraction + n_iter_no_change 监控验证分数自动停止，防止过拟合

• plot_gradient_boosting_oob.py：subsample<1.0 时启用 OOB 评分，利用袋外样本估计泛化误差

GBDT 进阶：分位数回归与正则化工程化

• plot_gradient_boosting_quantile.py：loss='quantile', alpha=0.05/0.95 构建 90% 预测区间，捕捉异方差

• plot_gradient_boosting_regularization.py：learning_rate、subsample、max_depth、min_samples_leaf 四大正则化手段联合控制模型复杂度

直方图梯度提升 (HistGB)：现代高效 GBDT 引擎

• plot_hgbt_regression.py 与 plot_gradient_boosting_categorical.py：分箱离散化连续特征，原生支持分类特征 (无需 One-Hot)，训练速度提升 10-100 倍

• plot_forest_hist_grad_boosting_comparison.py：随机森林 vs HistGB 在准确率、训练速度、推理速度、内存占用全维度对比

我们从代码开始，先看 AdaBoost 的核心逻辑。下面代码来自 plot_adaboost_twoclass.py，演示了二分类 AdaBoost 的工作原理。

源码路径：examples/ensemble/plot_adaboost_twoclass.py - __main__（1-80行）

# 第 71 章 —— 构建数据集：两个高斯分布的叠加，形成非线性可分数据
X1, y1 = make_gaussian_quantiles(
    cov=2.0, n_samples=200, n_features=2, n_classes=2, random_state=1
)
X2, y2 = make_gaussian_quantiles(
    mean=(3, 3), cov=1.5, n_samples=300, n_features=2, n_classes=2, random_state=1
)
X = np.concatenate((X1, X2))
y = np.concatenate((y1, -y2 + 1))  # 第二类标签取反，使两类呈交错分布

# 第 71 章 —— 创建并训练 AdaBoost 分类器，基学习器为深度为1的决策树（决策桩）
bdt = AdaBoostClassifier(DecisionTreeClassifier(max_depth=1), n_estimators=200)
bdt.fit(X, y)

# 第 71 章 —— 绘制决策边界
plot_colors = "br"
plot_step = 0.02
class_names = "AB"

plt.figure(figsize=(10, 5))
ax = plt.subplot(121)
disp = DecisionBoundaryDisplay.from_estimator(
    bdt,
    X,
    cmap=plt.cm.Paired,
    response_method="predict",
    ax=ax,
    xlabel="x",
    ylabel="y",
)
x_min, x_max = disp.xx0.min(), disp.xx0.max()
y_min, y_max = disp.xx1.min(), disp.xx1.max()
plt.axis("tight")

# 第 71 章 —— 绘制训练点
for i, n, c in zip(range(2), class_names, plot_colors):
    idx = (y == i).nonzero()
    plt.scatter(
        X[idx, 0],
        X[idx, 1],
        c=c,
        s=20,
        edgecolor="k",
        label="Class %s" % n,
    )
plt.xlim(x_min, x_max)
plt.ylim(y_min, y_max)
plt.legend(loc="upper right")
plt.title("Decision Boundary")

# 第 71 章 —— 绘制决策得分分布直方图
twoclass_output = bdt.decision_function(X)
plot_range = (twoclass_output.min(), twoclass_output.max())
plt.subplot(122)
for i, n, c in zip(range(2), class_names, plot_colors):
    plt.hist(
        twoclass_output[y == i],
        bins=10,
        range=plot_range,
        facecolor=c,
        label="Class %s" % n,
        alpha=0.5,
        edgecolor="k",
    )
x1, x2, y1, y2 = plt.axis()
plt.axis((x1, x2, y1, y2 * 1.2))
plt.legend(loc="upper right")
plt.ylabel("Samples")
plt.xlabel("Score")
plt.title("Decision Scores")

plt.tight_layout()
plt.subplots_adjust(wspace=0.35)
plt.show()

这段代码定义了一个二分类 AdaBoost 分类器，可视化了其决策边界和决策得分分布。通过迭代调整样本权重，使得被错误分类的样本获得更高权重，后续弱学习器聚焦于难分样本，从而逐步提升模型性能。

继续探索梯度提升回归的核心机制。下面代码来自 plot_gradient_boosting_regression.py，演示了糖尿病数据集上的 GBDT 回归流程。

源码路径：examples/ensemble/plot_gradient_boosting_regression.py - __main__（1-120行）

# 第 71 章 —— 加载糖尿病数据集
diabetes = datasets.load_diabetes()
X, y = diabetes.data, diabetes.target

# 第 71 章 —— 数据预处理：划分训练集和测试集
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.1, random_state=13
)

# 第 71 章 —— 设置梯度提升回归器参数
params = {
    "n_estimators": 500,
    "max_depth": 4,
    "min_samples_split": 5,
    "learning_rate": 0.01,
    "loss": "squared_error",
}

# 第 71 章 —— 训练模型
reg = ensemble.GradientBoostingRegressor(**params)
reg.fit(X_train, y_train)

# 第 71 章 —— 计算测试集均方误差
mse = mean_squared_error(y_test, reg.predict(X_test))
print("The mean squared error (MSE) on test set: {:.4f}".format(mse))

# 第 71 章 —— 可视化训练和测试偏差随迭代的变化
test_score = np.zeros((params["n_estimators"],), dtype=np.float64)
for i, y_pred in enumerate(reg.staged_predict(X_test)):
    test_score[i] = mean_squared_error(y_test, y_pred)

fig = plt.figure(figsize=(6, 6))
plt.subplot(1, 1, 1)
plt.title("Deviance")
plt.plot(
    np.arange(params["n_estimators"]) + 1,
    reg.train_score_,
    "b-",
    label="Training Set Deviance",
)
plt.plot(
    np.arange(params["n_estimators"]) + 1, test_score, "r-", label="Test Set Deviance"
)
plt.legend(loc="upper right")
plt.xlabel("Boosting Iterations")
plt.ylabel("Deviance")
fig.tight_layout()
plt.show()

# 第 71 章 —— 绘制特征重要性（基于不纯度下降和置换重要性）
feature_importance = reg.feature_importances_
sorted_idx = np.argsort(feature_importance)
pos = np.arange(sorted_idx.shape[0]) + 0.5
fig = plt.figure(figsize=(12, 6))
plt.subplot(1, 2, 1)
plt.barh(pos, feature_importance[sorted_idx], align="center")
plt.yticks(pos, np.array(diabetes.feature_names)[sorted_idx])
plt.title("Feature Importance (MDI)")

result = permutation_importance(
    reg, X_test, y_test, n_repeats=10, random_state=42, n_jobs=2
)
sorted_idx = result.importances_mean.argsort()
plt.subplot(1, 2, 2)
tick_labels_parameter_name = (
    "tick_labels"
    if parse_version(matplotlib.__version__) >= parse_version("3.9")
    else "labels"
)
tick_labels_dict = {
    tick_labels_parameter_name: np.array(diabetes.feature_names)[sorted_idx]
}
plt.boxplot(result.importances[sorted_idx].T, vert=False, **tick_labels_dict)
plt.title("Permutation Importance (test set)")
fig.tight_layout()
plt.show()

这段代码实现了梯度提升回归器，可视化了训练和测试偏差随迭代的变化，并对比了基于不纯度下降（MDI）和置换重要性的特征重要性。模型通过迭代拟合残差来逐步改进预测，每棵树专注于前一棵树的误差。

现在我们看看直方图梯度提升如何通过分箱加速训练。下面代码来自 plot_hgbt_regression.py，展示了 HistGB 在澳大利亚电力数据集上的应用。

源码路径：examples/ensemble/plot_hgbt_regression.py - __main__（1-100行）

# 第 71 章 —— 加载澳大利亚电力数据集
electricity = fetch_openml(
    name="electricity", version=1, as_frame=True, parser="pandas"
)
df = electricity.frame

# 第 71 章 —— 移除前17,760个样本（步长常量目标）
df = electricity.frame.iloc[17_760:]
X = df.drop(columns=["transfer", "class"])
y = df["transfer"]

# 第 71 章 —— 划分训练集和测试集（按时间顺序，不洗牌）
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]}")

# 第 71 章 —— 训练不同树数量的 HistGB 模型以观察收敛行为
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()

# 第 71 章 —— 启用早停机制自动确定最优迭代次数
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_})",
)

# 第 71 章 —— 根据早停结果调整 max_iter 以避免过度计算
import math
common_params["max_iter"] = math.ceil(hgbt.n_iter_ / 100) * 100
common_params["early_stopping"] = False
hgbt = HistGradientBoostingRegressor(**common_params)

这段代码演示了 HistGB 模型的训练过程，包括如何通过早停机制自动确定最优迭代次数，以及如何处理缺失值和分类特征。模型通过分箱离散化连续特征并原生支持分类特征，显著提升了训练速度。

接下来，我们通过时序图进一步理解 Boosting 家族的工作机制。

sequenceDiagram participant AdaBoost participant GBDT participant HistGB AdaBoost->>AdaBoost: 初始化样本权重均匀 AdaBoost->>AdaBoost: 训练弱学习器 AdaBoost->>AdaBoost: 根据错误增加难分样本权重 AdaBoost->>AdaBoost: 计算弱学习器权重（准确率越高权重越大） AdaBoost->>AdaBoost: 加权组合弱学习器输出 GBDT->>GBDT: 初始化预测为0 GBDT->>GBDT: 计算当前预测的残差（负梯度） GBDT->>GBDT: 拟合弱学习器到残差上 GBDT->>GBDT: 更新预测 += 学习率 × 弱学习器输出 HistGB->>HistGB: 对连续特征进行分箱离散化 HistGB->>HistGB: 构建直方图以加速分裂点搜索 HistGB->>HistGB: 原生处理分类特征和缺失值 HistGB->>HistGB: 使用梯度提升框架迭代优化

通过上述分析，我们理解了 Boosting 家族如何通过迭代纠错机制提升模型性能，其中 HistGB 通过分箱和原生特征处理显著提升了效率。

71.5 随机森林与装袋 —— 多样性的"群体智慧"

随机森林核心：Bootstrap + 特征随机化 + 并行训练

• plot_forest_iris.py：决策边界可视化，展示集成如何平滑单棵树的阶跃边界

• plot_forest_importances.py：MDI (Mean Decrease Impurity) vs 排列重要性，前者偏向高基数特征，后者更可靠

• plot_ensemble_oob.py：oob_score=True 利用 ~37% 袋外样本计算无偏泛化误差，无需额外验证集

偏差-方差分解：集成学习为何有效的数学本质

• plot_bias_variance.py：单棵树高方差 (过拟合) → Bagging 平均降方差 → Boosting 串行降偏差

• 随机森林 = Bagging (降方差) + 弱学习器加深 (适度降偏差)，达到 Bias-Variance 最优平衡

随机森林进阶应用：嵌入、多输出与特征变换

• plot_random_forest_embedding.py：RandomTreesEmbedding 将树叶子节点索引作为高维稀疏特征，配合线性模型处理非线性

• plot_random_forest_regression_multioutput.py：单模型同时预测多个目标变量，捕捉目标间相关性

• plot_feature_transformation.py：树嵌入作为非线性特征工程，替代 PolynomialFeatures/Kernel 方法

我们从代码开始，先看随机森林的 OOB 误差估计机制。下面代码来自 plot_ensemble_oob.py。

源码路径：examples/ensemble/plot_ensemble_oob.py - __main__（1-100行）

# 第 71 章 —— 生成二分类数据集
X, y = make_classification(
    n_samples=500,
    n_features=25,
    n_clusters_per_class=1,
    n_informative=15,
    random_state=RANDOM_STATE,
)

# 第 71 章 —— 定义不同 max_features 设置下的随机森林分类器（启用 warm_start 以追踪 OOB 错误）
ensemble_clfs = [
    (
        "RandomForestClassifier, max_features='sqrt'",
        RandomForestClassifier(
            warm_start=True,
            oob_score=True,
            max_features="sqrt",
            random_state=RANDOM_STATE,
        ),
    ),
    (
        "RandomForestClassifier, max_features='log2'",
        RandomForestClassifier(
            warm_start=True,
            max_features="log2",
            oob_score=True,
            random_state=RANDOM_STATE,
        ),
    ),
    (
        "RandomForestClassifier, max_features=None",
        RandomForestClassifier(
            warm_start=True,
            max_features=None,
            oob_score=True,
            random_state=RANDOM_STATE,
        ),
    ),
]

# 第 71 章 —— 记录不同树数量下的 OOB 错误率
error_rate = OrderedDict((label, []) for label, _ in ensemble_clfs)
min_estimators = 15
max_estimators = 150

for label, clf in ensemble_clfs:
    for i in range(min_estimators, max_estimators + 1, 5):
        clf.set_params(n_estimators=i)
        clf.fit(X, y)

        # 记录 OOB 错误率：1 - OOB 分数
        oob_error = 1 - clf.oob_score_
        error_rate[label].append((i, oob_error))

# 第 71 章 —— 绘制 OOB 错误率随树数量的变化
for label, clf_err in error_rate.items():
    xs, ys = zip(*clf_err)
    plt.plot(xs, ys, label=label)

plt.xlim(min_estimators, max_estimators)
plt.xlabel("n_estimators")
plt.ylabel("OOB error rate")
plt.legend(loc="upper right")
plt.show()

这段代码演示了如何通过 OOB（Out-of-Bag）误差估计随机森林的泛化性能。OOB 误差利用每棵树在自助采样过程中未被选中的样本（约37%）进行评估，无需额外验证集即可获得无偏泛化误差估计。

现在我们看看随机森林如何通过偏差-方差分解解释其有效性。下面代码来自 plot_bias_variance.py。

源码路径：examples/ensemble/plot_bias_variance.py - __main__（1-100行）

# 第 71 章 —— 设置实验参数
n_repeat = 50  # 重复实验次数以计算期望
n_train = 50   # 训练集大小
n_test = 1000  # 测试集大小
noise = 0.1    # 噪声标准差
np.random.seed(0)

# 第 71 章 —— 比较单棵决策树和装袋决策树的偏差-方差分解
estimators = [
    ("Tree", DecisionTreeRegressor()),
    ("Bagging(Tree)", BaggingRegressor(DecisionTreeRegressor())),
]

n_estimators = len(estimators)


# 第 71 章 —— 生成数据函数：f(x) = exp(-x²) + 1.5·exp(-(x-2)²)
def f(x):
    x = x.ravel()
    return np.exp(-(x**2)) + 1.5 * np.exp(-((x - 2) ** 2))


def generate(n_samples, noise, n_repeat=1):
    X = np.random.rand(n_samples) * 10 - 5
    X = np.sort(X)

    if n_repeat == 1:
        y = f(X) + np.random.normal(0.0, noise, n_samples)
    else:
        y = np.zeros((n_samples, n_repeat))

        for i in range(n_repeat):
            y[:, i] = f(X) + np.random.normal(0.0, noise, n_samples)

    X = X.reshape((n_samples, 1))
    return X, y


# 第 71 章 —— 生成训练集和测试集
X_train = []
y_train = []

for i in range(n_repeat):
    X, y = generate(n_samples=n_train, noise=noise)
    X_train.append(X)
    y_train.append(y)

X_test, y_test = generate(n_samples=n_test, noise=noise, n_repeat=n_repeat)

plt.figure(figsize=(10, 8))

# 第 71 章 —— 循环比较每个估计器
for n, (name, estimator) in enumerate(estimators):
    # 计算预测
    y_predict = np.zeros((n_test, n_repeat))

    for i in range(n_repeat):
        estimator.fit(X_train[i], y_train[i])
        y_predict[:, i] = estimator.predict(X_test)

    # 均方误差分解：误差 = 偏差² + 方差 + 噪声
    y_error = np.zeros(n_test)

    for i in range(n_repeat):
        for j in range(n_repeat):
            y_error += (y_test[:, j] - y_predict[:, i]) ** 2

    y_error /= n_repeat * n_repeat

    y_noise = np.var(y_test, axis=1)
    y_bias = (f(X_test) - np.mean(y_predict, axis=1)) ** 2
    y_var = np.var(y_predict, axis=1)

    print(
        "{0}: {1:.4f} (error) = {2:.4f} (bias^2) "
        " + {3:.4f} (var) + {4:.4f} (noise)".format(
            name, np.mean(y_error), np.mean(y_bias), np.mean(y_var), np.mean(y_noise)
        )
    )

    # 绘图：左上角显示预测，右下角显示误差分解
    plt.subplot(2, n_estimators, n + 1)
    plt.plot(X_test, f(X_test), "b", label="$f(x)$")
    plt.plot(X_train[0], y_train[0], ".b", label="LS ~ $y = f(x)+noise$")

    for i in range(n_repeat):
        if i == 0:
            plt.plot(X_test, y_predict[:, i], "r", label=r"$\^y(x)$")
        else:
            plt.plot(X_test, y_predict[:, i], "r", alpha=0.05)

    plt.plot(X_test, np.mean(y_predict, axis=1), "c", label=r"$\mathbb{E}_{LS} \^y(x)$")

    plt.xlim([-5, 5])
    plt.title(name)

    if n == n_estimators - 1:
        plt.legend(loc=(1.1, 0.5))

    plt.subplot(2, n_estimators, n_estimators + n + 1)
    plt.plot(X_test, y_error, "r", label="$error(x)$")
    plt.plot(X_test, y_bias, "b", label="$bias^2(x)$")
    plt.plot(X_test, y_var, "g", label="$variance(x)$")
    plt.plot(X_test, y_noise, "c", label="$noise(x)$")

    plt.xlim([-5, 5])
    plt.ylim([0, 0.1])

    if n == n_estimators - 1:
        plt.legend(loc=(1.1, 0.5))

plt.subplots_adjust(right=0.75)
plt.show()

这段代码通过实验验证了偏差-方差分解理论：单棵决策树具有低偏差但高方差（过拟合），而随机森林通过平均多棵自助采样树显著降低方差，同时仅略微增加偏差，从而降低总体误差。

接下来，我们通过架构分层图进一步理解随机森林的工作机制。

graph TD A[训练数据] --> B[自助采样生成多个子集] B --> C[每个子集训练一棵决策树] C --> D[特征随机化：每个节点只考虑部分特征] D --> E[并行训练所有决策树] E --> F[预测阶段：多数投票/平均] F --> G[最终输出] style B fill:#f9f,stroke:#333 style D fill:#9f9,stroke:#333

通过上述分析，我们理解了随机森林如何通过自助采样和特征随机化引入多样性，从而在保持偏差可控的同时显著降低方差。

71.6 堆叠与投票 —— 异构模型的"联合执政"

Stacking：两层架构的元学习融合

• plot_stack_predictors.py：一级基学习器 (RF, SVM, Ridge等) 通过 CV 生成元特征，二级元学习器 (RidgeCV) 学习最优权重

• 关键点：cv=5 避免数据泄露，passthrough=True 保留原始特征，元学习器通常选简单线性模型防过拟合

Voting：无需训练的直接民主制

• plot_voting_decision_regions.py：硬投票 (类别众数) vs 软投票 (概率平均)，软投票利用置信度信息边界更平滑

• plot_voting_regressor.py：VotingRegressor 加权平均多个回归器预测，权重可基于 CV 分数设定

对比启示：何时用 Stacking、何时用 Voting

• Stacking 适合基学习器异构、性能差异大、有足够数据训练元学习器的场景

• Voting 适合基学习器同构/性能接近、追求极简部署、需在线推理低延迟的场景

我们从代码开始，先看 Stacking 的工作机制。下面代码来自 plot_stack_predictors.py。

源码路径：examples/ensemble/plot_stack_predictors.py - __main__（1-120行）

# 第 71 章 —— 生成带有突变的合成数据：正弦波 + 线性趋势 + 突然下降 + 异方差噪声
rng = np.random.RandomState(42)
X = rng.uniform(-3, 3, size=500)
trend = 2.4 * X
seasonal = 3.1 * np.sin(3.2 * X)
drop = 10.0 * (X > 2).astype(float)
sigma = 0.75 + 0.75 * X**2
y = trend + seasonal - drop + rng.normal(loc=0.0, scale=np.sqrt(sigma))

df = pd.DataFrame({"X": X, "y": y})
_ = df.plot.scatter(x="X", y="y")

# 第 71 章 —— 定义基学习器管道
linear_ridge = make_pipeline(StandardScaler(), RidgeCV())
spline_ridge = make_pipeline(
    SplineTransformer(n_knots=6, degree=3),
    PolynomialFeatures(interaction_only=True),
    RidgeCV(),
)
hgbt = HistGradientBoostingRegressor(random_state=0)

estimators = [
    ("Linear Ridge", linear_ridge),
    ("Spline Ridge", spline_ridge),
    ("HGBT", hgbt),
]

# 第 71 章 —— 创建 Stacking 回归器，使用 RidgeCV 作为元学习器
stacking_regressor = StackingRegressor(estimators=estimators, final_estimator=RidgeCV())
stacking_regressor

# 第 71 章 —— 训练所有模型并进行交叉验证评估
X = X.reshape(-1, 1)
linear_ridge.fit(X, y)
spline_ridge.fit(X, y)
hgbt.fit(X, y)
stacking_regressor.fit(X, y)

x_plot = np.linspace(X.min() - 0.1, X.max() + 0.1, 500).reshape(-1, 1)
preds = {
    "Linear Ridge": linear_ridge.predict(x_plot),
    "Spline Ridge": spline_ridge.predict(x_plot),
    "HGBT": hgbt.predict(x_plot),
    "Stacking (Ridge final estimator)": stacking_regressor.predict(x_plot),
}

fig, axes = plt.subplots(2, 2, figsize=(10, 8), sharex=True, sharey=True)
axes = axes.ravel()
for ax, (name, y_pred) in zip(axes, preds.items()):
    ax.scatter(
        X[:, 0],
        y,
        s=6,
        alpha=0.35,
        linewidths=0,
        label="observed (sample)",
    )

    ax.plot(x_plot.ravel(), y_pred, linewidth=2, alpha=0.9, label=name)
    ax.set_title(name)
    ax.set_xlabel("x")
    ax.set_ylabel("y")
    ax.legend(loc="lower right")

plt.suptitle("Base Models Predictions versus Stacked Predictions", y=1)
plt.tight_layout()
plt.show()

# 第 71 章 —— 评估预测误差
fig, axs = plt.subplots(2, 2, figsize=(9, 7))
axs = np.ravel(axs)

for ax, (name, est) in zip(
    axs, estimators + [("Stacking Regressor", stacking_regressor)]
):
    scorers = {r"$R^2$": "r2", "MAE": "neg_mean_absolute_error"}

    start_time = time.time()
    scores = cross_validate(est, X, y, scoring=list(scorers.values()), n_jobs=-1)
    elapsed_time = time.time() - start_time

    y_pred = cross_val_predict(est, X, y, n_jobs=-1)
    scores = {
        key: (
            f"{np.abs(np.mean(scores[f'test_{value}'])):.2f}"
            r" $\pm$ "
            f"{np.std(scores[f'test_{value}']):.2f}"
        )
        for key, value in scorers.items()
    }

    display = PredictionErrorDisplay.from_predictions(
        y_true=y,
        y_pred=y_pred,
        kind="actual_vs_predicted",
        ax=ax,
        scatter_kwargs={"alpha": 0.2, "color": "tab:blue"},
        line_kwargs={"color": "tab:red"},
    )
    ax.set_title(f"{name}\nEvaluation in {elapsed_time:.2f} seconds")

    for name, score in scores.items():
        ax.plot([], [], " ", label=f"{name}: {score}")
    ax.legend(loc="upper left")

plt.suptitle("Prediction Errors of Base versus Stacked Predictors", y=1)
plt.tight_layout()
plt.subplots_adjust(top=0.9)
plt.show()

# 第 71 章 —— 检查元学习器的系数（即基学习器的权重）
stacking_regressor.fit(X, y)
stacking_regressor.final_estimator_.coef_

这段代码演示了 Stacking 的两层架构：一级基学习器（线性ridge、样条ridge、HistGB）通过交叉验证生成元特征，二级元学习器（RidgeCV）学习如何最优组合这些预测。可以看到 HistGB 模型在最终集成中贡献最大。

现在我们看看 Voting 的工作机制。下面代码来自 plot_voting_decision_regions.py。

源码路径：examples/ensemble/plot_voting_decision_regions.py - __main__（1-80行）

# 第 71 章 —— 生成 noisy XOR 数据集：二分类非线性可分问题
n_samples = 500
rng = np.random.default_rng(0)
feature_names = ["Feature #0", "Feature #1"]
common_scatter_plot_params = dict(
    cmap=ListedColormap(["tab:red", "tab:blue"]),
    edgecolor="white",
    linewidth=1,
)

xor = pd.DataFrame(
    np.random.RandomState(0).uniform(low=-1, high=1, size=(n_samples, 2)),
    columns=feature_names,
)
noise = rng.normal(loc=0, scale=0.1, size=(n_samples, 2))
target_xor = np.logical_xor(
    xor["Feature #0"] + noise[:, 0] > 0, xor["Feature #1"] + noise[:, 1] > 0
)

X = xor[feature_names]
y = target_xor.astype(np.int32)

fig, ax = plt.subplots()
ax.scatter(X["Feature #0"], X["Feature #1"], c=y, **common_scatter_plot_params)
ax.set_title("The XOR dataset")
plt.show()

# 第 71 章 —— 定义三个基分类器管道
clf1 = make_pipeline(
    SplineTransformer(degree=2, n_knots=2),
    PolynomialFeatures(interaction_only=True),
    LogisticRegression(C=10),
)
clf2 = make_pipeline(
    SplineTransformer(
        degree=2,
        n_knots=4,
        extrapolation="periodic",
        include_bias=True,
    ),
    PolynomialFeatures(interaction_only=True),
    LogisticRegression(C=10),
)
clf3 = make_pipeline(
    StandardScaler(),
    Nystroem(gamma=2, random_state=0),
    LogisticRegression(C=10),
)
weights = [2, 1, 3]
eclf = VotingClassifier(
    estimators=[
        ("constant splines model", clf1),
        ("periodic splines model", clf2),
        ("nystroem model", clf3),
    ],
    voting="soft",
    weights=weights,
)

# 第 71 章 —— 训练所有分类器
clf1.fit(X, y)
clf2.fit(X, y)
clf3.fit(X, y)
eclf.fit(X, y)

# 第 71 章 —— 使用 DecisionBoundaryDisplay 可视化预测概率
from itertools import product

from sklearn.inspection import DecisionBoundaryDisplay

fig, axarr = plt.subplots(2, 2, sharex="col", sharey="row", figsize=(10, 8))
for idx, clf, title in zip(
    product([0, 1], [0, 1]),
    [clf1, clf2, clf3, eclf],
    [
        "Splines with\nconstant extrapolation",
        "Splines with\nperiodic extrapolation",
        "RBF Nystroem",
        "Soft Voting",
    ],
):
    disp = DecisionBoundaryDisplay.from_estimator(
        clf,
        X,
        response_method="predict_proba",
        plot_method="pcolormesh",
        cmap="RdBu",
        alpha=0.8,
        ax=axarr[idx[0], idx[1]],
    )
    axarr[idx[0], idx[1]].scatter(
        X["Feature #0"],
        X["Feature #1"],
        c=y,
        **common_scatter_plot_params,
    )
    axarr[idx[0], idx[1]].set_title(title)
    fig.colorbar(disp.surface_, ax=axarr[idx[0], idx[1]], label="Probability estimate")

plt.show()

# 第 71 章 —— 验证软投票是加权平均
test_sample = pd.DataFrame({"Feature #0": [-0.5], "Feature #1": [1.5]})
predict_probas = [est.predict_proba(test_sample).ravel() for est in eclf.estimators_]
for (est_name, _), est_probas in zip(eclf.estimators, predict_probas):
    print(f"{est_name}'s predicted probabilities: {est_probas}")

print(
    "Weighted average of soft-predictions: "
    f"{np.dot(weights, predict_probas) / np.sum(weights)}"
)

print(
    "Predicted probability of VotingClassifier: "
    f"{eclf.predict_proba(test_sample).ravel()}"
)

print(
    "Class with the highest weighted average of soft-predictions: "
    f"{np.argmax(np.dot(weights, predict_probas) / np.sum(weights))}"
)

print(f"Predicted class of VotingClassifier: {eclf.predict(test_sample).ravel()}")

这段代码演示了 Voting 机制：多个分类器（常数外推样条、周期外推样条、Nystroem RBF）通过软投票（加权平均概率）结合预测。可以看到软投票相比硬投票能够利用置信度信息，产生更平滑的决策边界。

接下来，我们通过数据流图进一步理解 Stacking 和 Voting 的工作机制。

flowchart TD A[输入数据] --> B{集成方式} B -->|Stacking| C[一级基学习器] C --> D[交叉验证生成元特征] D --> E[二级元学习器] E --> F[最终预测] B -->|Voting| G[所有基学习器并行预测] G --> H[硬投票：类别众数<br/>软投票：概率平均] H --> I[最终预测] style C fill:#f9f,stroke:#333 style E fill:#9f9,stroke:#333 style G fill:#ff9,stroke:#333

通过上述分析，我们理解了 Stacking 通过交叉验证避免数据泄露并学习基学习器之间的最优组合，而 Voting 则提供了一种简单高效的无需额外训练的模型融合方式。

71.7 异常检测与单调约束 —— 集成的"边界守护者"

IsolationForest：基于路径长度的异常隔离机制

• plot_isolation_forest.py：随机选择特征+分裂值构建隔离树，异常样本因易被孤立而拥有更短平均路径长度

• n_estimators 森林聚合降低方差，contamination 控制异常判定阈值，max_samples 限制树规模加速训练

• 无监督、线性时间复杂度、适合高维数据，但对密度变化敏感、全局异常检测优于局部

单调约束：将领域知识硬编码进树分裂逻辑

• plot_monotonic_constraints.py：monotonic_cst=[1, -1, 0] 强制特征 0 单增、特征 1 单减、特征 2 无约束

• 实现原理：分裂增益计算时过滤违反单调性的分裂点，叶子值约束满足单调性

• 典型场景：信用评分 (收入↑违约率↓)、保险定价 (车龄↑保费↑)、医疗风险 (年龄↑风险↑)

工程权衡：约束 vs 灵活性

• 单调约束可能导致欠拟合 (无法拟合真实非单调关系) → 可放宽约束或仅对核心特征施加

• IsolationForest 在数据存在明显密度簇时易误报 → 可结合 LOF/OneClassSVM 多模型集成

我们从代码开始，先看 IsolationForest 的工作机制。下面代码来自 plot_isolation_forest.py。

源码路径：examples/ensemble/plot_isolation_forest.py - __main__（1-100行）

# 第 71 章 —— 生成数据：两个高斯簇作为内点，均匀分布作为离群点
n_samples, n_outliers = 120, 40
rng = np.random.RandomState(0)
covariance = np.array([[0.5, -0.1], [0.7, 0.4]])
cluster_1 = 0.4 * rng.randn(n_samples, 2) @ covariance + np.array([2, 2])  # general
cluster_2 = 0.3 * rng.randn(n_samples, 2) + np.array([-2, -2])  # spherical
outliers = rng.uniform(low=-4, high=4, size=(n_outliers, 2))

X = np.concatenate([cluster_1, cluster_2, outliers])
y = np.concatenate(
    [np.ones((2 * n_samples), dtype=int), -np.ones((n_outliers), dtype=int)]
)

X_train, X_test, y_train, y_test = train_test_split(X, y, stratify=y, random_state=42)

# 第 71 章 —— 可视化数据分布
scatter = plt.scatter(X[:, 0], X[:, 1], c=y, s=20, edgecolor="k")
handles, labels = scatter.legend_elements()
plt.axis("square")
plt.legend(handles=handles, labels=["outliers", "inliers"], title="true class")
plt.title("Gaussian inliers with \nuniformly distributed outliers")
plt.show()

# 第 71 章 —— 训练 Isolation Forest 模型
clf = IsolationForest(max_samples=100, random_state=0)
clf.fit(X_train)

# 第 71 章 —— 可视化离散决策边界（背景色表示是否被预测为离群点）
disp = DecisionBoundaryDisplay.from_estimator(
    clf,
    X,
    response_method="predict",
    alpha=0.5,
)
disp.ax_.scatter(X[:, 0], X[:, 1], c=y, s=20, edgecolor="k")
disp.ax_.set_title("Binary decision boundary \nof IsolationForest")
plt.axis("square")
plt.legend(handles=handles, labels=["outliers", "inliers"], title="true class")
plt.show()

# 第 71 章 —— 可视化路径长度决策边界（背景色表示异常分数）
disp = DecisionBoundaryDisplay.from_estimator(
    clf,
    X,
    response_method="decision_function",
    alpha=0.5,
)
disp.ax_.scatter(X[:, 0], X[:, 1], c=y, s=20, edgecolor="k")
disp.ax_.set_title("Path length decision boundary \nof IsolationForest")
plt.axis("square")
plt.legend(handles=handles, labels=["outliers", "inliers"], title="true class")
plt.colorbar(disp.ax_.collections[1])
plt.show()

这段代码演示了 Isolation Forest 的工作原理：通过随机选择特征和分裂值构建隔离树，容易被隔离（路径短）的样本被视为异常点。模型无需标签即可工作，因为其核心假设是异常点较少且属性值异构，因而更易被随机分割隔离。

现在我们看看单调约束如何注入领域知识。下面代码来自 plot_monotonic_constraints.py。

源码路径：examples/ensemble/plot_monotonic_constraints.py - __main__（1-100行）

# 第 71 章 —— 生成人工数据：特征0正相关，特征1负相关，带噪声
rng = np.random.RandomState(0)

n_samples = 1000
f_0 = rng.rand(n_samples)
f_1 = rng.rand(n_samples)
X = np.c_[f_0, f_1]
noise = rng.normal(loc=0.0, scale=0.01, size=n_samples)

# 第 71 章 —— y 正相关于 f_0，负相关于 f_1
y = 5 * f_0 + np.sin(10 * np.pi * f_0) - 5 * f_1 - np.cos(10 * np.pi * f_1) + noise


# 第 71 章 —— 训练无约束模型
gbdt_no_cst = HistGradientBoostingRegressor()
gbdt_no_cst.fit(X, y)

# 第 71 章 —— 训练带单调约束的模型：特征0单增，特征1单减
gbdt_with_monotonic_cst = HistGradientBoostingRegressor(monotonic_cst=[1, -1])
gbdt_with_monotonic_cst.fit(X, y)


# 第 71 章 —— 可视化特征对目标的部分依赖
fig, ax = plt.subplots()
disp = PartialDependenceDisplay.from_estimator(
    gbdt_no_cst,
    X,
    features=[0, 1],
    feature_names=(
        "First feature",
        "Second feature",
    ),
    line_kw={"linewidth": 4, "label": "unconstrained", "color": "tab:blue"},
    ax=ax,
)
PartialDependenceDisplay.from_estimator(
    gbdt_with_monotonic_cst,
    X,
    features=[0, 1],
    line_kw={"linewidth": 4, "label": "constrained", "color": "tab:orange"},
    ax=disp.axes_,
)

for f_idx in (0, 1):
    disp.axes_[0, f_idx].plot(
        X[:, f_idx], y, "o", alpha=0.3, zorder=-1, color="tab:green"
    )
    disp.axes_[0, f_idx].set_ylim(-6, 6)

plt.legend()
fig.suptitle("Monotonic constraints effect on partial dependences")
plt.show()