注意
转到结尾 下载完整的示例代码。或通过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()
我们将比较SVC估计器的性能，这些估计器在其kernel参数上有所不同，以确定此超参数的哪个选择最能预测我们的模拟数据。我们将使用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个分区上进行评估，这增加了模型性能之间的相关性。由于某些数据分区可能使所有模型都特别容易或难以区分类别，因此模型分数将共同变化。
让我们通过绘制每个fold中所有模型的性能，并计算模型在folds之间的相关性来检查这种分区效应。
# 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
我们可以观察到，模型的性能高度依赖于fold。
因此，如果我们假设样本之间是独立的，我们将低估在统计检验中计算的方差，从而增加假阳性错误的数量（即，在不存在这种差异时检测到模型之间存在显著差异）[1]。
已经为这些情况开发了几种方差校正统计检验。在本例中，我们将展示如何在两种不同的统计框架下实现其中一种（所谓的Nadeau和Bengio校正t检验）：频率论和贝叶斯。
比较两个模型：频率论方法#
我们可以从以下问题开始：“根据mean_test_score排名，第一个模型是否显著优于第二个模型？”
为了使用频率论方法回答这个问题，我们可以运行配对t检验并计算p值。这在预测文献中也称为Diebold-Mariano检验[5]。已经开发了许多此类t检验的变体来解决上一节中描述的“样本非独立性问题”。我们将使用经证明具有最高可重复性分数（评估相同数据集的不同随机分区时模型的性能有多相似）同时保持低假阳性和假阴性率的一种方法：Nadeau和Bengio校正t检验[2]，它使用10次重复的10倍交叉验证[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(list(cv.split(X, y))[0][0])
n_test = len(list(cv.split(X, y))[0][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，我们观察到未修正的 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
与频率论方法相比，我们可以计算一个模型优于另一个模型的概率。
请注意，我们获得的结果与频率论方法中的结果相似。鉴于我们选择的先验，我们实际上执行的是相同的计算，但我们可以做出不同的断言。
实际等效区域#
有时我们感兴趣的是确定我们的模型具有等效性能的概率，其中“等效”以实际的方式定义。一种简单的方法 [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%（alpha 值为 5%）？在这种情况下，我们可以得出结论，这两个模型在实践中是等效的。
贝叶斯估计方法还允许我们计算我们对差异估计的不确定性。这可以使用可信区间来计算。对于给定的概率，它们显示了估计量（在本例中为性能差异的平均值）可以取的值的范围。例如，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%。
与使用频率论方法获得的结论类似，所有模型比'2_poly'更好的概率为100%，并且没有模型与后者具有实际等效的性能。
要点#
性能指标的微小差异很可能仅仅是偶然的，而不是因为一个模型比另一个模型系统地预测得更好。如本例所示，统计数据可以告诉你这种情况发生的可能性。
当统计比较在GridSearchCV中评估的两个模型的性能时，需要校正计算出的方差，因为模型的得分彼此不独立，可能会低估方差。
使用（方差校正）配对t检验的频率论方法可以告诉我们，一个模型的性能比另一个模型更好的确定性程度是否高于偶然性。
贝叶斯方法可以提供一个模型比另一个模型更好、更差或实际上等效的概率。它还可以告诉我们，我们对知道模型的真实差异落在某个特定值范围内有多大的信心。
如果统计比较多个模型，则在使用频率论方法时需要进行多重比较校正。
参考文献
脚本总运行时间：（0分钟1.644秒）
相关示例
