注意
转到页面底部 下载完整示例代码，或通过 JupyterLite 或 Binder 在浏览器中运行此示例。
使用网格搜索进行模型的统计比较#
本示例展示了如何对使用 GridSearchCV 训练和评估的模型性能进行统计比较。
# Authors: The scikit-learn developers
# SPDX-License-Identifier: BSD-3-Clause
我们将从模拟月牙形数据（类与类之间的理想分隔是非线性的）开始，并向其中添加适量的噪声。数据点属于两个类别中的一个，需根据两个特征进行预测。我们将为每个类别模拟 50 个样本。
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.datasets import make_moons
X, y = make_moons(noise=0.352, random_state=1, n_samples=100)
sns.scatterplot(
x=X[:, 0], y=X[:, 1], hue=y, marker="o", s=25, edgecolor="k", legend=False
).set_title("Data")
plt.show()
我们将比较不同 kernel 参数的 SVC 估计器的性能，以确定哪个超参数选择最能预测我们的模拟数据。我们将使用 RepeatedStratifiedKFold 评估模型性能，即进行 10 次 10 折分层交叉验证，且每次重复使用不同的数据随机化方式。性能将使用 roc_auc_score 进行评估。
from sklearn.model_selection import GridSearchCV, RepeatedStratifiedKFold
from sklearn.svm import SVC
param_grid = [
{"kernel": ["linear"]},
{"kernel": ["poly"], "degree": [2, 3]},
{"kernel": ["rbf"]},
]
svc = SVC(random_state=0)
cv = RepeatedStratifiedKFold(n_splits=10, n_repeats=10, random_state=0)
search = GridSearchCV(estimator=svc, param_grid=param_grid, scoring="roc_auc", cv=cv)
search.fit(X, y)
现在我们可以检查搜索结果，按其 mean_test_score 进行排序
import pandas as pd
results_df = pd.DataFrame(search.cv_results_)
results_df = results_df.sort_values(by=["rank_test_score"])
results_df = results_df.set_index(
results_df["params"].apply(lambda x: "_".join(str(val) for val in x.values()))
).rename_axis("kernel")
results_df[["params", "rank_test_score", "mean_test_score", "std_test_score"]]
我们可以看到使用 'rbf' 内核的估计器表现最好，其次是 'linear'。两个使用 'poly' 内核的估计器表现较差，其中使用二次多项式的估计器表现远低于所有其他模型。
通常分析到此就结束了，但故事只讲了一半。GridSearchCV 的输出并未提供模型间差异确定性的信息。我们不知道这些差异是否具有统计学意义。为了进行评估，我们需要进行统计检验。具体来说，要对比两个模型的性能，我们应该统计比较它们的 AUC 分数。由于我们进行了 10 次 10 折交叉验证，每个模型有 100 个样本（AUC 分数）。
然而，模型的分数并不是独立的：所有模型都是在相同的 100 个分区上进行评估的，这增加了模型性能之间的相关性。由于数据的某些分区可能使所有模型都特别容易或难以区分类别，模型的分数会产生协变。
让我们通过绘制所有模型在每个折叠中的性能，并计算各折叠间模型的相关性来检查这种分区效应。
# create df of model scores ordered by performance
model_scores = results_df.filter(regex=r"split\d*_test_score")
# plot 30 examples of dependency between cv fold and AUC scores
fig, ax = plt.subplots()
sns.lineplot(
data=model_scores.transpose().iloc[:30],
dashes=False,
palette="Set1",
marker="o",
alpha=0.5,
ax=ax,
)
ax.set_xlabel("CV test fold", size=12, labelpad=10)
ax.set_ylabel("Model AUC", size=12)
ax.tick_params(bottom=True, labelbottom=False)
plt.show()
# print correlation of AUC scores across folds
print(f"Correlation of models:\n {model_scores.transpose().corr()}")
Correlation of models:
kernel rbf linear 3_poly 2_poly
kernel
rbf 1.000000 0.882561 0.783392 0.351390
linear 0.882561 1.000000 0.746492 0.298688
3_poly 0.783392 0.746492 1.000000 0.355440
2_poly 0.351390 0.298688 0.355440 1.000000
我们可以观察到模型的性能高度依赖于折叠。
因此，如果我们假设样本之间是独立的，我们将低估统计检验中计算的方差，从而增加假阳性错误的数量（即在不存在显著差异时检测出显著差异）[1]。
针对这些情况，已经开发出几种方差校正统计检验方法。在本示例中，我们将展示如何在两种不同的统计框架下实现其中一种（所谓的 Nadeau 和 Bengio 校正 t 检验）：频率派和贝叶斯派。
比较两个模型：频率派方法#
我们可以首先问：“第一个模型（按 mean_test_score 排名）是否明显优于第二个模型？”
要使用频率派方法回答这个问题，我们可以运行配对 t 检验并计算 p 值。这在预测文献中也被称为 Diebold-Mariano 检验 [5]。为了解决上一节描述的“样本非独立性问题”，已经开发了许多此类 t 检验的变体。我们将使用被证明能获得最高可重复性分数（衡量模型在评估同一数据集的不同随机分区时性能相似度）同时保持低假阳性和假阴性率的方法：使用 10 次重复 10 折交叉验证的 Nadeau 和 Bengio 校正 t 检验 [2] [3]。
此校正配对 t 检验计算公式为
其中 \(k\) 是折数，\(r\) 是交叉验证中的重复次数，\(x\) 是模型性能的差异，\(n_{test}\) 是用于测试的样本数，\(n_{train}\) 是用于训练的样本数，\(\hat{\sigma}^2\) 代表观察到的差异的方差。
让我们实现一个校正的右尾配对 t 检验，以评估第一个模型的性能是否显著优于第二个模型。我们的零假设是第二个模型的表现至少与第一个模型一样好。
import numpy as np
from scipy.stats import t
def corrected_std(differences, n_train, n_test):
"""Corrects standard deviation using Nadeau and Bengio's approach.
Parameters
----------
differences : ndarray of shape (n_samples,)
Vector containing the differences in the score metrics of two models.
n_train : int
Number of samples in the training set.
n_test : int
Number of samples in the testing set.
Returns
-------
corrected_std : float
Variance-corrected standard deviation of the set of differences.
"""
# kr = k times r, r times repeated k-fold crossvalidation,
# kr equals the number of times the model was evaluated
kr = len(differences)
corrected_var = np.var(differences, ddof=1) * (1 / kr + n_test / n_train)
corrected_std = np.sqrt(corrected_var)
return corrected_std
def compute_corrected_ttest(differences, df, n_train, n_test):
"""Computes right-tailed paired t-test with corrected variance.
Parameters
----------
differences : array-like of shape (n_samples,)
Vector containing the differences in the score metrics of two models.
df : int
Degrees of freedom.
n_train : int
Number of samples in the training set.
n_test : int
Number of samples in the testing set.
Returns
-------
t_stat : float
Variance-corrected t-statistic.
p_val : float
Variance-corrected p-value.
"""
mean = np.mean(differences)
std = corrected_std(differences, n_train, n_test)
t_stat = mean / std
p_val = t.sf(np.abs(t_stat), df) # right-tailed t-test
return t_stat, p_val
model_1_scores = model_scores.iloc[0].values # scores of the best model
model_2_scores = model_scores.iloc[1].values # scores of the second-best model
differences = model_1_scores - model_2_scores
n = differences.shape[0] # number of test sets
df = n - 1
n_train = len(next(iter(cv.split(X, y)))[0])
n_test = len(next(iter(cv.split(X, y)))[1])
t_stat, p_val = compute_corrected_ttest(differences, df, n_train, n_test)
print(f"Corrected t-value: {t_stat:.3f}\nCorrected p-value: {p_val:.3f}")
Corrected t-value: 0.750
Corrected p-value: 0.227
我们可以将校正后的 t 值和 p 值与未校正的值进行比较。
Uncorrected t-value: 2.611
Uncorrected p-value: 0.005
使用常规的 p=0.05 显著性 alpha 水平，我们观察到未校正的 t 检验得出结论认为第一个模型显著优于第二个模型。
相比之下，使用校正方法，我们无法检测到这种差异。
然而在后者情况下，频率派方法无法让我们得出第一个和第二个模型性能相等的结论。如果我们要做出这种断言，就需要使用贝叶斯方法。
比较两个模型：贝叶斯方法#
我们可以使用贝叶斯估计来计算第一个模型优于第二个模型的概率。贝叶斯估计将输出一个分布，该分布跟随两个模型性能差异的均值 \(\mu\)。
为了获得后验分布，我们需要定义一个先验（建模我们在查看数据前对均值分布的信念），并将其与似然函数相乘（该函数计算在给定均值可能的取值下，观察到的差异有多大可能）。
贝叶斯估计可以通过多种形式来回答我们的问题，但在本示例中，我们将实现 Benavoli 及其同事建议的方法 [4]。
使用闭式表达式定义后验的一种方法是选择与似然函数共轭的先验。Benavoli 及其同事 [4] 表明，在比较两个分类器的性能时，我们可以将先验建模为正态-伽马分布（均值和方差均未知），该分布与正态似然函数共轭，从而将后验表示为正态分布。从该正态后验中边缘化方差，我们可以将均值参数的后验定义为学生 t 分布。具体为
其中 \(n\) 是样本总数，\(\overline{x}\) 代表分数的平均差异，\(n_{test}\) 是用于测试的样本数，\(n_{train}\) 是用于训练的样本数，\(\hat{\sigma}^2\) 代表观察到的差异的方差。
注意，我们在贝叶斯方法中也使用了 Nadeau 和 Bengio 的校正方差。
让我们计算并绘制后验分布。
让我们绘制后验分布。
x = np.linspace(t_post.ppf(0.001), t_post.ppf(0.999), 100)
plt.plot(x, t_post.pdf(x))
plt.xticks(np.arange(-0.04, 0.06, 0.01))
plt.fill_between(x, t_post.pdf(x), 0, facecolor="blue", alpha=0.2)
plt.ylabel("Probability density")
plt.xlabel(r"Mean difference ($\mu$)")
plt.title("Posterior distribution")
plt.show()
我们可以通过计算后验分布曲线从零到无穷大的面积来计算第一个模型优于第二个模型的概率。反之亦然：我们可以通过计算从负无穷大到零的曲线面积来计算第二个模型优于第一个模型的概率。
better_prob = 1 - t_post.cdf(0)
print(
f"Probability of {model_scores.index[0]} being more accurate than "
f"{model_scores.index[1]}: {better_prob:.3f}"
)
print(
f"Probability of {model_scores.index[1]} being more accurate than "
f"{model_scores.index[0]}: {1 - better_prob:.3f}"
)
Probability of rbf being more accurate than linear: 0.773
Probability of linear being more accurate than rbf: 0.227
与频率派方法相反，我们可以计算一个模型优于另一个模型的概率。
请注意，我们得到的结果与频率派方法中的结果相似。鉴于我们对先验的选择，我们本质上是在进行相同的计算，但可以做出不同的断言。
实际等效区域 (ROPE)#
有时我们有兴趣确定模型性能相等的概率，其中“相等”是在实际意义上定义的。一种简单的做法 [4] 是将估计器定义为在准确率差异小于 1% 时具有实际等效性。但我们也可以结合所解决的问题来定义这种实际等效性。例如，准确率 5% 的差异可能意味着销售额增加 1000 美元，我们认为任何高于此的量对我们的业务都有意义。
在本示例中，我们将实际等效区域 (ROPE) 定义为 \([-0.01, 0.01]\)。也就是说，如果两个模型的性能差异小于 1%，我们将认为它们在实际应用中是等效的。
为了计算分类器在实际中等效的概率，我们计算 ROPE 区间内后验曲线下的面积。
rope_interval = [-0.01, 0.01]
rope_prob = t_post.cdf(rope_interval[1]) - t_post.cdf(rope_interval[0])
print(
f"Probability of {model_scores.index[0]} and {model_scores.index[1]} "
f"being practically equivalent: {rope_prob:.3f}"
)
Probability of rbf and linear being practically equivalent: 0.432
我们可以绘制后验分布在 ROPE 区间上的分布情况。
x_rope = np.linspace(rope_interval[0], rope_interval[1], 100)
plt.plot(x, t_post.pdf(x))
plt.xticks(np.arange(-0.04, 0.06, 0.01))
plt.vlines([-0.01, 0.01], ymin=0, ymax=(np.max(t_post.pdf(x)) + 1))
plt.fill_between(x_rope, t_post.pdf(x_rope), 0, facecolor="blue", alpha=0.2)
plt.ylabel("Probability density")
plt.xlabel(r"Mean difference ($\mu$)")
plt.title("Posterior distribution under the ROPE")
plt.show()
正如 [4] 中所建议的那样，我们可以使用与频率派方法相同的准则进一步解释这些概率：落在 ROPE 内的概率是否大于 95%（5% 的 alpha 值）？如果是，我们可以得出结论，两个模型在实际应用中是等效的。
贝叶斯估计方法还允许我们计算对差异估计的不确定性程度。这可以使用可信区间来计算。对于给定的概率，它们显示了估计量（在我们的例子中是性能的平均差异）可以取的数值范围。例如，50% 的可信区间 [x, y] 告诉我们，模型间性能的真实（平均）差异落在 x 和 y 之间的概率为 50%。
让我们使用 50%、75% 和 95% 的置信度确定我们数据的可信区间。
cred_intervals = []
intervals = [0.5, 0.75, 0.95]
for interval in intervals:
cred_interval = list(t_post.interval(interval))
cred_intervals.append([interval, cred_interval[0], cred_interval[1]])
cred_int_df = pd.DataFrame(
cred_intervals, columns=["interval", "lower value", "upper value"]
).set_index("interval")
cred_int_df
如表所示，模型间真实平均差异在 0.000977 和 0.019023 之间的概率为 50%，在 -0.005422 和 0.025422 之间的概率为 70%，在 -0.016445 和 0.036445 之间的概率为 95%。
所有模型的成对比较：频率派方法#
我们也可能对比较所有使用 GridSearchCV 评估的模型性能感兴趣。在这种情况下，我们将多次运行统计检验，这会导致 多重比较问题。
解决这个问题有很多可能的方法，但标准方法是应用 Bonferroni 校正。Bonferroni 可以通过将 p 值乘以我们测试的比较次数来计算。
让我们使用校正后的 t 检验比较模型的性能。
from itertools import combinations
from math import factorial
n_comparisons = factorial(len(model_scores)) / (
factorial(2) * factorial(len(model_scores) - 2)
)
pairwise_t_test = []
for model_i, model_k in combinations(range(len(model_scores)), 2):
model_i_scores = model_scores.iloc[model_i].values
model_k_scores = model_scores.iloc[model_k].values
differences = model_i_scores - model_k_scores
t_stat, p_val = compute_corrected_ttest(differences, df, n_train, n_test)
p_val *= n_comparisons # implement Bonferroni correction
# Bonferroni can output p-values higher than 1
p_val = 1 if p_val > 1 else p_val
pairwise_t_test.append(
[model_scores.index[model_i], model_scores.index[model_k], t_stat, p_val]
)
pairwise_comp_df = pd.DataFrame(
pairwise_t_test, columns=["model_1", "model_2", "t_stat", "p_val"]
).round(3)
pairwise_comp_df
我们观察到，在校正了多重比较后，唯一与其他模型显著不同的模型是 '2_poly'。'rbf'（被 GridSearchCV 排名第一的模型）与 'linear' 或 '3_poly' 没有显著差异。
所有模型的成对比较：贝叶斯方法#
当使用贝叶斯估计来比较多个模型时，我们不需要进行多重比较校正（原因请参阅 [4]）。
我们可以像第一节那样进行成对比较。
pairwise_bayesian = []
for model_i, model_k in combinations(range(len(model_scores)), 2):
model_i_scores = model_scores.iloc[model_i].values
model_k_scores = model_scores.iloc[model_k].values
differences = model_i_scores - model_k_scores
t_post = t(
df, loc=np.mean(differences), scale=corrected_std(differences, n_train, n_test)
)
worse_prob = t_post.cdf(rope_interval[0])
better_prob = 1 - t_post.cdf(rope_interval[1])
rope_prob = t_post.cdf(rope_interval[1]) - t_post.cdf(rope_interval[0])
pairwise_bayesian.append([worse_prob, better_prob, rope_prob])
pairwise_bayesian_df = pd.DataFrame(
pairwise_bayesian, columns=["worse_prob", "better_prob", "rope_prob"]
).round(3)
pairwise_comp_df = pairwise_comp_df.join(pairwise_bayesian_df)
pairwise_comp_df
使用贝叶斯方法，我们可以计算一个模型表现得比另一个模型更好、更差或实际等效的概率。
结果表明，被 GridSearchCV 排名第一的 'rbf' 模型，比 'linear' 更差的概率约为 6.8%，比 '3_poly' 更差的概率约为 1.8%。'rbf' 和 'linear' 有 43% 的概率在实际应用中等效，而 'rbf' 和 '3_poly' 有 10% 的概率如此。
与频率派方法得出的结论类似，所有模型都有 100% 的概率比 '2_poly' 表现更好，且没有模型与后者具有实际等效的性能。
核心要点#
性能指标上的细微差异很容易仅仅是偶然造成的，而不是因为一个模型系统地预测得比另一个好。正如本示例所示，统计学可以告诉你这种情况的可能性。
在统计比较 GridSearchCV 评估的两个模型性能时，有必要校正计算出的方差，因为模型的分数之间并非独立，可能导致方差被低估。
使用（方差校正）配对 t 检验的频率派方法可以告诉我们一个模型的性能是否以高于随机性的确定程度优于另一个模型。
贝叶斯方法可以提供一个模型比另一个模型更好、更差或在实际中等效的概率。它还可以告诉我们对模型差异的真实取值落在某个数值范围内的确信程度。
如果要统计比较多个模型，在使用频率派方法时需要进行多重比较校正。
References
脚本运行总耗时： (0 分钟 2.407 秒)
相关示例
