Sklearn-源码解析-书-v1-0-十四-
Sklearn 源码解析（书）v1.0（十四）
（这段代码定义了 Ledoit-Wolf 收缩估计器的公共接口：它内部创建一个 LedoitWolf 估计器实例，调用其 fit 方法，然后返回收缩后的协方差矩阵和收缩系数。）
源码路径：sklearn/covariance/_graph_lasso.py - GraphicalLasso.fit()（200-350行）
// 逐行注释解释
@_fit_context(prefer_skip_nested_validation=True)
def fit(self, X, y=None):
"""Fit the GraphicalLasso model to X.
Parameters
----------
X : array-like of shape (n_samples, n_features)
Data from which to compute the covariance estimate.
y : Ignored
Not used, present for API consistency by convention.
Returns
-------
self : object
Returns the instance itself.
"""
# Covariance does not make sense for a single feature
X = validate_data(self, X, ensure_min_features=2, ensure_min_samples=2)
if self.covariance == "precomputed":
emp_cov = X.copy()
self.location_ = np.zeros(X.shape[1])
else:
emp_cov = empirical_covariance(X, assume_centered=self.assume_centered)
if self.assume_centered:
self.location_ = np.zeros(X.shape[1])
else:
self.location_ = X.mean(0)
self.covariance_, self.precision_, self.costs_, self.n_iter_ = _graphical_lasso(
emp_cov,
alpha=self.alpha,
cov_init=None,
mode=self.mode,
tol=self.tol,
enet_tol=self.enet_tol,
max_iter=self.max_iter,
verbose=self.verbose,
eps=self.eps,
)
return self
（这段代码实现了 GraphicalLasso 的 fit 方法：先验证输入数据（确保至少两个特征和两个样本），然后根据是否使用预计算协方差来决定输入是否为协方差矩阵；接着计算经验协方差（如果需要）；最后调用底层的 _graphical_lasso 函数（实现坐标下降或 LARS 算法）来求解 L1 正则化的精度矩阵估计问题。）
源码路径：sklearn/covariance/_robust_covariance.py - MinCovDet.fit()（300-450行）
// 逐行注释解释
@_fit_context(prefer_skip_nested_validation=True)
def fit(self, X, y=None):
"""Fit a Minimum Covariance Determinant with the FastMCD algorithm.
Parameters
----------
X : array-like of shape (n_samples, n_features)
Training data, where `n_samples` is the number of samples
and `n_features` is the number of features.
y : Ignored
Not used, present for API consistency by convention.
Returns
-------
self : object
Returns the instance itself.
"""
X = validate_data(self, X, ensure_min_samples=2, estimator="MinCovDet")
random_state = check_random_state(self.random_state)
n_samples, n_features = X.shape
# check that the empirical covariance is full rank
if (linalg.svdvals(np.dot(X.T, X)) > 1e-8).sum() != n_features:
warnings.warn(
"The covariance matrix associated to your dataset is not full rank"
)
# compute and store raw estimates
raw_location, raw_covariance, raw_support, raw_dist = fast_mcd(
X,
support_fraction=self.support_fraction,
cov_computation_method=self._nonrobust_covariance,
random_state=random_state,
)
if self.assume_centered:
raw_location = np.zeros(n_features)
raw_covariance = self._nonrobust_covariance(
X[raw_support], assume_centered=True
)
# get precision matrix in an optimized way
precision = linalg.pinvh(raw_covariance)
raw_dist = np.sum(np.dot(X, precision) * X, 1)
self.raw_location_ = raw_location
self.raw_covariance_ = raw_covariance
self.raw_support_ = raw_support
self.location_ = raw_location
self.support_ = raw_support
self.dist_ = raw_dist
# obtain consistency at normal models
self.correct_covariance(X)
# re-weight estimator
self.reweight_covariance(X)
return self
（这段代码实现了 MinCovDet 的 fit 方法：先验证输入数据，然后调用 fast_mcd 函数计算原始的稳健估计（位置、协方差、支持掩码和马氏距离）；如果假设数据以零中心，则重新计算以零为中心的协方差；接着应用一致性校正因子（使得在高斯假设下协方差无偏）；最后进行重新加权步骤（基于卡方分位数，进一步降低离群点的影响），得到最终的稳健位置和协方差估计。）
源码路径：sklearn/covariance/_elliptic_envelope.py - EllipticEnvelope.fit()（100-200行）
// 逐行注释解释
@_fit_context(prefer_skip_nested_validation=True)
def fit(self, X, y=None):
"""Fit the EllipticEnvelope model.
Parameters
----------
X : array-like of shape (n_samples, n_features)
Training data.
y : Ignored
Not used, present for API consistency by convention.
Returns
-------
self : object
Returns the instance itself.
"""
super().fit(X)
self.offset_ = np.percentile(-self.dist_, 100.0 * self.contamination)
return self
（这段代码实现了 EllipticEnvelope 的 fit 方法：先调用父类 MinCovDet 的 fit 方法来获取稳健位置、协方差和马氏距离；然后根据假设的污染比例（contamination）计算一个偏移量，使得决策函数为负的样本比例恰好等于该污染比例——也就是说，我们将马氏距离取负后，再减去这个偏移量，使得阈值为零。）
源码路径：sklearn/covariance/_elliptic_envelope.py - EllipticEnvelope.decision_function()（200-250行）
// 逐行注释解释
def decision_function(self, X):
"""Compute the decision function of the given observations.
Parameters
----------
X : array-like of shape (n_samples, n_features)
The data matrix.
Returns
-------
decision : ndarray of shape (n_samples,)
Decision function of the samples.
It is equal to the shifted Mahalanobis distances.
The threshold for being an outlier is 0, which ensures a
compatibility with other outlier detection algorithms.
"""
check_is_fitted(self)
negative_mahal_dist = self.score_samples(X)
return negative_mahal_dist - self.offset_
（这段代码实现了 decision_function 方法：先检查模型是否已拟合，然后调用 score_samples 方法得到负的马氏距离（即 -D²），再减去之前在 fit 中计算的偏移量，得到最终的决策函数值。当这个值小于 0 时，样本被判定为异常。）
27.6.2 完整数据流/流程图
27.7 设计中的取舍
-
为什么不用全协方差在所有情况下？
全协方差虽然最灵活，能够捕捉任意形状的高斯分量，但它参数众多（每个分量有 d(d+1)/2 个自由度），在高维或样本不足时易导过拟合和数值不稳（协方差矩阵奇异）。因此我们提供了 tied、diag、spherical 三种受限形式，以换取更好的泛化能力和数值稳定性，特别是在实际应用中特征往往存在冗余或近似独立的情况下。
-
贝叶斯 GMM 与标准 GMM 在组件数选择上的 trade-off 是什么？
标准 GMM 需要依赖准则如 BIC/AIC 来事后选择组件数，这需要反复训练多个不同分量数的模型，计算开销大；而贝叶斯 GMM 通过变分推断和 Dirichlet 过程先验，在一次训练过程中就能自动驱动无效组件的权重趋近于零，从而“激活”出真正需要的组件数，既节省了计算又提供了不确定性估计，但代价是模型假设更复杂（需接受变分近似）和收敛对初始化更敏感。
-
为什么 EllipticEnvelope 默认使用 MCD 而不是经验协方差？
因为经验协方差对离群点极其敏感：即使少量离群点也会显著扭曲协方差估计，导致马氏距离失真，从而使得基于经验协方差的椭圆包络既可能误将正常点标记为异常（过敏），又可能遗漏真实异常（过保守）。MCD 拥有高达 50% 的破坏点，能够在存在大量离群点时仍能可靠地估计内生数据的协方差结构，是构建稳健异常检测器的理想基础。
-
标准化 vs 归一化的选择策略是什么？
StandardScaler（标准化）假设数据近似高斯分布，对离群点敏感；RobustScaler（鲁棒缩放）使用中位数和四分位距，对离群点不敏感；MinMaxScaler（最小-最大缩放）将数据压缩到固定区间，保留稀疏性但受离群点影响大；MaxAbsScaler（最大绝对值缩放）仅缩放不平移，保留稀疏矩阵结构。选择时需考虑数据分布形态、离群点存在性以及下游模型对尺度的敏感度。
-
OneHotEncoder 与 OrdinalEncoder 的适用边界在哪里？
OneHotEncoder 适合名义类别（无序），生成稀疏矩阵避免引入虚假序数关系，但高基数特征会导致维度爆炸；OrdinalEncoder 适合有序类别或树模型（能处理整数编码的序数关系），输出单列整数，内存高效。对于高基数名义特征，可考虑 TargetEncoder 或 HashingVectorizer 等替代方案。
-
KBinsDiscretizer 的三种策略如何权衡？
uniform（等宽）简单快速但对离群点敏感；quantile（等频）使每箱样本数均衡，鲁棒于分布偏态；kmeans（聚类）基于数据密度分箱，能捕捉自然簇结构但计算开销大。选择时需考虑数据分布特征、离群点影响及计算预算。
-
FunctionTransformer 的 validate 参数为何默认 False？
默认 validate=False 允许任意可调用对象处理任意数据结构（如 DataFrame、列表），提供最大灵活性；设为 True 时强制转为二维数组，适合标准数值管道但限制了输入类型。这是权衡通用性与类型安全的结果。
-
ColumnTransformer 如何平衡异构列的并行处理与输出一致性？
并行拟合各列变换器（通过 n_jobs），再水平堆叠结果；通过 sparse_threshold 控制稀疏/稠密输出策略；通过 verbose_feature_names_out 自动添加变换器前缀避免特�名冲突。设计核心是：输入列互斥分区、变换器独立、输出拼接。
-
TransformedTargetRegressor 为何分离目标变换与回归器？
解耦目标空间变换（如 log、Box-Cox）与回归建模，使回归器始终在变换后的高斯近似空间工作，预测时自动逆变换。这种设计避免了手动管道拼接的错误，并支持元数据路由传递 sample_weight 等参数。
-
文本向量化器三种策略的核心权衡是什么？
CountVectorizer（精确词频）内存随词汇表增长，适合中小规模；TfidfVectorizer（TF-IDF）引入 IDF 抑制高频词，提升区分度但需两遍扫描；HashingVectorizer（哈希技巧）固定维度、无状态、可流式处理，但不可逆、有碰撞风险。选择取决于语料规模、内存预算、是否需逆向映射。
-
缺失值插补策略的层级递进逻辑是什么？
SimpleImputer（单变量统计量）最快但忽略特征间关系；IterativeImputer（多变量回归）建模特征依赖但计算重、需指定回归器；KNNImputer（邻域平均）非参数、利用局部结构但高维下距离失效、预测慢。实践中常先用 SimpleImputer 基线，再尝试迭代或 KNN 提升。
27.8 动手练习
-
阅读 GMM 的 EM 核心实现
-
文件：
sklearn/mixture/_gaussian_mixture.py -
方法：
_e_step、_m_step、_estimate_gaussian_parameters -
问题：E 步如何利用 log-sum-exp 技巧避免数值下溢？M 步中不同 covariance_type 的协方差更新公式有何区别？
-
-
分析贝叶斯 GMM 的变分推断与稳健协方差估计
-
文件：
sklearn/mixture/_bayesian_mixture.py、sklearn/covariance/_robust_covariance.py -
方法：
_compute_lower_bound、_estimate_weights、fast_mcd、_c_step -
问题：变分下界 ELBO 的各项物理意义是什么？FastMCD 的 C-步为何能单调减小行列式？
-
-
探索图拉索与椭圆包络异常检测
-
文件：
sklearn/covariance/_graph_lasso.py、sklearn/covariance/_elliptic_envelope.py -
方法：
_graphical_lasso（坐标下降内循环）、decision_function -
问题：坐标下降如何逐列更新精度矩阵？EllipticEnvelope 的 offset_ 如何由 contamination 参数决定？
-
-
实践数据预处理：标准化与归一化对比
-
文件：
sklearn/preprocessing/_data.py -
类：
StandardScaler、MinMaxScaler、RobustScaler、MaxAbsScaler -
问题：对含离群点的数据分别应用四种缩放器，观察变换后分布差异；RobustScaler 为何用 IQR 而非标准差？
-
-
深入类别编码器：OneHotEncoder 与 TargetEncoder
-
文件：
sklearn/preprocessing/_encoders.py、sklearn/preprocessing/_target_encoder.py -
方法：
OneHotEncoder.fit_transform、TargetEncoder.fit_transform（交叉拟合逻辑） -
问题：OneHotEncoder 如何处理 handle_unknown='infrequent_if_exist'？TargetEncoder 的交叉拟合如何防止目标泄露？
-
-
文本向量化器对比：CountVectorizer vs TfidfVectorizer vs HashingVectorizer
-
文件：
sklearn/feature_extraction/text.py -
方法：
_count_vocab、fit_transform（TF-IDF 两遍扫描）、_get_hasher -
问题：在固定语料上对比三者输出维度、内存占用、可逆性；HashingVectorizer 的 alternate_sign 如何近似保持内积？
-
-
缺失值插补策略对比：SimpleImputer vs IterativeImputer vs KNNImputer
-
文件：
sklearn/impute/_base.py、sklearn/impute/_iterative.py、sklearn/impute/_knn.py -
方法：
SimpleImputer._dense_fit、IterativeImputer._impute_one_feature、KNNImputer._calc_impute -
问题：构造含不同缺失模式（MCAR/MAR/MNAR）的数据，对比三种插补器的 RMSE；IterativeImputer 为何需要初始插补？
-
-
组合工具实战：ColumnTransformer 与 TransformedTargetRegressor
-
文件：
sklearn/compose/_column_transformer.py、sklearn/compose/_target.py -
方法：
ColumnTransformer._call_func_on_transformers、TransformedTargetRegressor._fit_transformer -
问题：ColumnTransformer 如何并行执行多个变换器？TransformedTargetRegressor 如何确保 predict 时自动应用 inverse_transform？
-
27.9 本章小结
这一章中我们学习了高斯混合模型及其贝叶斯变体如何通过 EM 算法和变分推断来聚类和密度估计，理解了不同协方差约束如何适应数据形状；深入探讨了稳健协方差估计方法（如 Ledoit-Wolf 收缩、图拉索、MCD）如何应对离群点和小样本问题，以及它们如何驱动异常检测器 EllipticEnvelope；系统梳理了特征缩放、编码、离散化、分布变换等核心预处理器的实现细节，掌握了类别特征如何通过 OneHotEncoder、OrdinalEncoder 等编码器转换为机器可处理的形式；学习了如何使用 PolynomialFeatures 和 SplineTransformer 构造非线性特征，以及如何用 KBinsDiscretizer 和 FunctionTransformer 进行离散化和自定义变换；掌握了 ColumnTransformer 如何实现异构列的并行处理，以及 TransformedTargetRegressor 如何通过目标变换来增强回归模型的灵活性；最后，我们理解了文本特征如何通过 CountVectorizer、TfidfVectorizer 和 HashingVectorizer 从原始文本转化为数值向量，为后续的文本挖掘奠定基础。
下一章中，我们将学习模型解释与高斯过程 —— 打开黑盒的“X 光机”。
第 28 章 —— 模型解释与高斯过程 —— 打开黑盒的“X 光机”
28.1 学习目标
-
难度：★★★☆☆（3/5）
-
预备知识：Python 基础、面向对象编程与 Markdown/代码阅读基础
-
理解部分依赖图（PDP）与个体条件期望（ICE）的原理及两种计算方法（递归法与暴力法）的区别
-
掌握排列重要性的计算流程与并行化实现机制
-
理解决策边界可视化的核心步骤与多分类场景下的颜色映射策略
-
理解高斯过程内核的设计模式（组合、运算符重载、超参数管理）及常用内核（RBF、Matern、RationalQuadratic）的数学形式
-
理解高斯过程回归（GPR）的后验推断、对数边际似然优化与不确定性量化预测
-
理解高斯过程分类（GPC）中拉普拉斯近似的原理、二分类后验模式求解及多分类扩展策略
28.2 生活类比
想象模型解释工具是一台“透视仪”，帮助我们看穿黑盒模型的内部结构。部分依赖图（PDP）就像一张“平均透视切片”：固定某个特征，平均其他特征的影响，揭示该特征对模型预测的全局趋势，就像医学CT扫描的平均切片图。而个体条件期望（ICE）则是“个体透视切片”：为每条样本单独绘制一条特征变化下的预测曲线，揭示不同样本对同一特征的个体差异，就像为每位病人单独建模。在检验特征重要性时，排列重要性采用“打乱测试法”：像考试抽查，随机打乱某个特征的取值（相当于抽走答案卡），观察模型总分下降多少，下降越大说明该特征越重要，且支持“抽样考生”（子采样）来加速计算。决策边界可视化则是“地形图绘制”：在二维特征平面上铺设密集网格，对每个网格点预测类别或概率，连接相同预测的点形成等高线，就像绘制山脉的等高线图。多分类问题下，不同类别用不同颜色填充，形成彩色地形图。高斯过程内核就像“相似度乐高积木”：RBF、Matern等基础积木定义了任意两点之间的协方差（相似度），通过加法（+）、�法（*）、幂运算（**）组合搭建更复杂的协方差结构，超参数是积木的“规格参数”（如长度尺），支持梯度优化自动调节。高斯过程回归（GPR）是“带误差棒的插值器”：不仅给出预测均值（最佳猜测），还提供方差（不确定性），通过最大化对数边际似然自动调节核的超参数。高斯过程分类（GPC）中的拉普拉斯近似则是“非高斯后验的高斯拟合”：分类似然本非高斯分布，用高斯分布（二阶泰勒展开）近似它，先用牛顿法迭代寻找后验模式（峰值），再用误差函数（erf）近似积分得到概率预测。
28.3 源码地图
sklearn/inspection/_partial_dependence.py
├── _grid_from_X() # 生成特征网格（支持分位数、分类、自定义值）
├── _partial_dependence_recursion() # 树模型专用递归法计算PDP
├── _partial_dependence_brute() # 通用暴力法计算PDP/ICE
├── partial_dependence() # 主入口：参数校验、方法选择、结果聚合
sklearn/inspection/_pd_utils.py
├── _check_feature_names() # 特征名校验与默认生成
├── _get_feature_index() # 特征名转索引
sklearn/inspection/_plot/partial_dependence.py
├── PartialDependenceDisplay # PDP/ICE可视化类
│ ├── init() # 初始化存储计算结果与绘图参数
│ ├── from_estimator() # 类方法：计算并绘制（并行调度partial_dependence）
│ ├── plot() # 绘制入口：网格布局、1D/2D分发
│ ├── _get_sample_count() # 计算ICE子采样数量
│ ├── _plot_one_way_partial_dependence() # 1D曲线/柱状图绘制
│ ├── _plot_two_way_partial_dependence() # 2D等高线/热力图绘制
│ ├── _plot_ice_lines() # ICE曲线绘制（支持采样）
│ └── _plot_average_dependence() # 平均PD曲线/柱状图绘制
sklearn/inspection/_permutation_importance.py
├── _weights_scorer() # 带样本权重的评分包装
├── _calculate_permutation_scores() # 单特征置换计算核心（支持子采样）
├── _create_importances_bunch() # 重要性统计量聚合（均值/标准差）
├── permutation_importance() # 主入口：并行调度、基线评分、结果聚合
sklearn/inspection/_plot/decision_boundary.py
├── _check_boundary_response_method() # 响应方法校验与推断
├── DecisionBoundaryDisplay # 决策边界可视化类
│ ├── init() # 存储网格、响应、类别数等
│ ├── plot() # 绘制入口：contourf/contour/pcolormesh
│ └── from_estimator() # 类方法：网格生成、预测、响应后处理
sklearn/gaussian_process/kernels.py
├── Hyperparameter # 超参数规范（namedtuple子类）
├── Kernel (ABC) # 基类：参数管理、theta/bounds属性、运算符重载
│ ├── get_params/set_params/clone_with_theta
│ ├── add/mul/pow # 运算符重载 -> Sum/Product/Exponentiation
├── CompoundKernel # 多核堆叠（用于多输出/多分类）
├── KernelOperator # 二元运算符基类
│ ├── Sum # k1 + k2
│ ├── Product # k1 * k2
│ └── Exponentiation # k ** p
├── StationaryKernelMixin/NormalizedKernelMixin/GenericKernelMixin
├── ConstantKernel/WhiteKernel # 基础核
├── RBF # 径向基核（含各向异性、梯度计算）
├── Matern(RBF) # Matern核（nu=0.5/1.5/2.5/inf解析式）
├── RationalQuadratic # 有理二次核
├── ExpSineSquared # 周期核
├── DotProduct # 非平稳核
└── PairwiseKernel # sklearn.metrics.pairwise封装
sklearn/gaussian_process/_gpr.py
├── GaussianProcessRegressor # GPR主类
│ ├── init() # 初始化核、优化器、正则化参数等
│ ├── fit() # 训练：超参数优化、Cholesky分解、alpha计算
│ ├── predict() # 预测：后验均值/方差/协方差
│ ├── sample_y() # 后验采样
│ ├── log_marginal_likelihood() # 边际似然及梯度（含多输出）
│ ├── _constrained_optimization() # L-BFGS-B优化封装
│ └── sklearn_tags() # 标记无需fit即可预测
sklearn/gaussian_process/_gpc.py
├── _BinaryGaussianProcessClassifierLaplace # 二分类核心实现
│ ├── init() # 初始化核、优化器、预测迭代上限等
│ ├── fit() # 训练：超参数优化、后验模式求解
│ ├── predict/predict_proba() # 预测：MAP决策 / 误差函数近似积分
│ ├── log_marginal_likelihood() # 边际似然及梯度（Algorithm 5.1）
│ ├── latent_mean_and_variance() # 潜变量均值/方差（Algorithm 3.2）
│ ├── _posterior_mode() # 牛顿法求后验模式（Algorithm 3.1）
│ └── _constrained_optimization() # 优化封装
└── GaussianProcessClassifier # 多分类包装器
├── init() # 初始化多类策略、并行数等
├── fit() # OvR/OvO策略分发
├── predict/predict_proba() # 委托基估计器
├── kernel_ (property) # CompoundKernel聚合
├── log_marginal_likelihood() # 多类别平均似然
└── latent_mean_and_variance() # 仅支持二分类
28.4 部分依赖与 ICE —— 特征影响的“沙盘推演”
部分依赖图（PDP）与个体条件期望（ICE）是理解“黑盒”模型如何响应特征变化的关键工具。PDP 展示特征对模型预测的平均影响，揭示全局趋势；ICE 展示每个样本对特征变化的个体响应，揭示异质性。两者结合提供模型行为的全面解读，帮助理解特征作用机制。PDP 可通过递归（针对树模型）或暴力（通用）方法高效计算。
我们从源码入手，先看 _grid_from_X 函数如何生成计算网格。
源码路径：sklearn/inspection/_partial_dependence.py - _grid_from_X()（40-110行）
def _grid_from_X(X, percentiles, is_categorical, grid_resolution, custom_values):
"""Generate a grid of points based on the percentiles of X."""
if not isinstance(percentiles, Iterable) or len(percentiles) != 2:
raise ValueError("'percentiles' must be a sequence of 2 elements.")
if not all(0 <= x <= 1 for x in percentiles):
raise ValueError("'percentiles' values must be in [0, 1].")
if percentiles[0] >= percentiles[1]:
raise ValueError("percentiles[0] must be strictly less than percentiles[1].")
if grid_resolution <= 1:
raise ValueError("'grid_resolution' must be strictly greater than 1.")
def _convert_custom_values(values):
# Convert custom types such that object types are always used for string arrays
dtype = object if any(isinstance(v, str) for v in values) else None
return np.asarray(values, dtype=dtype)
custom_values = {k: _convert_custom_values(v) for k, v in custom_values.items()}
if any(v.ndim != 1 for v in custom_values.values()):
error_string = ", ".join(
f"Feature {k}: {v.ndim} dimensions"
for k, v in custom_values.items()
if v.ndim != 1
)
raise ValueError(
"The custom grid for some features is not a one-dimensional array. "
f"{error_string}"
)
values = []
for feature, is_cat in enumerate(is_categorical):
if feature in custom_values:
axis = custom_values[feature]
else:
try:
uniques = np.unique(_safe_indexing(X, feature, axis=1))
except TypeError as exc:
raise ValueError(
f"The column #{feature} contains mixed data types. Finding unique "
"categories fail due to sorting. It usually means that the column "
"contains `np.nan` values together with `str` categories. Such use "
"case is not yet supported in scikit-learn."
) from exc
if is_cat or uniques.shape[0] < grid_resolution:
axis = uniques
else:
emp_percentiles = mquantiles(
_safe_indexing(X, feature, axis=1), prob=percentiles, axis=0
)
if np.allclose(emp_percentiles[0], emp_percentiles[1]):
raise ValueError(
"percentiles are too close to each other, "
"unable to build the grid. Please choose percentiles "
"that are further apart."
)
axis = np.linspace(
emp_percentiles[0],
emp_percentiles[1],
num=grid_resolution,
endpoint=True,
)
values.append(axis)
return cartesian(values), values
这段代码定义了特征网格的生成逻辑：根据输入数据的分位数、唯一值或自定义值构建每个特征的取值轴，然后通过笛卡尔积生成多维网格。它支持三种网格生成方式：分位数法（默认）、唯一值法（用于类别特征或低基数特征）和自定义值法。
接下来，我们看树模型专用的递归法 _partial_dependence_recursion。
源码路径：sklearn/inspection/_partial_dependence.py - _partial_dependence_recursion()（112-138行）
def _partial_dependence_recursion(est, grid, features):
"""Calculate partial dependence via the recursion method."""
averaged_predictions = est._compute_partial_dependence_recursion(grid, features)
if averaged_predictions.ndim == 1:
averaged_predictions = averaged_predictions.reshape(1, -1)
return averaged_predictions
这个函数实际上是一个包装器，真正的计算委托给了估计器自身的 _compute_partial_dependence_recursion 方法。递归法的核心思想是：对于树模型，遍历树的每个节点，如果节点分裂特征是我们感兴趣的特征，则只进入对应分支；否则按训练样本在该节点的比例加权进入左右两个分支。最后用所有访问到的叶子节点的预测值求加权平均。这种方法避免了对每个网格点重新预测全部样本，因此对树模型非常高效。
然而，递归法有局限：它不支持个体条件期望（ICE），因为它内在地计算了所有样本的平均值（通过加权遍历隐式求了ICE的平均）。要得到ICE，需要暴力法。
我们现在看通用暴力法 _partial_dependence_brute。
源码路径：sklearn/inspection/_partial_dependence.py - _partial_dependence_brute()（140-210行）
def _partial_dependence_brute(
est, grid, features, X, response_method, sample_weight=None
):
"""Calculate partial dependence via the brute force method."""
predictions = []
averaged_predictions = []
if response_method == "auto":
response_method = (
"predict" if is_regressor(est) else ["predict_proba", "decision_function"]
)
X_eval = X.copy()
for new_values in grid:
for i, variable in enumerate(features):
_safe_assign(X_eval, new_values[i], column_indexer=variable)
pred, _ = _get_response_values(est, X_eval, response_method=response_method)
predictions.append(pred)
averaged_predictions.append(np.average(pred, axis=0, weights=sample_weight))
n_samples = X.shape[0]
predictions = np.array(predictions).T
if is_regressor(est) and predictions.ndim == 2:
predictions = predictions.reshape(n_samples, -1)
elif is_classifier(est) and predictions.shape[0] == 2:
predictions = predictions[1]
predictions = predictions.reshape(n_samples, -1)
averaged_predictions = np.array(averaged_predictions).T
if averaged_predictions.ndim == 1:
averaged_predictions = averaged_predictions.reshape(1, -1)
return averaged_predictions, predictions
这段代码实现了暴力法：为网格中的每个点，复制原始数据 X，将目标特征列替换为网格值，然后对修改后的数据进行预测，最后对所有样本的预测求平均（支持样本权重）。暴力法的优点是通用——任何实现 predict、predict_proba 或 decision_function 的估计器都可以使用；缺点是计算开销大，因为对于网格中的每个点，都需要对全部样本重新预测一次。它支持 ICE（返回所有样本的预测矩阵）和样本权重。
主入口函数 partial_dependence 负责统一调度。
源码路径：sklearn/inspection/_partial_dependence.py - partial_dependence()（212-350行）
@validate_params(
{
"estimator": [
HasMethods(["fit", "predict"]),
HasMethods(["fit", "predict_proba"]),
HasMethods(["fit", "decision_function"]),
],
"X": ["array-like", "sparse matrix"],
"features": ["array-like", Integral, str],
"sample_weight": ["array-like", None],
"categorical_features": ["array-like", None],
"feature_names": ["array-like", None],
"response_method": [StrOptions({"auto", "predict_proba", "decision_function"})],
"percentiles": [tuple],
"grid_resolution": [Interval(Integral, 1, None, closed="left")],
"method": [StrOptions({"auto", "recursion", "brute"})],
"kind": [StrOptions({"average", "individual", "both"})],
"custom_values": [dict, None],
},
prefer_skip_nested_validation=True,
)
def partial_dependence(
estimator,
X,
features,
*,
sample_weight=None,
categorical_features=None,
feature_names=None,
response_method="auto",
percentiles=(0.05, 0.95),
grid_resolution=100,
custom_values=None,
method="auto",
kind="average",
):
"""Partial dependence of ``features``."""
check_is_fitted(estimator)
if not (is_classifier(estimator) or is_regressor(estimator)):
raise ValueError("'estimator' must be a fitted regressor or classifier.")
if is_classifier(estimator) and isinstance(estimator.classes_[0], np.ndarray):
raise ValueError("Multiclass-multioutput estimators are not supported")
if not (hasattr(X, "__array__") or sparse.issparse(X)):
X = check_array(X, ensure_all_finite="allow-nan", dtype=object)
if is_regressor(estimator) and response_method != "auto":
raise ValueError(
"The response_method parameter is ignored for regressors and "
"must be 'auto'."
)
if kind != "average":
if method == "recursion":
raise ValueError(
"The 'recursion' method only applies when 'kind' is set to 'average'"
)
method = "brute"
if method == "recursion" and sample_weight is not None:
raise ValueError(
"The 'recursion' method can only be applied when sample_weight is None."
)
if method == "auto":
if sample_weight is not None:
method = "brute"
elif isinstance(estimator, BaseGradientBoosting) and estimator.init is None:
method = "recursion"
elif isinstance(
estimator,
(BaseHistGradientBoosting, DecisionTreeRegressor, RandomForestRegressor),
):
method = "recursion"
else:
method = "brute"
if method == "recursion":
if not isinstance(
estimator,
(
BaseGradientBoosting,
BaseHistGradientBoosting,
DecisionTreeRegressor,
RandomForestRegressor,
),
):
supported_classes_recursion = (
"GradientBoostingClassifier",
"GradientBoostingRegressor",
"HistGradientBoostingClassifier",
"HistGradientBoostingRegressor",
"HistGradientBoostingRegressor",
"DecisionTreeRegressor",
"RandomForestRegressor",
)
raise ValueError(
"Only the following estimators support the 'recursion' "
"method: {}. Try using method='brute'.".format(
", ".join(supported_classes_recursion)
)
)
if response_method == "auto":
response_method = "decision_function"
if response_method != "decision_function":
raise ValueError(
"With the 'recursion' method, the response_method must be "
"'decision_function'. Got {}.".format(response_method)
)
if sample_weight is not None:
sample_weight = _check_sample_weight(sample_weight, X)
if _determine_key_type(features, accept_slice=False) == "int":
if np.any(np.less(features, 0)):
raise ValueError("all features must be in [0, {}]".format(X.shape[1] - 1))
features_indices = np.asarray(
_get_column_indices(X, features), dtype=np.intp, order="C"
).ravel()
feature_names = _check_feature_names(X, feature_names)
n_features = X.shape[1]
if categorical_features is None:
is_categorical = [False] * len(features_indices)
else:
categorical_features = np.asarray(categorical_features)
if categorical_features.size == 0:
raise ValueError(
"Passing an empty list (`[]`) to `categorical_features` is not "
"supported. Use `None` instead to indicate that there are no "
"categorical features."
)
if categorical_features.dtype.kind == "b":
is_categorical = [categorical_features[idx] for idx in features_indices]
elif categorical_features.dtype.kind in ("i", "O", "U"):
categorical_features_idx = [
_get_feature_index(cat, feature_names=feature_names)
for cat in categorical_features
]
is_categorical = [
idx in categorical_features_idx for idx in features_indices
]
else:
raise ValueError(
"Expected `categorical_features` to be an array-like of boolean,"
f" integer, or string. Got {categorical_features.dtype} instead."
)
custom_values = custom_values or {}
if isinstance(features, (str, int)):
features = [features]
for feature_idx, feature, is_cat in zip(features_indices, features, is_categorical):
if is_cat:
continue
if _safe_indexing(X, feature_idx, axis=1).dtype.kind in "iu":
warnings.warn(
f"The column {feature!r} contains integer data. Partial "
"dependence plots are not supported for integer data: this "
"can lead to implicit rounding with NumPy arrays or even errors "
"with newer pandas versions. Please convert numerical features"
"to floating point dtypes ahead of time to avoid problems. "
"This will raise ValueError in scikit-learn 1.9.",
FutureWarning,
)
break
X_subset = _safe_indexing(X, features_indices, axis=1)
custom_values_for_X_subset = {
index: custom_values.get(feature)
for index, feature in enumerate(features)
if feature in custom_values
}
grid, values = _grid_from_X(
X_subset,
percentiles,
is_categorical,
grid_resolution,
custom_values_for_X_subset,
)
if method == "brute":
averaged_predictions, predictions = _partial_dependence_brute(
estimator, grid, features_indices, X, response_method, sample_weight
)
predictions = predictions.reshape(
-1, X.shape[0], *[val.shape[0] for val in values]
)
else:
averaged_predictions = _partial_dependence_recursion(
estimator, grid, features_indices
)
averaged_predictions = averaged_predictions.reshape(
-1, *[val.shape[0] for val in values]
)
pdp_results = Bunch(grid_values=values)
if kind == "average":
pdp_results["average"] = averaged_predictions
elif kind == "individual":
pdp_results["individual"] = predictions
else: # kind='both'
pdp_results["average"] = averaged_predictions
pdp_results["individual"] = predictions
return pdp_results
这个函数是部分依赖计算的总控台：它先做参数校验和预处理（如特征名处理、分类特征识别），然后根据 method 和 kind 参数选择计算路径（递归法或暴力法），最后将结果包装成 Bunch 对象返回。关键设计包括：递归法仅用于特定树模型且要求 kind='average' 和 sample_weight=None；当需要 ICE 或样本权重时，自动降级到暴力法。
为了更直观地理解这两种方法的区别，我们可以看一个简单的数据流图。
不同于暴力法对每个网格点重新预测全部样本，递归法通过一次遍历树结构就能得到所有网格点的平均预测，这正是它高效的原因。
下面我们使用 Mermaid 展示完整的数据流，从网格生成到最终结果。
有了计算基础，我们现在看如何将这些结果可视化。PartialDependenceDisplay 类负责将计算结果转化为图表。
源码路径：sklearn/inspection/_plot/partial_dependence.py - PartialDependenceDisplay.__init__()（60-120行）
def __init__(
self,
pd_results,
*,
features,
feature_names,
target_idx,
deciles,
kind="average",
subsample=1000,
random_state=None,
is_categorical=None,
):
self.pd_results = pd_results
self.features = features
self.feature_names = feature_names
self.target_idx = target_idx
self.deciles = deciles
self.kind = kind
self.subsample = subsample
self.random_state = random_state
self.is_categorical = is_categorical
这个构造函数简单地存储了所有必要的元数据：计算结果、features、feature_names、target_idx（用于多分类/多输出）、分位数用于绘图参考线、kind 决定绘制平均线还是ICE线，以及用于ICE子采样的 subsample 和 random_state。
真正的绘图发生在 plot 方法中，它会根据特征数量（1D或2D）、特征类型（连续或分类）和 kind 参数分发到不同的绘图函数。
源码路径：sklearn/inspection/_plot/partial_dependence.py - PartialDependenceDisplay.plot()（720-900行，简化关键逻辑）
def plot(self, *, ax=None, n_cols=3, ...):
# ... 参数预处理和中心化处理 ...
if not centered:
pd_results_ = self.pd_results
else:
# 中心化处理：使ICE和PD线在y=0处对齐
pd_results_ = []
for kind_plot, pd_result in zip(kind, self.pd_results):
current_results = {"grid_values": pd_result["grid_values"]}
if kind_plot in ("individual", "both"):
preds = pd_result.individual
preds = preds - preds[self.target_idx, :, 0, None]
current_results["individual"] = preds
if kind_plot in ("average", "both"):
avg_preds = pd_result.average
avg_preds = avg_preds - avg_preds[self.target_idx, 0, None]
current_results["average"] = avg_preds
pd_results_.append(Bunch(**current_results))
# ... 处理pdp_lim（全局y轴范围） ...
# 设置图表布局
if isinstance(ax, plt.Axes):
# 单个axes情况：创建网格布局
# ... 创建子图网格 ...
else:
# 数组axes情况：直接使用传入的axes
# ... 验证axes数量 ...
# 遍历每个特征进行绘图
for pd_plot_idx, (axi, feature_idx, cat, pd_result, kind_plot) in enumerate(
zip(self.axes_.ravel(), self.features, is_categorical, pd_results_, kind)
):
# 提取数据
avg_preds = None
preds = None
feature_values = pd_result["grid_values"]
if kind_plot == "individual":
preds = pd_result.individual
elif kind_plot == "average":
avg_preds = pd_result.average
else: # kind_plot == 'both'
avg_preds = pd_result.average
preds = pd_result.individual
if len(feature_values) == 1:
# 1D情况：绘制曲线或柱状图
# ... 设置线型/柱状图默认参数 ...
# ICE线条绘制（如果需要）
if kind_plot in ("individual", "both"):
self._plot_ice_lines(preds[self.target_idx], feature_values, ...)
# 平均依赖线/柱状图绘制（如果需要）
if kind_plot in ("average", "both"):
self._plot_average_dependence(avg_preds[self.target_idx].ravel(), feature_values, axi, ...)
# ... 设置轴标签、图例等 ...
else:
# 2D情况：绘制等高线图或热力图
self._plot_two_way_partial_dependence(avg_preds, feature_values, feature_idx, axi, ...)
return self
1D情况下，_plot_one_way_partial_dependence 负责绘制；2D情况下，_plot_two_way_partial_dependence 负责绘制等高线图或热力图。
对于1D图（连续特征），它可以绘制ICE线条（半透明、采样）、平均PD线（实线）、以及分位数标记；对于分类特征，它绘制柱状图。
对于2D图，如果两个特征都是连续的，它绘制填充等高线图（contourf）；如果任意一个是分类的，它绘制热力图（imshow），并在每个单元格中标注数值。
ICE线条的绘制会进行子采样以避免图形过于密集。
源码路径：sklearn/inspection/_plot/partial_dependence.py - PartialDependenceDisplay._plot_ice_lines()（510-540行）
def _plot_ice_lines(
self,
preds,
feature_values,
n_ice_to_plot,
ax,
pd_plot_idx,
n_total_lines_by_plot,
individual_line_kw,
):
"""Plot the ICE lines."""
rng = check_random_state(self.random_state)
ice_lines_idx = rng.choice(
preds.shape[0],
n_ice_to_plot,
replace=False,
)
ice_lines_subsampled = preds[ice_lines_idx, :]
for ice_idx, ice in enumerate(ice_lines_subsampled):
line_idx = np.unravel_index(
pd_plot_idx * n_total_lines_by_plot + ice_idx, self.lines_.shape
)
self.lines_[line_idx] = ax.plot(
feature_values, ice.ravel(), **individual_line_kw
)[0]
这个函数从所有样本中随机选择 n_ice_to_plot 条ICE曲线进行绘制，使用透明度和细线宽以避免遮挡平均线。
平均依赖线的绘制则更为直接。
源码路径：sklearn/inspection/_plot/partial_dependence.py - PartialDependenceDisplay._plot_average_dependence()（545-580行）
def _plot_average_dependence(
self,
avg_preds,
feature_values,
ax,
pd_line_idx,
line_kw,
categorical,
bar_kw,
):
"""Plot the average partial dependence."""
if categorical:
bar_idx = np.unravel_index(pd_line_idx, self.bars_.shape)
self.bars_[bar_idx] = ax.bar(feature_values, avg_preds, **bar_kw)[0]
ax.tick_params(axis="x", rotation=90)
else:
line_idx = np.unravel_index(pd_line_idx, self.lines_.shape)
self.lines_[line_idx] = ax.plot(
feature_values,
avg_preds,
**line_kw,
)[0]
对于连续特征，它绘制一条线；对于分类特征，它绘制柱状图（并旋转x轴标签以避免重叠）。
理解了PDP和ICE，我们现在转向另一个重要的模型解释技术：排列重要性。
28.5 排列重要性与决策边界 —— 特征重要性的“置换检验”与“边界可视化”
排列重要性通过随机打乱特征列来衡量其对模型性能的贡献：基线得分减去打乱后得分得到原始重要性分数。它支持多次重复以减少方差，并可通过 max_samples 控制计算开销。决策边界可视化则通过在特征空间上构建密集网格并预测每个点的类别来绘制决策边界。
我们先看排列重要性的主函数 permutation_importance。
源码路径：sklearn/inspection/_permutation_importance.py - permutation_importance()（220-290行）
@validate_params(
{
"estimator": [HasMethods(["fit"])],
"X": ["array-like"],
"y": ["array-like", None],
"scoring": [
StrOptions(set(get_scorer_names())),
callable,
list,
tuple,
dict,
None,
],
"n_repeats": [Interval(Integral, 1, None, closed="left")],
"n_jobs": [Integral, None],
"random_state": ["random_state"],
"sample_weight": ["array-like", None],
"max_samples": [
Interval(Integral, 1, None, closed="left"),
Interval(RealNotInt, 0, 1, closed="right"),
],
},
prefer_skip_nested_validation=True,
)
def permutation_importance(
estimator,
X,
y,
*,
scoring=None,
n_repeats=5,
n_jobs=None,
random_state=None,
sample_weight=None,
max_samples=1.0,
):
"""Permutation importance for feature evaluation [BRE]_."""
if not hasattr(X, "iloc"):
X = check_array(X, ensure_all_finite="allow-nan", dtype=None)
random_state = check_random_state(random_state)
random_seed = random_state.randint(np.iinfo(np.int32).max + 1)
if not isinstance(max_samples, numbers.Integral):
max_samples = int(max_samples * X.shape[0])
elif max_samples > X.shape[0]:
raise ValueError("max_samples must be <= n_samples")
scorer = check_scoring(estimator, scoring=scoring)
baseline_score = _weights_scorer(scorer, estimator, X, y, sample_weight)
scores = Parallel(n_jobs=n_jobs)(
delayed(_calculate_permutation_scores)(
estimator,
X,
y,
sample_weight,
col_idx,
random_seed,
n_repeats,
scorer,
max_samples,
)
for col_idx in range(X.shape[1])
)
if isinstance(baseline_score, dict):
return {
name: _create_importances_bunch(
baseline_score[name],
np.array([scores[col_idx][name] for col_idx in range(X.shape[1])]),
)
for name in baseline_score
}
else:
return _create_importances_bunch(baseline_score, np.array(scores))
这个函数的总体流程是：
-
预处理输入数据
X（转换为NumPy数组，允许NaN） -
处理
max_samples参数：如果是比例则转换为绝对数量，并检查是否超过样本数 -
创建评分器
scorer（处理字符串、可调用对象或字典等多种输入形式） -
计算基线得分
baseline_score（在原始数据上评估） -
并行计算每个特征的置换得分：对于每个特征列，复制数据，打乱该列，重新评估，重复
n_repeats次 -
聚合结果：如果基线得分是字典（多指标情况），则为每个指标创建一个
Bunch；否则直接创建一个Bunch
并行化的粒度是以特征列为单位：外层循环遍历每个特征索引 col_idx，内层的 _calculate_permutation_scores 函数处理单个特征的所有重复。
为了保证多进程/多线程下的可复现性，函数在并行开始前从主随机状态中生成一个整数种子 random_seed，并将其传递给每个并行作业。每个作业内部再次使用 check_random_state 确保拥有独立的随机状态实例，但所有作业共享相同的初始种子，从而保证在相同输入下结果可复现。
我们现在看单特征置换的核心函数 _calculate_permutation_scores。
源码路径：sklearn/inspection/_permutation_importance.py - _calculate_permutation_scores()（130-170行）
def _calculate_permutation_scores(
estimator,
X,
y,
sample_weight,
col_idx,
random_state,
n_repeats,
scorer,
max_samples,
):
"""Calculate score when `col_idx` is permuted."""
random_state = check_random_state(random_state)
if max_samples < X.shape[0]:
row_indices = _generate_indices(
random_state=random_state,
bootstrap=False,
n_population=X.shape[0],
n_samples=max_samples,
)
X_permuted = _safe_indexing(X, row_indices, axis=0)
y = _safe_indexing(y, row_indices, axis=0)
if sample_weight is not None:
sample_weight = _safe_indexing(sample_weight, row_indices, axis=0)
else:
X_permuted = X.copy()
scores = []
shuffling_idx = np.arange(X_permuted.shape[0])
for _ in range(n_repeats):
random_state.shuffle(shuffling_idx)
if hasattr(X_permuted, "iloc"):
col = X_permuted.iloc[shuffling_idx, col_idx]
col.index = X_permuted.index
X_permuted[X_permuted.columns[col_idx]] = col
else:
X_permuted[:, col_idx] = X_permuted[shuffling_idx, col_idx]
scores.append(_weights_scorer(scorer, estimator, X_permuted, y, sample_weight))
if isinstance(scores[0], dict):
scores = _aggregate_score_dicts(scores)
else:
scores = np.array(scores)
return scores
这个函数实现了单个特征的置换检验：
-
如果
max_samples < X.shape[0]，则进行无放回子采样：使用_generate_indices生成行索引，然后用_safe_indexing从X、y和sample_weight中采样对应的行 -
否则，直接复制整个
X作为X_permuted -
初始化一个索引数组
shuffling_idx = np.arange(X_permuted.shape[0]) -
重复
n_repeats次：a. 随机打乱
shuffling_idxb. 如果
X_permuted是DataFrame（有iloc属性），则使用.iloc进行列赋值以保持索引对齐；否则直接用NumPy索引X_permuted[:, col_idx] = X_permuted[shuffling_idx, col_idx]进行就地置换c. 在置换后的数据上评估得分，使用
_weights_scorer处理样本权重 -
如果得分是字典（多指标情况），则使用
_aggregate_score_dicts聚合；否则转换为NumPy数组返回
子采样的实现依赖于两个工具函数：_generate_indices 生成无放回样本索引，_safe_indexing 按这些索引安全地索引数据（支持NumPy数组、稀疏矩阵、pandas等）。
对于多指标情况（例如scoring返回字典如{'roc_auc': 0.8, 'f1': 0.75}），聚合过程是这样的：
-
_calculate_permutation_scores中，如果scores[0]是字典，则调用_aggregate_score_dicts(scores)将列表中的字典按键聚合（例如，所有重复的'roc_auc'得分聚合成一个数组） -
在
permutation_importance主函数中，如果baseline_score是字典，则为每个指标名称name创建一个结果：使用baseline_score[name]作为基线得分，并从并行结果中提取对应指标的得分数组np.array([scores[col_idx][name] for col_idx in range(X.shape[1])])
这种设计使得排列重要性既支持单指标也支持多指标评估，并且能够高效地重用预测结果（在多指标情况下避免重复计算）。
我们现在看决策边界可视化的核心类 DecisionBoundaryDisplay。
源码路径：sklearn/inspection/_plot/decision_boundary.py - DecisionBoundaryDisplay.__init__()（60-85行）
def __init__(
self,
*,
xx0,
xx1,
n_classes,
response,
multiclass_colors=None,
xlabel=None,
ylabel=None,
):
self.xx0 = xx0
self.xx1 = xx1
self.n_classes = n_classes
self.response = response
self.multiclass_colors = multiclass_colors
self.xlabel = xlabel
self.ylabel = ylabel
这个构造函数存储了绘图所需的所有数据：
-
xx0,xx1：由np.meshgrid生成的网格坐标矩阵 -
n_classes：期望的类别数量（用于多分类颜色映射） -
response：在网格点上预测的响应值（可以是类别标签、概率或决策函数值） -
multiclass_colors：多分类问题的颜色映射策略 -
xlabel,ylabel：坐标轴标签
真正的绘图发生在 plot 方法中。
源码路径：sklearn/inspection/_plot/decision_boundary.py - DecisionBoundaryDisplay.plot()（87-150行，简化关键逻辑）
def plot(self, plot_method="contourf", ax=None, xlabel=None, ylabel=None, **kwargs):
check_matplotlib_support("DecisionBoundaryDisplay.plot")
import matplotlib as mpl
import matplotlib.pyplot as plt
if plot_method not in ("contourf", "contour", "pcolormesh"):
raise ValueError(
"plot_method must be 'contourf', 'contour', or 'pcolormesh'. "
f"Got {plot_method} instead."
)
if ax is None:
_, ax = plt.subplots()
plot_func = getattr(ax, plot_method)
if self.n_classes == 2:
# 二分类情况：直接绘制响应值
self.surface_ = plot_func(self.xx0, self.xx1, self.response, **kwargs)
else: # 多分类情况
# ... 处理颜色映射警告 ...
if self.multiclass_colors is None or isinstance(self.multiclass_colors, str):
# ... 创建颜色映射 ...
if isinstance(self.multiclass_colors, list):
colors = [mpl.colors.to_rgba(color) for color in self.multiclass_colors]
else:
# ... 处理字符串颜色映射 ...
self.multiclass_colors_ = colors
if self.response.ndim == 2: # predict 方法返回类别标签
# 使用ListedColormap绘制离散类别
cmap = mpl.colors.ListedColormap(colors)
self.surface_ = plot_func(self.xx0, self.xx1, self.response, cmap=cmap, **kwargs)
else:
# predict_proba 或 decision_function 返回概率/决策值
if plot_method == "contour":
# 只绘制整数类别值（通过argmax得到硬分类）
self.surface_ = plot_func(
self.xx0,
self.xx1,
self.response.argmax(axis=2),
colors=colors,
**kwargs,
)
else:
# 为每个类别创建单独的colormap（从白色到类别色）
multiclass_cmaps = [
mpl.colors.LinearSegmentedColormap.from_list(
f"colormap_{class_idx}",
[(1.0, 1.0, 1.0, 1.0), (r, g, b, 1.0)],
)
for class_idx, (r, g, b, _) in enumerate(colors)
]
self.surface_ = []
for class_idx, cmap in enumerate(multiclass_cmaps):
# 掩码数组：只显示当前类别的概率，其他类别设为掩码
response = np.ma.array(
self.response[:, :, class_idx],
mask=(self.response.argmax(axis=2) != class_idx),
)
self.surface_.append(
plot_func(self.xx0, self.xx1, response, cmap=cmap, **kwargs)
)
# ... 设置坐标轴标签 ...
self.ax_ = ax
self.figure_ = ax.figure
return self
这个方法实现了决策边界的绘制逻辑：
-
参数校验：确保
plot_method是有效的Matplotlib绘图方法 -
创建图表：如果未提供
ax，则创建新的figure和axes -
二分类情况（
n_classes == 2）：直接使用指定的绘图方法（如contourf）绘制响应值矩阵 -
多分类情况（
n_classes > 2）：a. 处理颜色映射：如果未提供
multiclass_colors，则根据类别数选择默认 colormap（10类以下用'tab10'，超过10类用'gist_rainbow'）；如果是字符串则验证其是否为有效colormap；如果是列表则验证长度和颜色有效性b. 根据响应值的维度分情况处理：
-
如果是类别标签（
response.ndim == 2，形状为网格形状）：直接绘制，使用离散colormap -
如果是概率或决策函数值（
response.ndim == 3，形状为网格形状 + 类别数）：-
如果是
contour方法：只绘制硬分类结果（使用argmax(axis=2)得到预测类别） -
否则（
contourf或pcolormesh）：为每个类别创建一个从白色到该类别颜色的线性colormap，然后绘制每个类别的概率/决策值，并使用掩码数组确保每个类别的绘图只在该类别为最高概率的区域可见
-
-
-
设置坐标轴标签和属性存储
对于多分类概率的可视化，这种方法通过绘制每个类别的“贡献度”来工作：每个类别的概率图只在该类别具有最高预测概率的区域可见，并且透明度从白色（零概率）渐变到纯色（高概率）。当多个类别的概率图叠加时，可见的颜色表示具有最高概率的类别，而颜色的深浅表示该概率的置信度。
下面我们使用时序图展示从估计器到决策边界图的完整调用链。
这种设计使得决策边界可视化既支持硬分类（通过predict），也支持软分类（通过predict_proba或decision_function展示概率不确定性），并且在多分类情况下能够直观地显示每个类别的主导区域和概率分布。
理解了这些模型解释技术，我们现在转向高斯过程，这是一种提供不确定性估计的强大建模方法。
28.6 高斯过程内核 —— 函数空间的“协方差积木”
高斯过程内核定义了函数空间中任意两点之间的协方差，即相似度。内核函数 k(x, x') 不仅给出了两点的协方差，还通过其形式和超参数刻画了函数空间的几何结构。常见内核如 RBF、Matern、RationalQuadratic 通过超参数（如 length_scale）控制平滑度和影响范围。内核可通过加法（Sum）和乘法（Product）组合，构建更复杂的协方差结构。平稳内核仅依赖于输入差异，如 RBF 和 Matern；非平稳如 DotProduct 依赖绝对位置。内核支持解析梯度， enabling 高效的超参数通过对数边际似然最大化进行优化。
我们从内核的基类 Kernel 开始看它如何管理超参数并支持组合操作。
源码路径：sklearn/gaussian_process/kernels.py - Kernel.get_params、set_params、clone_with_theta（150-240行，简化关键部分）
def get_params(self, deep=True):
params = dict()
init_sign = signature(self.__class__.__init__)
for parameter in init_sign.parameters.values():
if parameter.kind != parameter.VAR_KEYWORD and parameter.name != "self":
args.append(parameter.name)
for arg in args:
params[arg] = getattr(self, arg)
return params
def set_params(self, **params):
if not params:
return self
valid_params = self.get_params(deep=True)
for key, value in params.items():
split = key.split("__", 1)
if len(split) > 1:
name, sub_name = split
sub_object = valid_params[name]
sub_object.set_params(**{sub_name: value})
else:
setattr(self, key, value)
return self
def clone_with_theta(self, theta):
cloned = clone(self)
cloned.theta = theta
return cloned
这些方法共同实现了内核的参数管理框架：
-
get_params：内省__init__方法的签名，提取所有非变长关键字参数（即超参数），并获取它们的当前值 -
set_params：支持嵌套参数的设置（如k1__length_scale用于设置组合内核中子内核的参数），通过分割键名并递归设置 -
clone_with_theta：克隆当前内核实例并设置其theta属性（对数空间的超参数向量），这在超参数优化过程中非常有用，因为我们需要在保持内核结构不变的同时更新超参数
theta 和 bounds 属性是超参数优化的核心：它们以对数空间表示非固定超参数，使得搜索空间更加规则（特别是长度尺等参数天然在对数尺度上更合适）。
源码路径：sklearn/gaussian_process/kernels.py - Kernel.theta、Kernel.bounds（260-315行）
@property
def theta(self):
theta = []
params = self.get_params()
for hyperparameter in self.hyperparameters:
if not hyperparameter.fixed:
theta.append(params[hyperparameter.name])
if len(theta) > 0:
return np.log(np.hstack(theta))
else:
return np.array([])
@theta.setter
def theta(self, theta):
params = self.get_params()
i = 0
for hyperparameter in self.hyperparameters:
if hyperparameter.fixed:
continue
if hyperparameter.n_elements > 1:
params[hyperparameter.name] = np.exp(theta[i : i + hyperparameter.n_elements])
i += hyperparameter.n_elements
else:
params[hyperparameter.name] = np.exp(theta[i])
i += 1
if i != len(theta):
raise ValueError(
"theta has not the correct number of entries."
" Should be %d; given are %d" % (i, len(theta))
)
self.set_params(**params)
@property
def bounds(self):
bounds = [
hyperparameter.bounds
for hyperparameter in self.hyperparameters
if not hyperparameter.fixed
]
if len(bounds) > 0:
return np.log(np.vstack(bounds))
else:
return np.array([])
这些属性实现了：
-
theta：返回所有非固定超参数的对数值，连接成一个向量（用于梯度优化） -
theta.setter：接受对数空间的超参数向量，将其转换回原始空间并设置到相应的参数上（标量参数取单个元素，向量参数取连续多个元素） -
bounds：返回所有非固定超参数的对数边界，连接成一个矩阵（用于优化器的约束）
内核的组合操作通过运算符重载实现：__add__、__mul__、__pow__ 分别返回 Sum、Product、Exponentiation 实例。
源码路径：sklearn/gaussian_process/kernels.py - Kernel.__add__、__mul__、__pow__（320-350行）
def __add__(self, b):
if not isinstance(b, Kernel):
return Sum(self, ConstantKernel(b))
return Sum(self, b)
def __rmul__(self, b):
if not isinstance(b, Kernel):
return Product(ConstantKernel(b), self)
return Product(b, self)
def __pow__(self, b):
return Exponentiation(self, b)
这些重载方法使得内核组合变得直观：
-
k1 + k2创建一个Sum内核 -
k1 * k2创建一个Product内核 -
k1 ** p创建一个Exponentiation内核（将基内核提升到幂 p）
我们现在看这些复合内核如何实现它们的 __call__ 方法，特别是梯度的计算。
源码路径：sklearn/gaussian_process/kernels.py - Sum.__call__、Product.__call__（360-390行、410-440行）
def __call__(self, X, Y=None, eval_gradient=False):
if eval_gradient:
K1, K1_gradient = self.k1(X, Y, eval_gradient=True)
K2, K2_gradient = self.k2(X, Y, eval_gradient=True)
return K1 + K2, np.dstack((K1_gradient, K2_gradient))
else:
return self.k1(X, Y) + self.k2(X, Y)
def __call__(self, X, Y=None, eval_gradient=False):
if eval_gradient:
K1, K1_gradient = self.k1(X, Y, eval_gradient=True)
K2, K2_gradient = self.k2(X, Y, eval_gradient=True)
return K1 * K2, np.dstack(
(K1_gradient * K2[:, :, np.newaxis], K2_gradient * K1[:, :, np.newaxis])
)
else:
return self.k1(X, Y) * self.k2(X, Y)
这些方法展示了复合内核梯度传播的链式法则：
-
对于
Sum内核：梯度是各子内核梯度在最后一个维度上的堆叠（np.dstack((K1_gradient, K2_gradient))），因为(k1+k2)的梯度对theta是dk1/dtheta + dk2/dtheta -
对于
Product内核：梯度遵循乘积法则-
对于k1的梯度：
dk1/dtheta * k2（即K1_gradient * K2[:, :, np.newaxis]，其中新增的维度用于广播以匹配k2的形状） -
对于k2的梯度：
k1 * dk2/dtheta（即K2_gradient * K1[:, :, np.newaxis]） -
然后将这两个梯度在最后一个维度上堆叠
-
这种梯度计算方式确保了当我们使用梯度-based优化器（如L-BFGS-B）时，可以高效地计算复合内核超参数的梯度。
我们现在看超参数规范类 Hyperparameter 如何设计以支撑固定/可优化超参数的统一管理。
源码路径：sklearn/gaussian_process/kernels.py - Hyperparameter.__new__（100-140行）
def __new__(cls, name, value_type, bounds, n_elements=1, fixed=None):
if not isinstance(bounds, str) or bounds != "fixed":
bounds = np.atleast_2d(bounds)
if n_elements > 1:
if bounds.shape[0] == 1:
bounds = np.repeat(bounds, n_elements, 0)
elif bounds.shape[0] != n_elements:
raise ValueError(
"Bounds on %s should have either 1 or "
"%d dimensions. Given are %d"
% (name, n_elements, bounds.shape[0])
)
if fixed is None:
fixed = isinstance(bounds, str) and bounds == "fixed"
return super().__new__(cls, name, value_type, bounds, n_elements, fixed)
这个构造方法实现了：
-
边界处理：如果边界不是字符串"fixed"，则转换为至少二维的数组；如果是向量超参数（n_elements>1）且边界只有一行，则复制该行以匹配维度
-
固定状态推断：如果未显式提供
fixed参数，则根据边界是否为字符串"fixed"来推断 -
使用
namedtuple的__new__方法创建实例，并利用__slots__ = ()避免实例字典以节省内存
这种设计使得超参数可以统一表示：无论是标量还是向量，无论是固定还是可优化，都通过 Hyperparameter 实例来描述。内核类可以通过属性如 hyperparameter_length_scale 自动生成这些规范（利用属性名前缀匹配）。
下面我们看几个常用内核的数学实现，以 RBF 为例。
源码路径：sklearn/gaussian_process/kernels.py - RBF.__call__（470-510行，简化关键部分）
def __call__(self, X, Y=None, eval_gradient=False):
X = np.atleast_2d(X)
length_scale = _check_length_scale(X, self.length_scale)
if Y is None:
dists = pdist(X / length_scale, metric="sqeuclidean")
K = np.exp(-0.5 * dists)
K = squareform(K)
np.fill_diagonal(K, 1)
else:
if eval_gradient:
raise ValueError("Gradient can only be evaluated when Y is None.")
dists = cdist(X / length_scale, Y / length_scale, metric="sqeuclidean")
K = np.exp(-0.5 * dists)
if eval_gradient:
if self.hyperparameter_length_scale.fixed:
return K, np.empty((X.shape[0], X.shape[0], 0))
elif not self.anisotropic or length_scale.shape[0] == 1:
K_gradient = (K * squareform(dists))[:, :, np.newaxis]
return K, K_gradient
elif self.anisotropic:
K_gradient = (X[:, np.newaxis, :] - X[np.newaxis, :, :]) ** 2 / (
length_scale**2
)
K_gradient *= K[..., np.newaxis]
return K, K_gradient
else:
return K
这个实现展示了：
-
核心计算：对于输入X，计算标准化后的平方欧氏距离矩阵，然后应用RBF核公式
exp(-0.5 * d^2) -
各向异性处理：通过检查
length_scale是否为可迭代且长度>1来判断；各向异性情况下，长度向量的每个元素对应一个特征维度 -
梯度计算（当
eval_gradient=True且 Y=None）：-
各向同性情况：梯度是
K * d^2在最后一个维度上的展开，因为d/dl [exp(-d^2/(2l^2))] = exp(-d^2/(2l^2)) * (d^2/l^3) = K * (d^2/l^2) * (1/l)，但代码中利用了squareform(dists)来重构完整距离矩阵 -
各向异性情况：梯度张量的形状是
(n_samples, n_samples, n_features)，因为对于每对样本点(i,j)和每个特征维度d，我们有∂k/∂l_d = k * ( (x_i,d - x_j,d)^2 / l_d^3 )，代码中通过(X[:, np.newaxis, :] - X[np.newaxis, :, :]) ** 2 / (length_scale**2)计算了所有样本对和特征维度的平方差除以长度平方，然后乘以核值K[..., np.newaxis]得到梯度
-
这种梯度张量的形状 (n, n, n_features) 正是为了后续在复合内核中进行梯度堆叠而设计的：当我们有多个内核时，它们的梯度张量可以在最后一个维度上连接，形成 (n, n, total_n_dims) 的总梯度张量。
同样，我们看Matern核在特殊nu值下的解析式，以避免昂贵的Bessel函数计算。
源码路径：sklearn/gaussian_process/kernels.py - Matern.__call__ 的特殊nu值处理（520-580行的核心部分）
def __call__(self, X, Y=None, eval_gradient=False):
# ... 长度尺检查 ...
if Y is None:
dists = pdist(X / length_scale, metric="euclidean")
else:
if eval_gradient:
raise ValueError("Gradient can only be evaluated when Y is None.")
dists = cdist(X / length_scale, Y / length_scale, metric="euclidean")
if self.nu == 0.5:
K = np.exp(-dists)
elif self.nu == 1.5:
K = dists * math.sqrt(3)
K = (1.0 + K) * np.exp(-K)
elif self.nu == 2.5:
K = dists * math.sqrt(5)
K = (1.0 + K + K**2 / 3.0) * np.exp(-K)
elif self.nu == np.inf:
K = np.exp(-(dists**2) / 2.0)
else: # 一般情况：昂贵的Bessel函数
# ... 使用修改Bessel函数kv和伽马函数 ...
# ... 转换为方阵并填充对角线 ...
if eval_gradient:
# ... 特殊nu值的梯度公式 ...
if self.nu == 0.5:
# ... 梯度涉及归一化距离 ...
elif self.nu == 1.5:
K_gradient = 3 * D * np.exp(-np.sqrt(3 * D.sum(-1)))[..., np.newaxis]
elif self.nu == 2.5:
tmp = np.sqrt(5 * D.sum(-1))[..., np.newaxis]
K_gradient = 5.0 / 3.0 * D * (tmp + 1) * np.exp(-tmp)
elif self.nu == np.inf:
K_gradient = D * K[..., np.newaxis]
else:
# ... 数值近似梯度 ...
else:
return K
这个实现展示了Matern核如何在nu=0.5, 1.5, 2.5, inf时提供解析式：
-
nu=0.5: 绝对指数核
exp(-d)，对应Once differentiable -
nu=1.5:
(1 + sqrt(3)d) * exp(-sqrt(3)d) -
nu=2.5:
(1 + sqrt(5)d + (5/3)d^2) * exp(-sqrt(5)d) -
nu=inf: 高斯核
exp(-d^2/2)，即RBF核（不含长度尺因子，因为这里的d已经是标准化距离） -
其他nu值: 需要修改Bessel函数
kv，计算成本高约10倍
对应的梯度公式也被推导出来以保持一致性：
-
nu=1.5: 梯度比例为
3 * D * exp(-sqrt(3 * sum(D)))，其中D是平方距离矩阵 -
nu=2.5: 梯度比例为
(5/3) * D * (sqrt(5 * sum(D)) + 1) * exp(-sqrt(5 * sum(D))) -
nu=inf: 梯度是
D * K，与RBF核的一致（因为nu=inf时Matern退化为RBF）
这种设计使得在常用平滑度设置下，Matern核的计算既快又准确，同时保持了与梯度-based超参数优化的兼容性。
下面我们看有理二次核（RationalQuadratic）作为RBF尺度的混合。
源码路径：sklearn/gaussian_process/kernels.py - RationalQuadratic.__call__（590-650行，简化关键部分）
def __call__(self, X, Y=None, eval_gradient=False):
X = np.atleast_2d(X)
if Y is None:
dists = squareform(pdist(X, metric="sqeuclidean"))
tmp = dists / (2 * self.alpha * self.length_scale**2)
base = 1 + tmp
K = base**-self.alpha
np.fill_diagonal(K, 1)
else:
if eval_gradient:
raise ValueError("Gradient can only be evaluated when Y is None.")
dists = cdist(X, Y, metric="sqeuclidean")
K = (1 + dists / (2 * self.alpha * self.length_scale**2)) ** -self.alpha
if eval_gradient:
# 长度尺梯度
if not self.hyperparameter_length_scale.fixed:
length_scale_gradient = dists * K / (self.length_scale**2 * base)
length_scale_gradient = length_scale_gradient[:, :, np.newaxis]
else:
length_scale_gradient = np.empty((K.shape[0], K.shape[1], 0))
# alpha梯度
if not self.hyperparameter_alpha.fixed:
alpha_gradient = K * (
-self.alpha * np.log(base)
+ dists / (2 * self.length_scale**2 * base)
)
alpha_gradient = alpha_gradient[:, :, np.newaxis]
else:
alpha_gradient = np.empty((K.shape[0], K.shape[1], 0))
return K, np.dstack((alpha_gradient, length_scale_gradient))
else:
return K
这个实现展示了有理二次核如何作为RBF的尺度混合：
-
核心公式：
(1 + d^2/(2αl^2))^(-α)可以看作是无限个RBF核的均匀混合，其中每个RBF核的长度尺服从逆伽马分布 -
当α→∞时，有理二次核退化为RBF核（使用极限：
(1 + x/α)^(-α) → exp(-x)） -
当α→0+时，核变得更不平滑，对异常值更敏感
-
梯度计算分为两部分：
-
长度尺梯度：涉及
d^2 * K / (l^2 * base)，反映了改变长度尺对核值的影响 -
alpha梯度：涉及
K * (-α * log(base) + d^2/(2l^2 * base))，来自于对公式求导
-
这种内核在建模具有多尺度特征的函数时非常有用，因为它能够自适应地捕捉不同长度尺的相关性。
理解了高斯过程内核，我们现在看高斯过程回归如何利用这些内核进行不确定性量化预测。
28.7 高斯过程回归与分类 —— 概率预测的“前哨站”
高斯过程回归 (GPR) 通过后验高斯过程预测均值和方差，提供完整的不确定性估计。均值来自后验均值，方差包含数据噪声和模型不确定性两部分。对数边际似然用于超参数优化，支持梯度-based 优化器如 L-BFGS-B。预测时可返回标准差或协方差，支持风险敏感决策。
高斯过程分类 (GPC) 使用 Laplace 近似将非高斯后验转换为高斯，以便进行解析推断。二分类中，后验模式通过牛顿法迭代求得，方差通过 Hessian 逆近似。多分类通过 OvR 或 OvO 策略分解为多个二分类问题。预测概率通过误差函数的线性组合近似 logistic sigmoid 积分。
我们先看GPR的核心：训练过程中的超参数优化和后验计算。
源码路径：sklearn/gaussian_process/_gpr.py - GaussianProcessRegressor.fit()（150-210行，简化关键部分）
def fit(self, X, y):
if self.kernel is None:
self.kernel_ = C(1.0, constant_value_bounds="fixed") * RBF(
1.0, length_scale_bounds="fixed"
)
else:
self.kernel_ = clone(self.kernel)
self._rng = check_random_state(self.random_state)
if self.kernel_.requires_vector_input:
dtype, ensure_2d = "numeric", True
else:
dtype, ensure_2d = None, False
X, y = validate_data(
self,
X,
y,
multi_output=True,
y_numeric=True,
ensure_2d=ensure_2d,
dtype=dtype,
)
if self.normalize_y:
self._y_train_mean = np.mean(y, axis=0)
self._y_train_std = _handle_zeros_in_scale(np.std(y, axis=0), copy=False)
y = (y - self._y_train_mean) / self._y_train_std
self.X_train_ = np.copy(X) if self.copy_X_train else X
self.y_train_ = np.copy(y) if self.copy_X_train else y
if self.optimizer is not None and self.kernel_.n_dims > 0:
def obj_func(theta, eval_gradient=True):
if eval_gradient:
lml, grad = self.log_marginal_likelihood(
theta, eval_gradient=True, clone_kernel=False
)
return -lml, -grad
else:
return -self.log_marginal_likelihood(theta, clone_kernel=False)
optima = [self._constrained_optimization(obj_func, self.kernel_.theta, self.kernel_.bounds)]
if self.n_restarts_optimizer > 0:
if not np.isfinite(self.kernel_.bounds).all():
raise ValueError(
"Multiple optimizer restarts (n_restarts_optimizer>0) "
"requires that all bounds are finite."
)
bounds = self.kernel_.bounds
for iteration in range(self.n_restarts_optimizer):
theta_initial = self._rng.uniform(bounds[:, 0], bounds[:, 1])
optima.append(
self._constrained_optimization(obj_func, theta_initial, bounds)
)
lml_values = list(map(itemgetter(1), optima))
self.kernel_.theta = optima[np.argmin(lml_values)][0]
self.kernel_._check_bounds_params()
self.log_marginal_likelihood_value_ = -np.min(lml_values)
else:
self.log_marginal_likelihood_value_ = self.log_marginal_likelihood(
self.kernel_.theta, clone_kernel=False
)
K = self.kernel_(self.X_train_)
K[np.diag_indices_from(K)] += self.alpha
try:
self.L_ = cholesky(K, lower=GPR_CHOLESKY_LOWER, check_finite=False)
except np.linalg.LinAlgError as exc:
exc.args = (
(
f"The kernel, {self.kernel_}, is not returning a positive "
"definite matrix. Try gradually increasing the 'alpha' "
"parameter of your GaussianProcessRegressor estimator."
),
) + exc.args
raise
self.alpha_ = cho_solve(
(self.L_, GPR_CHOLESKY_LOWER),
self.y_train_,
check_finite=False,
)
return self
这个训练过程实现了高斯过程回归的核心步骤：
-
核初始化：如果未提供核，则使用默认的
ConstantKernel * RBF组合（相当于在RBF核上加一个常数项，建模非零均值过程） -
数据验证和预处理：根据核的要求设置数据类型和维度；如果
normalize_y=True，则将目标值标准化为零均值单位方差（训练后会反向变换） -
超参数优化：如果核有可优化参数且优化器未设为None
a. 定义目标函数：负对数边际似然（因为我们要最小化它）
b. 从核当前参数开始的优化
c. 如果
n_restarts_optimizer > 0，则从对数均匀分布中采样初始点进行重启（避免局部最优）d. 选择所有运行中负对数边际似然最小（即边际似然最大）的参数
e. 检查优化后的参数是否接近边界（可能表明边界太紧）
-
后验计算准备：
a. 计算训练数据上的核矩阵
K = kernel_(X_train_)b. 添加噪声方差到对角线：
K[diag] += alpha（这里的alpha建模观测噪声）c. 进行Cholesky分解：
L = cholesky(K)，得到下三角矩阵满足K = L L^Td. 计算双系数：
alpha = L^T \ (L \ y)，利用三角求解高效计算
这种分解使得预测变得非常高效，因为我们已经将问题分解为两个三角求解步骤。
我们现在看GPR的预测过程，这是它提供不确定性估计的关键。
源码路径：sklearn/gaussian_process/_gpr.py - GaussianProcessRegressor.predict()（260-320行，简化关键部分）
def predict(self, X, return_std=False, return_cov=False):
if return_std and return_cov:
raise RuntimeError(
"At most one of return_std or return_cov can be requested."
)
if self.kernel is None or self.kernel.requires_vector_input:
dtype, ensure_2d = "numeric", True
else:
dtype, ensure_2d = None, False
X = validate_data(self, X, ensure_2d=ensure_2d, dtype=dtype, reset=False)
if not hasattr(self, "X_train_"): # 未拟合：基于先验预测
if self.kernel is None:
kernel = C(1.0, constant_value_bounds="fixed") * RBF(
1.0, length_scale_bounds="fixed"
)
else:
kernel = self.kernel
y_mean = np.zeros(shape=(X.shape[0], n_targets)).squeeze()
if return_cov:
y_cov = kernel(X)
# ... 处理多输出 ...
return y_mean, y_cov
elif return_std:
y_var = kernel.diag(X)
# ... 处理多输出 ...
return y_mean, np.sqrt(y_var)
else:
return y_mean
else: # 已拟合：基于后验预测
K_trans = self.kernel_(X, self.X_train_) # K(X_test, X_train)
y_mean = K_trans @ self.alpha_ # K_test_train * alpha
y_mean = self._y_train_std * y_mean + self._y_train_mean # 反标准化
if y_mean.ndim > 1 and y_mean.shape[1] == 1:
y_mean = np.squeeze(y_mean, axis=1)
if not return_cov and not return_std:
return y_mean
V = solve_triangular(
self.L_, K_trans.T, lower=GPR_CHOLESKY_LOWER, check_finite=False
) # V = L^T \ K_test_train
if return_cov:
y_cov = self.kernel_(X) - V.T @ V # K_test_test - V^T V
y_cov = np.outer(y_cov, self._y_train_std**2).reshape(*y_cov.shape, -1)
if y_cov.shape[2] == 1:
y_cov = np.squeeze(y_cov, axis=2)
return y_mean, y_cov
else: # return_std
y_var = self.kernel_.diag(X).copy() # prior variance at test points
y_var -= np.einsum("ij,ji->i", V.T, V) # 减去 V^T V 的对角线（等价于 trace(V^T V)但更高效）
y_var_negative = y_var < 0
if np.any(y_var_negative):
warnings.warn(
"Predicted variances smaller than 0. "
"Setting those variances to 0."
)
y_var[y_var_negative] = 0.0
y_var = np.outer(y_var, self._y_train_std**2).reshape(*y_var.shape, -1)
if y_var.shape[1] == 1:
y_var = np.squeeze(y_var, axis=1)
return y_mean, np.sqrt(y_var)
这个预测过程清晰地区分了未拟合和已拟合两种情况：
-
未拟合时：直接使用高斯过程先验，预测均值为零（或先验均值），方差仅来自核函数
-
已拟合时：利用已经计算好的后验量
a. 均值预测：
y_mean = K(X_test, X_train) @ alpha，其中alpha包含了训练数据的信息b. 不确定性量化：
-
协方差预测：
y_cov = K(X_test, X_test) - V^T V，其中V = L^T \ K(X_test, X_train)^T这个公式来源于：后验协方差 = 先验协方差 - 由于观测数据减少的不确定性
其中
V^T V正是由于观测数据带来的不确定性减少量 -
方差预测（更高效）：先得到先验方差
kernel_.diag(X)，然后减去np.einsum("ij,ji->i", V.T, V)（这是V^T V的对角线之和，更高效因为不需要构建完整矩阵）这个技巧利用了：
diag(V^T V) = sum_j V_ij V_ji = sum_j V_ij^2（如果V是实矩阵），但更一般地，对于任何矩阵A，diag(A^T A)可以通过einsum("ij,ji->i", A, A)高效计算
-
这种方法不仅给出了预测均值（最佳猜测），还提供了方差或完整协方差，使得我们能够量化预测的不确定性——这正是高斯过程的核心优势之一。
对数边际似然在超参数优化中起着核心作用，它的计算和梯度如下所示。
源码路径：sklearn/gaussian_process/_gpr.py - GaussianProcessRegressor.log_marginal_likelihood()（340-400行，简化关键部分）
def log_marginal_likelihood(
self, theta=None, eval_gradient=False, clone_kernel=True
):
if theta is None:
if eval_gradient:
raise ValueError("Gradient can only be evaluated for theta!=None")
return self.log_marginal_likelihood_value_
if clone_kernel:
kernel = self.kernel_.clone_with_theta(theta)
else:
kernel = self.kernel_
kernel.theta = theta
if eval_gradient:
K, K_gradient = kernel(self.X_train_, eval_gradient=True)
else:
K = kernel(self.X_train_)
K[np.diag_indices_from(K)] += self.alpha
try:
L = cholesky(K, lower=GPR_CHOLESKY_LOWER, check_finite=False)
except np.linalg.LinAlgError:
return (-np.inf, np.zeros_like(theta)) if eval_gradient else -np.inf
y_train = self.y_train_
if y_train.ndim == 1:
y_train = y_train[:, np.newaxis]
alpha = cho_solve((L, GPR_CHOLESKY_LOWER), y_train, check_finite=False)
log_likelihood_dims = -0.5 * np.einsum("ik,ik->k", y_train, alpha)
log_likelihood_dims -= np.log(np.diag(L)).sum()
log_likelihood_dims -= K.shape[0] / 2 * np.log(2 * np.pi)
log_likelihood = log_likelihood_dims.sum(axis=-1)
if eval_gradient:
inner_term = np.einsum("ik,jk->ijk", alpha, alpha)
K_inv = cho_solve(
(L, GPR_CHOLESKY_LOWER), np.eye(K.shape[0]), check_finite=False
)
inner_term -= K_inv[..., np.newaxis]
log_likelihood_gradient_dims = 0.5 * np.einsum(
"ijl,jik->kl", inner_term, K_gradient
)
log_likelihood_gradient = log_likelihood_gradient_dims.sum(axis=-1)
if eval_gradient:
return log_likelihood, log_likelihood_gradient
else:
return log_likelihood
这个函数实现了对数边际似然及其梯度的高效计算：
-
核矩阵计算：
K = kernel(X_train_)（如果需要梯度则同时计算核对超参数的梯度） -
添加噪声：
K[diag] += alpha -
Cholesky分解：
L = cholesky(K)，若失败则返回负无穷（表示不是正定矩阵） -
计算双系数：
alpha = L^T \ (L \ y_train) -
对数边际似然（对数似然）：
-0.5 * y_train^T @ alpha：数据拟合项- sum(log(diag(L)))： Occam's razor项，惩罚复杂核- n_samples / 2 * log(2*pi)：归一化常数 -
梯度计算（如果需要）：
a. 计算
inner_term = alpha @ alpha^T（样本协方差矩阵的估计）b. 计算逆核矩阵：
K_inv = L^T \ (L \ I)通过 Cholesky 分解高效得到c. 计算
inner_term -= K_inv[..., np.newaxis]得到alpha @ alpha^T - K^{-1}d. 梯度是
0.5 * trace((alpha @ alpha^T - K^{-1}) @ K_gradient)，通过einsum("ijl,jik->kl", inner_term, K_gradient)高效计算（避免显式矩阵乘法）
这种实现不仅在数学上正确，而且在计算上非常高效，因为它充分利用了Cholesky分解的中间结果。
我们现在看二分类高斯过程分类中的拉普拉斯近似，这是如何处理非似然问题的关键。
源码路径：sklearn/gaussian_process/_gpc.py - _BinaryGaussianProcessClassifierLaplace.__init__() 和 fit()（100-140行，简化关键部分）
def __init__(
self,
kernel=None,
*,
optimizer="fmin_l_bfgs_b",
n_restarts_optimizer=0,
max_iter_predict=100,
warm_start=False,
copy_X_train=True,
random_state=None,
):
self.kernel = kernel
self.optimizer = optimizer
self.n_restarts_optimizer = n_restarts_optimizer
self.max_iter_predict = max_iter_predict
self.warm_start = warm_start
self.copy_X_train = copy_X_train
self.random_state = random_state
def fit(self, X, y):
if self.kernel is None:
self.kernel_ = C(1.0, constant_value_bounds="fixed") * RBF(
1.0, length_scale_bounds="fixed"
)
else:
self.kernel_ = clone(self.kernel)
self.rng = check_random_state(self.random_state)
self.X_train_ = np.copy(X) if self.copy_X_train else X
label_encoder = LabelEncoder()
self.y_train_ = label_encoder.fit_transform(y)
self.classes_ = label_encoder.classes_
if self.classes_.size > 2:
raise ValueError(
"%s supports only binary classification. y contains classes %s"
% (self.__class__.__name__, self.classes_)
)
elif self.classes_.size == 1:
raise ValueError(
"{0:s} requires 2 classes; got {1:d} class".format(
self.__class__.__name__, self.classes_.size
)
)
if self.optimizer is not None and self.kernel_.n_dims > 0:
def obj_func(theta, eval_gradient=True):
if eval_gradient:
lml, grad = self.log_marginal_likelihood(
theta, eval_gradient=True, clone_kernel=False
)
return -lml, -grad
else:
return -self.log_marginal_likelihood(theta, clone_kernel=False)
optima = [self._constrained_optimization(obj_func, self.kernel_.theta, self.kernel_.bounds)]
if self.n_restarts_optimizer > 0:
if not np.isfinite(self.kernel_.bounds).all():
raise ValueError(
"Multiple optimizer restarts (n_restarts_optimizer>0) "
"requires that all bounds are finite."
)
bounds = self.kernel_.bounds
for iteration in range(self.n_restarts_optimizer):
theta_initial = np.exp(self.rng.uniform(bounds[:, 0], bounds[:, 1]))
optima.append(
self._constrained_optimization(obj_func, theta_initial, bounds)
)
lml_values = list(map(itemgetter(1), optima))
self.kernel_.theta = optima[np.argmin(lml_values)][0]
self.kernel_._check_bounds_params()
self.log_marginal_likelihood_value_ = -np.min(lml_values)
else:
self.log_marginal_likelihood_value_ = self.log_marginal_likelihood(
self.kernel_.theta
)
K = self.kernel_(self.X_train_)
_, (self.pi_, self.W_sr_, self.L_, _, _) = self._posterior_mode(
K, return_temporaries=True
)
return self
这个训练过程与GPR类似，但有几个关键区别：
-
只支持二分类：通过标签编码并检查类别数
-
超参数优化目标仍然是对数边际似然，但计算方式不同（使用拉普拉斯近似）
-
训练结束后，不仅计算了超参数，还通过
_posterior_mode方法得到了：-
self.pi_：训练点的后验类别概率（正类的概率） -
self.W_sr_：Hessian矩阵的平方根（用于后续预测） -
self.L_：核矩阵加噪声的Cholesky分解
-
-
存储了对数边际似然的值用于后续参考
真正的拉普拉斯近似发生在 _posterior_mode 方法中，它使用牛顿法迭代寻找后验模式。
源码路径：sklearn/gaussian_process/_gpc.py - _BinaryGaussianProcessClassifierLaplace._posterior_mode()（简化核心迭代循环）
def _posterior_mode(self, K, return_temporaries=False):
if self.warm_start and hasattr(self, "f_cached") and self.f_cached.shape == self.y_train_.shape:
f = self.f_cached
else:
f = np.zeros_like(self.y_train_, dtype=np.float64)
log_marginal_likelihood = -np.inf
for _ in range(self.max_iter_predict):
pi = expit(f) # sigmoid函数
W = pi * (1 - pi) # 对数似然的二阶导数（权重）
W_sr = np.sqrt(W)
W_sr_K = W_sr[:, np.newaxis] * K
B = np.eye(W.shape[0]) + W_sr_K * W_sr
L = cholesky(B, lower=True)
b = W * f + (self.y_train_ - pi)
a = b - W_sr * cho_solve((L, True), W_sr_K.dot(b))
f = K.dot(a)
lml = (
-0.5 * a.T.dot(f)
- np.log1p(np.exp(-(self.y_train_ * 2 - 1) * f)).sum()
- np.log(np.diag(L)).sum()
)
if lml - log_marginal_likelihood < 1e-10:
break
log_marginal_likelihood = lml
self.f_cached = f
if return_temporaries:
return log_marginal_likelihood, (pi, W_sr, L, b, a)
else:
return log_marginal_likelihood
这个实现直接遵循了Rasmussen & Williams《高斯过程机器学习》一书中的Algorithm 3.1：
-
初始化 latent 函数向量
f（ warm_start 时使用上次解） -
牛顿法迭代：
a. 计算后验类别概率：
pi = sigmoid(f)这是连接函数的应用b. 计算权重：
W = pi * (1 - pi)这是对数似然在latent函数上的二阶导数c. 构建 Hessian 近似：
B = I + W_sr K W_sr其中W_sr = sqrt(W)d. Cholesky分解得到
L使得B = L L^Te. 计算修正项：
b = W f + (y - pi)和a = b - W_sr L^T \ (L \ (W_sr K b))f. 更新 latent 函数：
f = K ag. 计算当前对数边际似然作为收敛标准
-
迭代直到对数边际似然收敛（变化小于1e-10）或达到最大迭代次数
训练结束后，我们有了后验模式 f（存储在 self.f_cached 中），以及中间量 pi_、W_sr_、L_，这些将用于预测。
现在我们看预测过程，特别是如何得到概率预测。
源码路径：sklearn/gaussian_process/_gpc.py - _BinaryGaussianProcessClassifierLaplace.predict_proba()（220-260行，简化关键部分）
def predict_proba(self, X):
latent_mean, latent_var = self.latent_mean_and_variance(X)
alpha = 1 / (2 * latent_var)
gamma = LAMBDAS * latent_mean
integrals = (
np.sqrt(np.pi / alpha)
* erf(gamma * np.sqrt(alpha / (alpha + LAMBDAS**2)))
/ (2 * np.sqrt(latent_var * 2 * np.pi))
)
pi_star = (COEFS * integrals).sum(axis=0) + 0.5 * COEFS.sum()
return np.vstack((1 - pi_star, pi_star)).T
这个预测过程展示了拉普拉斯近似如何工作：
-
先计算latent函数在测试点上的均值和方差
-
latent_mean = K_star^T (y_train - pi_)其中K_star = K(X_train, X_test) -
latent_var = kernel_.diag(X) - diag(V^T V)其中V = L^T \ (W_sr[:, np.newaxis] * K_star)^T
-
-
然后近似计算
E[sigma(f*)]其中sigma是logistic sigmoid，f* ~ N(latent_mean, latent_var) -
关键近似：将logistic sigmoid近似为5个误差函数的线性组合
-
sigma(x) ≈ sum_{i=1}^5 COEFS[i] * erf(LAMBDAS[i] * x) -
其中
LAMBDAS和COEFS是通过最小二乘法预先计算好的常数（在文件顶部给出）
-
-
利用高斯积分的性质：对于
z ~ N(μ, σ^2)，有E[erf(a z)] = erf( a μ / sqrt(1 + 2 a^2 σ^2) )-
因此每个项的期望是：
COEFS[i] * erf( LAMBDAS[i] * latent_mean / sqrt(1 + 2 * LAMBDAS[i]^2 * latent_var) ) -
代码中通过重排得到等价形式：
sqrt(π/alpha) * erf(gamma * sqrt(alpha/(alpha + LAMBDAS^2))) / (2 * sqrt(latent_var * 2π))其中alpha = 1/(2*latent_var)和gamma = LAMBDAS * latent_mean
-
-
最后加上常数项
0.5 * sum(COEFS)来校正近似
这种方法的优势在于：
-
误差函数
erf有解析形式，且其积分在高斯分布下也有解析解 -
避免了需要数值积分来计算
E[sigma(f*)] -
提供了在整个实数域上都很准确的sigmoid近似（最大误差小于1e-5）
最后，我们看多分类是如何扩展的。
源码路径：sklearn/gaussian_process/_gpc.py - GaussianProcessClassifier.fit() 和 kernel_ 属性（400-490行，简化关键部分）
def __init__(
self,
kernel=None,
*,
optimizer="fmin_l_bfgs_b",
n_restarts_optimizer=0,
max_iter_predict=100,
warm_start=False,
copy_X_train=True,
random_state=None,
multi_class="one_vs_rest",
n_jobs=None,
):
self.kernel = kernel
self.optimizer = optimizer
self.n_restarts_optimizer = n_restarts_optimizer
self.max_iter_predict = max_iter_predict
self.warm_start = warm_start
self.copy_X_train = copy_X_train
self.random_state = random_state
self.multi_class = multi_class
self.n_jobs = n_jobs
def fit(self, X, y):
if isinstance(self.kernel, CompoundKernel):
raise ValueError("kernel cannot be a CompoundKernel")
if self.kernel is None or self.kernel.requires_vector_input:
X, y = validate_data(self, X, y, multi_output=False, ensure_2d=True, dtype="numeric")
else:
X, y = validate_data(self, X, y, multi_output=False, ensure_2d=False, dtype=None)
self.base_estimator_ = _BinaryGaussianProcessClassifierLaplace(
kernel=self.kernel,
optimizer=self.optimizer,
n_restarts_optimizer=self.n_restarts_optimizer,
max_iter_predict=self.max_iter_predict,
warm_start=self.warm_start,
copy_X_train=self.copy_X_train,
random_state=self.random_state,
)
self.classes_ = np.unique(y)
self.n_classes_ = self.classes_.size
if self.n_classes_ == 1:
raise ValueError(
"GaussianProcessClassifier requires 2 or more "
"distinct classes; got %d class (only class %s "
"is present)" % (self.n_classes_, self.classes_[0])
)
if self.n_classes_ > 2:
if self.multi_class == "one_vs_rest":
self.base_estimator_ = OneVsRestClassifier(
self.base_estimator_, n_jobs=self.n_jobs
)
elif self.multi_class == "one_vs_one":
self.base_estimator_ = OneVsOneClassifier(
self.base_estimator_, n_jobs=self.n_jobs
)
else:
raise ValueError("Unknown multi-class mode %s" % self.multi_class)
self.base_estimator_.fit(X, y)
if self.n_classes_ > 2:
self.log_marginal_likelihood_value_ = np.mean(
[
estimator.log_marginal_likelihood()
for estimator in self.base_estimator_.estimators_
]
)
else:
self.log_marginal_likelihood_value_ = (
self.base_estimator_.log_marginal_likelihood()
)
return self
@property
def kernel_(self):
if self.n_classes_ == 2:
return self.base_estimator_.kernel_
else:
return CompoundKernel(
[estimator.kernel_ for estimator in self.base_estimator_.estimators_]
)
这个多分类扩展采用了经典的一对剩余（OvR）或一对一（OvO）策略：
-
二分类直接使用底层的
_BinaryGaussianProcessClassifierLaplace -
多分类情况下：
-
OvR：为每个类别训练一个二分类器，将该类作为正类，其余所有类作为负类
-
OvO：为每对类别训练一个二分类器
-
-
超参数优化：在OvR/OvO框架下，每个二分类器独立优化其自身的超参数（通过最大化各自的对数边际似然）
-
预测：将所有二分类器的预测结果结合（OvR使用概率平均或投票，OvO使用投票）
-
注意：OvO不支持概率预测，因为缺乏一致的方式将 pairwise 概率转换为 multi-class 概率
-
核属性：在多分类情况下，
kernel_返回一个CompoundKernel，包含所有二分类器使用的核（这样外部代码仍然可以访问到所使用的核结构）
这种设计使得高斯过程分类既能够保持二分类情况下的理论优雅（拉普拉斯近似），又能够通过经典的多分类简化策略扩展到多分类问题，同时在可能的情况下保持概率输出。
28.8 设计中的取舍
在模型解释和高斯过程中，设计者们在易用性、性能和功能完整性之间做出了诸多权衡。以下是一些关键取舍的解释。
-
为什么部分依赖计算提供两种方法（递归法和暴力法）？
递归法仅针对树模型（如GBDT、随机森林）专门优化，通过一次遍历树结构就能计算所有网格点的平均预测，因此在支持的模型上快得多。然而，它有局限：不支持个体条件期望（ICE），因为它内在地计算了所有样本的平均值；也不支持样本权重，因为它使用的是训练时的样本分布。暴力法虽然计算开销大（对于每个网格点需要重新预测全部样本），但它是通用的——任何实现标准预测接口的估计器都可以使用，并且支持ICE和样本权重。这种设计让用户在需要速度时可以选择递归法（前提是模型支持且不需要ICE或权重），而在需要更灵活性时可以退回到暴力法。
-
高斯过程内核为什么选择使用运算符重载来组合内核？
运算符重载（如
+、*、**）使得内核组合变得直观和简洁，用户可以像处理数学表达式一样构建复杂的协方差结构（例如k1 + k2 * k3 ** 2）。这种设计降低了使用门槛，使得内核工程变得像搭积木一样简单。实现上，这些操作符返回特殊的复合内核类（如Sum、Product、Exponentiation），这些类负责在__call__方法中正确地组合子内核的结果和梯度。虽然这引入了一些额外的类和间接调用，但考虑到内核评估在高斯过程中的频繁发生（特别是在超参数优化期间），这种抽象带来的可用性提升远超过其微小的性能开销。 -
高斯过程分类为什么在多分类情况下使用一对剩余（OvR）而不是真正的多类拉普拉斯近似？
真正的多类拉普拉斯近似需要求解高维非线性方程组来找到后验模式，计算复杂度随着类别数的增加而显著增长。OvR方法将多分类问题分解为多个二分类问题，每个都可以使用已经高度优化的二分类拉普拉斯近似求解器。虽然这不是理论上最优的（因为它没有充分利用类别之间的相关性），但它在实践中非常有效，并且计算上可扩展到较大的类别数。此外，OvR保持了概率预测的能力（通过平均或投票二分类器的概率），而更精确的方法如一对一（OvO）则牺牲了这一特性。这种取舍在实用性和理论纯粹性之间倾向于前者，使得高斯过程分类在实际机器学习任务中成为一个可用的工具。
28.9 动手练习
-
对比PDP的两种计算方法
阅读
sklearn/inspection/_partial_dependence.py中_partial_dependence_recursion(112-138行) 与_partial_dependence_brute(140-210行)。回答问题：
-
递归法要求估计器必须实现什么方法？为什么它不支持
kind='individual'？ -
暴力法中
X_eval = X.copy()的副本用途是什么？为何每次循环都用_safe_assign修改副本？ -
两种方法在
sample_weight处理上有何根本区别？
-
-
排列重要性的并行化与子采样机制
阅读
sklearn/inspection/_permutation_importance.py中permutation_importance(220-290行) 与_calculate_permutation_scores(130-170行)。回答问题：
-
Parallel并行化的粒度是什么（以什么为单位并行）？random_seed如何保证多进程/多线程下的可复现性？ -
max_samples < X.shape[0]时，_generate_indices与_safe_indexing如何协作实现无放回子采样？ -
当
scoring返回字典（多指标）时，_aggregate_score_dicts与后续字典推导式如何聚合结果？
-
-
高斯过程内核的组合与梯度机制
阅读
sklearn/gaussian_process/kernels.py中Kernel.__add__(280-310行)、Sum.__call__(360-390行)、Product.__call__(410-440行) 及Hyperparameter(100-140行)。回答问题：
-
Sum与Product的eval_gradient=True时，梯度张量是如何沿最后一个轴拼接的？形状变化规律是什么？ -
Hyperparameter.__new__中bounds为'fixed'字符串时，如何影响fixed属性与后续theta属性的长度？ -
RBF.__call__中各向异性情况下eval_gradient的计算逻辑（K_gradient形状为何是(n, n, n_features)）？
-
-
GPR 与 GPC 的后验推断对比
对比阅读
sklearn/gaussian_process/_gpr.py中predict(260-320行) 与sklearn/gaussian_process/_gpc.py中_BinaryGaussianProcessClassifierLaplace.latent_mean_and_variance(140-190行) 及predict_proba(220-260行)。回答问题：
-
GPR 的后验均值/方差为何有解析解（Cholesky + 三角求解），而 GPC 需要牛顿法迭代求后验模式？
-
GPC 的
predict_proba为何使用erf(误差函数) 近似积分？LAMBDAS与COEFS代表什么？ -
GPR 的
log_marginal_likelihood梯度公式 (340-400行) 与 GPC 的_posterior_mode中对数边际似然计算 (140-190行) 的核心数学差异在哪里？
-
28.10 本章小结
这一章中我们学习了如何解释“黑盒”机器学习模型以及如何使用高斯过程进行不确定性量化预测。首先，我们探讨了部分依赖图（PDP）和个体条件期望（ICE）这两种互补的技术，它们分别揭示特征的全局趋势和个体异质性，并理解了递归法（专用于树模型的高效通道）和暴力法（通用但计算开销大）的区别。其次，我们掌握了排列重要性的置换检验原理及其并行化和子采样实现，学习了如何通过打乱特征来衡量其对模型性能的贡献。第三，我们了解了决策边界可视化如何通过在特征空间上绘制等高线或热力图来直观展示分类器行为，特别是在多分类情况下的颜色映射策略。最后，我们深入了高斯过程的核心：内核作为协方差积木的设计模式、超参数的对数空间管理以及运算符重载带来的直观组合方式，理解了GPR如何通过Cholesky分解和对数边际似然优化提供均值和方差的预测，以及GPC如何通过拉普拉斯近似和牛顿法求后验模式来处理非似然问题。
以下表格总结了本章的核心概念：
| 概念 | 解释 |
|------|------|
| _grid_from_X | 生成PDP/ICE计算网格，支持分位数、分类唯一值、自定义值 |
| _partial_dependence_recursion | 树模型专用加速算法，加权遍历树结构计算PDP，不支持ICE |
| _partial_dependence_brute | 通用暴力算法，替换特征值重复预测取平均，支持PDP/ICE/样本权重 |
| partial_dependence | 统一入口，自动选择方法，返回Bunch含average/individual/grid_values |
| PartialDependenceDisplay | 可视化类，支持1D曲线/ICE/柱状图、2D等高线/热力图、多分图网格 |
| permutation_importance | 置换重要性主函数，并行计算各特征打乱后得分下降，支持多指标/子采样 |
| _calculate_permutation_scores | 单特征置换核心，支持DataFrame/ndarray、子采样、多次重复 |
| DecisionBoundaryDisplay | 决策边界可视化，支持contourf/contour/pcolormesh、多分类颜色映射 |
| Kernel (基类) | 内核基类，管理超参数(theta/bounds)，运算符重载支持+/ */**组合 |
| RBF/Matern/RationalQuadratic | 常用平稳核，RBF无限光滑，Matern可控光滑度，RQ为RBF尺度混合 |
| Sum/Product/Exponentiation | 内核组合操作符，分别对应协方差加法、乘法、幂运算 |
| GaussianProcessRegressor | GPR回归，Cholesky分解求后验，边际似然优化超参数，预测含均值/方差/协方差 |
| log_marginal_likelihood (GPR) | 边际似然及解析梯度，基于Cholesky因子高效计算，支持多输出 |
| _BinaryGaussianProcessClassifierLaplace | 二分类GPC核心，拉普拉斯近似+牛顿法求后验模式，误差函数近似预测概率 |
| GaussianProcessClassifier | 多分类GPC包装器，支持OvR/OvO，聚合CompoundKernel，不支持OvO预测概率 |
下一章中，我们将学习 scikit-learn 的评估指标体系，了解它如何为分类、回归、聚类等任务提供全面的度量工具，帮助我们客观评估模型性能。
第 29 章 —— scikit-learn metrics 概览 —— 认识这座“模型评估的度量衡博物馆”
29.1 学习目标
-
难度：★★★☆☆（3/5）
-
预备知识：Python 基础、面向对象编程与 Markdown/代码阅读基础
-
理解 metrics 模块的核心定位与模块化组织结构
-
掌握分类评估指标的统一计算骨架：混淆矩阵、MCM 与 precision_recall_fscore_support
-
熟悉回归评估指标的统一目标验证与多输出聚合框架
-
理解 Tweedie 偏差族与 D² 分数族如何统一多种分布假设下的损失计算
-
掌握多标签、多分类场景下的平均策略与零除处理机制
-
了解 Array API 兼容层在指标计算中的应用与设计思想
29.2 生活类比
想象 scikit-learn 的 metrics 模块是一座模型评估的“度量衡博物馆”。这里收藏着各式各样的“测量仪器”，帮助我们量化机器学习模型的表现，就像博物馆里的展品用来衡量历史文物的价值一样。这座博物馆分为若干展厅：
-
分类指标展厅（“裁判组”）：混淆矩阵提供详细的录像回放，显示每个类别的预测正确与否；Precision、Recall 和 F1 分数则像是命中率、覆盖率和综合评分，从不同角度评估模型的准确性和完整性；而 ROC 曲线和 AUC 则相当于全阈值扫描雷达，在所有可能的决策阈值下考察模型的区分能力。
-
回归指标展厅（“尺子箱”）：MAE 和 MSE 提供直尺和卷尺般的绝对误差测量，R² 就像拟合优度仪，衡量模型解释目标变量方差的程度；Tweedie 偏差族更像是一把可调节的弹性尺，通过 power 参数在不同的分布假设之间切换，适用于正态、泊松、伽马等各种场景。
-
聚类指标展厅（“质检员”）：有监督聚类指标如 ARI 和 NMI 对照“标准答案”来评估聚类结果；而无监督指标如轮廓系数则靠“自查内部结构”，在没有真实标签的情况下判断聚类的紧密程度和分离度。
为了确保所有这些仪器都能准确使用，博物馆设有“安检门”：
_check_targets和_check_reg_targets负责统一输入格式、对齐形状并处理缺失值，确保每个待测样本都符合要求。
“计分规则手册”则由 _average_binary_score 和 _prf_divide 等函数构成，它们定义了：
-
如何在多标签和多分类场景下进行微观、宏观、加权平均
-
明确了零除时是返回 0、1、还是 np.nan，或者直接发出警告
就像博物馆里每件仪器都有统一的“使用说明书”（API）和“校准证书”（测试），metrics 模块为模型评估提供了标准化、可扩展且跨后端的度量工具集。它通过 Array API 兼容层像一个“万能转换插头”，让同一套测量逻辑能够插入 NumPy、CuPy、PyTorch 等不同的“电源”上运行，确保无论底层数组库如何变化，度量结果都能保持一致和可靠。
在接下来的章节中，我们将沿用这个类比，深入博物馆的各个展厅，拆解每件“测量仪器”的内部结构与工作原理。
29.3 源码地图
29.4 metrics 模块概览 —— 认识这座“模型评估的度量衡博物馆”
29.4.1 源码路径：sklearn/metrics/__init__.py、sklearn/metrics/cluster/__init__.py、sklearn/metrics/_base.py
29.4.2 metrics 模块的核心定位与设计哲学
我们将 metrics 模块视为模型评估的“度量衡博物馆”：它收藏了用于量化机器学习模型表现的各种“测量仪器”。所有度量函数遵循统一的函数签名 (y_true, y_pred, **kwargs)，支持多分类、多标签场景、样本权重 (sample_weight) 以及零除处理 (zero_division)。该模块采用高度模块化的组织方式，将不同任务类型的评估指标划分到独立的子模块中：_classification.py 处理分类任务，_regression.py 处理回归任务，_ranking.py 处理排序和曲线评估，cluster/ 目录存放聚类评估指标，pairwise.py 实现成对距离和核函数，_scorer.py 提供评分器包装机制，而 _plot/ 则负责可视化 Display 类。这种设计不仅使代码结构清晰，还便于用户根据需求灵活导入和使用特定的评估工具。
29.4.3 公共 API 与模块组织结构
sklearn/metrics/__init__.py 采用 re-export 模式，将上百个评估函数和可视化类进行聚合并统一暴露，构建了对外的统一接口。它从子模块中导入了分类指标（如 accuracy_score, f1_score）、回归指标（如 mean_squared_error, r2_score）、排序指标（如 roc_auc_score, average_precision_score）、聚类指标（如 adjusted_rand_score, silhouette_score）、成对距离函数（如 euclidean_distances, pairwise_distances）、评分器工具（如 make_scorer, check_scoring）以及可视化类（如 ConfusionMatrixDisplay, RocCurveDisplay）。同时，sklearn/metrics/cluster/__init__.py 单独负责导出聚类评估指标，覆盖有监督（如 adjusted_rand_score）、无监督（如 silhouette_score）和双聚类（如 consensus_score）三类场景。
29.4.4 基础工具函数：统一的预处理与平均策略
在度量计算的底层，有一套通用的预处理和平均策略确保了结果的一致性和正确性。
29.4.4.1 源码路径：sklearn/metrics/_base.py
29.4.4.1.1 _average_binary_score
该函数实现了二分类指标在多标签和多分类场景下的平均策略。它支持 micro、macro、weighted 和 samples 四种平均模式。
def _average_binary_score(binary_metric, y_true, y_score, average, sample_weight=None):
# 获取数组命名空间（NumPy/CuPy/PyTorch等）
xp, _, _device = get_namespace_and_device(y_true, y_score, sample_weight)
average_options = (None, "micro", "macro", "weighted", "samples")
if average not in average_options:
raise ValueError("average has to be one of {0}".format(average_options))
y_type = type_of_target(y_true)
if y_type not in ("binary", "multilabel-indicator"):
raise ValueError("{0} format is not supported".format(y_type))
if y_type == "binary":
return binary_metric(y_true, y_score, sample_weight=sample_weight)
check_consistent_length(y_true, y_score, sample_weight)
y_true = check_array(y_true)
y_score = check_array(y_score)
not_average_axis = 1
score_weight = sample_weight
average_weight = None
if average == "micro":
# Micro 平均：将所有标签展平，全局计算
if score_weight is not None:
score_weight = xp.repeat(score_weight, y_true.shape[1])
y_true = _ravel(y_true)
y_score = _ravel(y_score)
elif average == "weighted":
# Weighted 平均：按每个标签的支持度加权
if score_weight is not None:
y_true = xp.asarray(y_true, dtype=score_weight.dtype)
average_weight = xp.sum(
xp.multiply(y_true, xp.reshape(score_weight, (-1, 1))), axis=0
)
else:
average_weight = xp.sum(y_true, axis=0)
if xpx.isclose(
xp.sum(average_weight),
xp.asarray(0, dtype=average_weight.dtype, device=_device),
):
return 0
elif average == "samples":
# Samples 平均：按样本维度计算，再平均
average_weight = score_weight
score_weight = None
not_average_axis = 0
if y_true.ndim == 1:
y_true = xp.reshape(y_true, (-1, 1))
if y_score.ndim == 1:
y_score = xp.reshape(y_score, (-1, 1))
n_classes = y_score.shape[not_average_axis]
score = xp.zeros((n_classes,), device=_device)
for c in range(n_classes):
y_true_c = _ravel(
xp.take(y_true, xp.asarray([c], device=_device), axis=not_average_axis)
)
y_score_c = _ravel(
xp.take(y_score, xp.asarray([c], device=_device), axis=not_average_axis)
)
score[c] = binary_metric(y_true_c, y_score_c, sample_weight=score_weight)
# 最终聚合
if average is not None:
if average_weight is not None:
score[average_weight == 0] = 0
return float(_average(score, weights=average_weight, xp=xp))
else:
return score
代码解析：
该函数首先通过 get_namespace_and_device 获取统一的数组 API 命名空间 xp，这是实现跨后端兼容性的关键。对于 micro 平均，它将多标签矩阵展平为一维数组，相当于把所有标签的预测结果混在一起统一计算；对于 weighted 平均，它计算每个标签的支持度作为权重；对于 samples 平均，它交换了平均维度和权重维度，实现逐样本计算后再平均。循环遍历每个类别调用 binary_metric 计算单标签得分，最后根据平均策略进行聚合。这种设计使得同一个二分类指标函数能无缝扩展到多标签和多分类场景。
29.4.4.1.2 _average_multiclass_ovo_score
该函数基于 Hand & Till (2001) 算法实现了 One-vs-One 的多类别平均策略，用于 ROC AUC 等需要成对比较的指标。
def _average_multiclass_ovo_score(binary_metric, y_true, y_score, average="macro"):
check_consistent_length(y_true, y_score)
y_true_unique = np.unique(y_true)
n_classes = y_true_unique.shape[0]
n_pairs = n_classes * (n_classes - 1) // 2
pair_scores = np.empty(n_pairs)
is_weighted = average == "weighted"
prevalence = np.empty(n_pairs) if is_weighted else None
for ix, (a, b) in enumerate(combinations(y_true_unique, 2)):
a_mask = y_true == a
b_mask = y_true == b
ab_mask = np.logical_or(a_mask, b_mask)
if is_weighted:
prevalence[ix] = np.average(ab_mask)
a_true = a_mask[ab_mask]
b_true = b_mask[ab_mask]
a_true_score = binary_metric(a_true, y_score[ab_mask, a])
b_true_score = binary_metric(b_true, y_score[ab_mask, b])
pair_scores[ix] = (a_true_score + b_true_score) / 2
return np.average(pair_scores, weights=prevalence)
代码解析：
函数遍历所有类别对 (a, b)，构造二分类子问题，分别以 a 和 b 为正类计算二分类指标，再取平均得到该类别对的得分。macro 平均直接对所有类别对得分取均值，weighted 平均则按每个类别对包含的样本比例加权。这种实现避免了显式构造 OvO 编码矩阵，内存效率更高。
29.5 分类评估核心实现 —— 从混淆矩阵到 F-beta 的统一计算骨架
29.5.1 源码路径：sklearn/metrics/_classification.py
29.5.2 混淆矩阵：所有分类指标的“原子基石”
混淆矩阵是所有分类评估指标的基础，它以二维数组的形式记录了每个真实类别被预测为各个类别的数量。在二分类情况下，混淆矩阵的四个元素分别对应真阴性（TN）、假阳性（FP）、假阴性（FN）和真阳性（TP），这些基本计数量是构建 Precision、Recall、F1 等指标的原子材料。
29.5.2.1 confusion_matrix
@validate_params(
{
"y_true": ["array-like", "sparse matrix"],
"y_pred": ["array-like", "sparse matrix"],
"labels": ["array-like", None],
"sample_weight": ["array-like", None],
"normalize": [StrOptions({"true", "pred", "all"}), None],
},
prefer_skip_nested_validation=True,
)
def confusion_matrix(
y_true, y_pred, *, labels=None, sample_weight=None, normalize=None
):
xp, _, device_ = get_namespace_and_device(y_true, y_pred, labels, sample_weight)
y_true = check_array(
y_true,
dtype=None,
ensure_2d=False,
ensure_all_finite=False,
ensure_min_samples=0,
)
y_pred = check_array(
y_pred,
dtype=None,
ensure_2d=False,
ensure_all_finite=False,
ensure_min_samples=0,
)
# 转换为 NumPy 数组以利用 SciPy 的高效稀疏矩阵构造
y_true = _convert_to_numpy(y_true, xp)
y_pred = _convert_to_numpy(y_pred, xp)
if sample_weight is None:
sample_weight = np.ones(y_true.shape[0], dtype=np.int64)
else:
sample_weight = _convert_to_numpy(sample_weight, xp)
if len(sample_weight) > 0:
y_type, y_true, y_pred, sample_weight = _check_targets(
y_true, y_pred, sample_weight
)
else:
y_type, y_true, y_pred, _ = _check_targets(y_true, y_pred)
y_true, y_pred = attach_unique(y_true, y_pred)
if y_type not in ("binary", "multiclass"):
raise ValueError("%s is not supported" % y_type)
if labels is None:
labels = unique_labels(y_true, y_pred)
else:
labels = _convert_to_numpy(labels, xp)
n_labels = labels.size
if n_labels == 0:
raise ValueError("'labels' should contain at least one label.")
elif y_true.size == 0:
return np.zeros((n_labels, n_labels), dtype=int)
elif len(np.intersect1d(y_true, labels)) == 0:
raise ValueError("At least one label specified must be in y_true")
n_labels = labels.size
# 如果标签不是从 0 开始的连续整数，需要转换为索引形式
need_index_conversion = not (
labels.dtype.kind in {"i", "u", "b"}
and np.all(labels == np.arange(n_labels))
and y_true.min() >= 0
and y_pred.min() >= 0
)
if need_index_conversion:
label_to_ind = {label: index for index, label in enumerate(labels)}
y_pred = np.array([label_to_ind.get(label, n_labels + 1) for label in y_pred])
y_true = np.array([label_to_ind.get(label, n_labels + 1) for label in y_true])
# 过滤掉不在 labels 中的样本
ind = np.logical_and(y_pred < n_labels, y_true < n_labels)
if not np.all(ind):
y_pred = y_pred[ind]
y_true = y_true[ind]
sample_weight = sample_weight[ind]
# 选择累加器 dtype 以保证高精度
if sample_weight.dtype.kind in {"i", "u", "b"}:
dtype = np.int64
else:
dtype = np.float32 if str(device_).startswith("mps") else np.float64
# 核心：使用 scipy.sparse.coo_matrix 高效统计加权计数
cm = coo_matrix(
(sample_weight, (y_true, y_pred)),
shape=(n_labels, n_labels),
dtype=dtype,
).toarray()
with np.errstate(all="ignore"):
if normalize == "true":
cm = cm / cm.sum(axis=1, keepdims=True)
elif normalize == "pred":
cm = cm / cm.sum(axis=0, keepdims=True)
elif normalize == "all":
cm = cm / cm.sum()
cm = xpx.nan_to_num(cm)
if cm.shape == (1, 1):
warnings.warn(
(
"A single label was found in 'y_true' and 'y_pred'. For the confusion "
"matrix to have the correct shape, use the 'labels' parameter to pass "
"all known labels."
),
UserWarning,
)
return xp.asarray(cm, device=device_)
代码解析：
该函数是混淆矩阵的核心实现。关键设计点在于：将输入转换为 NumPy 数组后，利用 scipy.sparse.coo_matrix 的构造函数直接从 (row_indices, col_indices) 和 data（样本权重）构建稀疏矩阵，再转为稠密数组。这种方式避免了显式循环，极大提高了大规模数据下的计算效率。标签到索引的映射由 attach_unique 和字典查找完成，支持任意类型的标签。归一化选项 true/pred/all 分别按行、按列、按总和归一化，便于查看条件概率分布。
29.5.3 多标签混淆矩阵：类别/样本维度的扩展
在多标签学习中，每个样本可以同时属于多个类别，传统的混淆矩阵无法直接捕捉这种复杂性。multilabel_confusion_matrix 函数扩展了这一概念。
29.5.3.1 multilabel_confusion_matrix
@validate_params(
{
"y_true": ["array-like", "sparse matrix"],
"y_pred": ["array-like", "sparse matrix"],
"sample_weight": ["array-like", None],
"labels": ["array-like", None],
"samplewise": ["boolean"],
},
prefer_skip_nested_validation=True,
)
def multilabel_confusion_matrix(
y_true, y_pred, *, sample_weight=None, labels=None, samplewise=False
):
y_true, y_pred = attach_unique(y_true, y_pred)
xp, _, device_ = get_namespace_and_device(y_true, y_pred, sample_weight)
y_type, y_true, y_pred, sample_weight = _check_targets(
y_true, y_pred, sample_weight
)
if y_type not in ("binary", "multiclass", "multilabel-indicator"):
raise ValueError("%s is not supported" % y_type)
present_labels = unique_labels(y_true, y_pred)
if labels is None:
labels = present_labels
n_labels = None
else:
labels = xp.asarray(labels, device=device_)
n_labels = labels.shape[0]
labels = xp.concat(
[labels, xpx.setdiff1d(present_labels, labels, assume_unique=True, xp=xp)],
axis=-1,
)
if y_true.ndim == 1:
if samplewise:
raise ValueError(
"Samplewise metrics are not available outside of "
"multilabel classification."
)
le = LabelEncoder()
le.fit(labels)
y_true = le.transform(y_true)
y_pred = le.transform(y_pred)
sorted_labels = le.classes_
tp = y_true == y_pred
tp_bins = y_true[tp]
if sample_weight is not None:
tp_bins_weights = sample_weight[tp]
else:
tp_bins_weights = None
if tp_bins.shape[0]:
tp_sum = _bincount(
tp_bins, weights=tp_bins_weights, minlength=labels.shape[0], xp=xp
)
else:
true_sum = pred_sum = tp_sum = xp.zeros(labels.shape[0])
if y_pred.shape[0]:
pred_sum = _bincount(
y_pred, weights=sample_weight, minlength=labels.shape[0], xp=xp
)
if y_true.shape[0]:
true_sum = _bincount(
y_true, weights=sample_weight, minlength=labels.shape[0], xp=xp
)
indices = xp.searchsorted(sorted_labels, labels[:n_labels])
tp_sum = xp.take(tp_sum, indices, axis=0)
true_sum = xp.take(true_sum, indices, axis=0)
pred_sum = xp.take(pred_sum, indices, axis=0)
else:
sum_axis = 1 if samplewise else 0
if labels.shape != present_labels.shape or xp.any(
xp.not_equal(labels, present_labels)
):
if xp.max(labels) > xp.max(present_labels):
raise ValueError(
"All labels must be in [0, n labels) for "
"multilabel targets. "
"Got %d > %d" % (xp.max(labels), xp.max(present_labels))
)
if xp.min(labels) < 0:
raise ValueError(
"All labels must be in [0, n labels) for "
"multilabel targets. "
"Got %d < 0" % xp.min(labels)
)
if n_labels is not None:
y_true = y_true[:, labels[:n_labels]]
y_pred = y_pred[:, labels[:n_labels]]
if issparse(y_true) or issparse(y_pred):
true_and_pred = y_true.multiply(y_pred)
else:
true_and_pred = xp.multiply(y_true, y_pred)
tp_sum = _count_nonzero(
true_and_pred,
axis=sum_axis,
sample_weight=sample_weight,
xp=xp,
device=device_,
)
pred_sum = _count_nonzero(
y_pred,
axis=sum_axis,
sample_weight=sample_weight,
xp=xp,
device=device_,
)
true_sum = _count_nonzero(
y_true,
axis=sum_axis,
sample_weight=sample_weight,
xp=xp,
device=device_,
)
fp = pred_sum - tp_sum
fn = true_sum - tp_sum
tp = tp_sum
if sample_weight is not None and samplewise:
tp = xp.asarray(tp)
fp = xp.asarray(fp)
fn = xp.asarray(fn)
tn = sample_weight * y_true.shape[1] - tp - fp - fn
elif sample_weight is not None:
tn = xp.sum(sample_weight) - tp - fp - fn
elif samplewise:
tn = y_true.shape[1] - tp - fp - fn
else:
tn = y_true.shape[0] - tp - fp - fn
return xp.reshape(xp.stack([tn, fp, fn, tp]).T, (-1, 2, 2))
浙公网安备 33010602011771号