注意
转到末尾以下载完整示例代码或通过 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 个分区上进行评估，这增加了模型性能之间的相关性。由于某些数据分区可能使得所有模型都特别容易或难以找到类别区分，因此模型得分会共同变化。
让我们通过绘制每个折叠中所有模型的性能，并计算模型在折叠之间的相关性来检查这种分区效应
# 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 检验的许多变体。我们将使用被证明能获得最高可复制性分数（评估模型在同一数据集的不同随机分区上评估时的性能相似程度）同时保持较低假阳性率和假阴性率的变体：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(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
使用常规的显著性 alpha 水平 p=0.05，我们观察到未校正的 t 检验得出结论：第一个模型显着优于第二个模型。
相比之下，使用校正方法，我们未能检测到这种差异。
然而，在后一种情况下，频率派方法不允许我们得出第一个和第二个模型具有相同性能的结论。如果我们要做出这种断言，我们需要使用贝叶斯方法。
比较两个模型：贝叶斯方法#
我们可以使用贝叶斯估计来计算第一个模型优于第二个模型的概率。贝叶斯估计将输出一个分布，该分布的均值 \(\mu\) 是两个模型性能差异的均值。
为了获得后验分布，我们需要定义一个先验分布来建模我们在查看数据之前对均值分布的信念，并将其乘以一个似然函数，该函数计算给定均值可能取的值时，我们观察到的差异有多大可能性。
贝叶斯估计可以用多种形式来回答我们的问题，但在本例中，我们将实现 Benavoli 及其同事建议的方法 [4]。
使用闭式表达式定义后验的一种方法是选择一个与似然函数共轭的先验分布。Benavoli 及其同事 [4] 表明，在比较两个分类器的性能时，我们可以将先验建模为 Normal-Gamma 分布（均值和方差均未知），该分布与正态似然共轭，从而将后验表示为正态分布。通过对这个正态后验的方差进行边缘化，我们可以将均值参数的后验定义为 Student’s 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'。 GridSearchCV 排名第一的模型 'rbf' 与 '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'，有大约 6.8% 的可能性比 'linear' 差，有 1.8% 的可能性比 '3_poly' 差。 'rbf' 和 'linear' 有 43% 的概率实际等效，而 'rbf' 和 '3_poly' 有 10% 的概率实际等效。
与使用频率派方法得出的结论相似，所有模型都有 100% 的概率优于 '2_poly'，并且没有一个模型与后者具有实际等效的性能。
总结要点#
性能度量的微小差异很容易仅仅是偶然造成的，而不是因为一个模型系统性地优于另一个模型。如本例所示，统计数据可以告诉你这种可能性的高低。
当对在 GridSearchCV 中评估的两个模型的性能进行统计比较时，有必要校正计算出的方差，因为模型的得分不是相互独立的，这可能导致方差被低估。
使用（方差校正的）配对 t 检验的频率派方法可以告诉我们一个模型的性能是否以高于偶然的确定性优于另一个模型。
贝叶斯方法可以提供一个模型比另一个模型更好、更差或实际等效的概率。它还可以告诉我们对模型真实差异落在特定值范围内的确定程度。
如果对多个模型进行统计比较，使用频率派方法时需要进行多重比较校正。
References
脚本总运行时间： (0 minutes 1.443 seconds)
相关示例
