实现无需定制软件包的分类拟合预测
本博客文章受Chris Molner的书《使用Python进行符合预测介绍》的启发。 Chris在使新的机器学习技术对他人可及的方面非常出色。 我特别推荐他关于可解释机器学习的书籍。
完整代码的GitHub存储库可以在此处找到：符合预测。
什么是符合预测？
符合预测既是一种不确定性量化方法，也是一种分类实例的方法（可以对类别或子群进行微调）。 通过将分类结果转换为潜在类别的集合，来传达不确定性。
符合预测指定了覆盖度，即真实结果被预测区域覆盖的概率。 符合预测中预测区域的解释取决于任务。 对于分类，我们得到预测集，而对于回归，我们得到预测区间。
以下是“传统”分类（基于最可能类别）和符合预测（集合）之间的区别示例。
此方法的优点包括：
保证覆盖度：符合预测生成的预测集具有对真实结果的覆盖度保证，即它们将检测出您设定的最低目标覆盖度下的任何真实值。 符合预测不依赖于经过良好校准的模型-唯一重要的是，与所有机器学习一样，进行分类的新样本必须来自于与训练和校准数据相似的数据分布。 覆盖度还可以在类别或子群之间得到保证，尽管这需要在我们将要介绍的方法中增加一步。
- 易于使用：符合预测方法可以从头开始实现，并只需几行代码，就像我们在这里做的那样。
- 模型无关：符合预测适用于任何机器学习模型。 它使用您所选择的任何模型的正常输出。
- 无分布要求：符合预测不对数据的基础分布做出任何假设；它是一种非参数方法。
- 无需重新训练：符合预测可以在不重新训练模型的情况下使用。 这是另一种观察和使用模型输出的方式。
- 广泛应用：符合预测适用于表格数据分类、图像或时间序列分类、回归以及许多其他任务，尽管我们仅在这里演示分类。
为什么我们应该关心不确定性量化？
在许多情况下，不确定性量化至关重要：
- 当我们使用模型预测进行决策时。 我们对这些预测有多有把握？ 使用“最可能的类别”对于我们的任务是否足够好？
- 当我们想要将与预测相关的不确定性传达给利益相关方时，而不谈论概率、赔率甚至对数赔率！
符合预测中的alpha – 描述覆盖度
覆盖度是符合预测的关键。 在分类中，它是特定类别所属的数据正常范围。 覆盖度等同于敏感性或召回率；它是被分类集合中识别的观察值比例。 通过调整𝛼（覆盖度= 1-𝛼），我们可以缩小或放宽覆盖面积。
让我们编码吧！
导入软件包
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from sklearn.datasets import make_blobs
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split创建分类的合成数据
将使用SK-Learn的`make_blobs`方法生成示例数据。
n_classes = 3
# 创建训练和测试数据
X, y = make_blobs(n_samples=10000, n_features=2, centers=n_classes, cluster_std=3.75, random_state=42)
# 减小第一类的大小以创建不平衡数据集
# 设置numpy随机种子
np.random.seed(42)
# 获取y为类别0时的索引
class_0_idx = np.where(y == 0)[0]
# 获取0类索引的30%样本
class_0_idx = np.random.choice(class_0_idx, int(len(class_0_idx) * 0.3), replace=False)
# 获取其它类的索引
rest_idx = np.where(y != 0)[0]
# 合并索引
idx = np.concatenate([class_0_idx, rest_idx])
# 打乱索引
np.random.shuffle(idx)
# 拆分数据
X = X[idx]
y = y[idx]
# 拆分出用于模型训练的数据集
X_train, X_rest, y_train, y_rest = train_test_split(X, y, test_size=0.5, random_state=42)
# 将剩余数据拆分为校准集和测试集
X_Cal, X_test, y_cal, y_test = train_test_split(X_rest, y_rest, test_size=0.5, random_state=42)
# 设置类别标签
class_labels = ['blue', 'orange', 'green']
# 绘制数据
fig = plt.subplots(figsize=(5, 5))
ax = plt.subplot(111)
for i in range(n_classes):
ax.scatter(X_test[y_test == i, 0], X_test[y_test == i, 1], label=class_labels[i], alpha=0.5, s=10)
legend = ax.legend()
legend.set_title("类别")
ax.set_xlabel("特征1")
ax.set_ylabel("特征2")
plt.show()
构建分类器
我们将在这里使用一个简单的逻辑回归模型，但该方法适用于基于表格数据的任何模型，从简单的基于表格数据的逻辑回归模型到用于图像分类的3D ConvNets。
# 构建并训练分类器
classifier = LogisticRegression(random_state=42)
classifier.fit(X_train, y_train)
# 测试分类器
y_pred = classifier.predict(X_test)
accuracy = np.mean(y_pred == y_test)
print(f"准确率: {accuracy:0.3f}")
# 测试每个类别的召回率
for i in range(n_classes):
recall = np.mean(y_pred[y_test == i] == y_test[y_test == i])
print(f"类别 {class_labels[i]} 的召回率: {recall:0.3f}")
准确率: 0.930
类别 blue 的召回率: 0.772
类别 orange 的召回率: 0.938
类别 green 的召回率: 0.969注意，少数类别的召回率低于其他类别。召回率，也称为敏感性，是分类器正确识别出某类别的数量。
S_i ，或 非符合分数 分数
在符合性预测中，非符合分数通常用 s_i 表示，它是一个衡量新实例与训练集中现有实例的差异程度的指标。它用于确定新实例是否属于特定类别。
在分类的上下文中，最常见的非符合度量是 1 减去给定标签的预测类别概率。所以，如果新实例属于某类别的预测概率较高，则非符合分数将较低，反之亦然。
对于一致预测，我们获得所有类别的s_i分数（注意：我们只看实例的真实类别的模型输出，即使它在其他类别中的预测概率更高）。然后，我们找到一个包含（或覆盖）95％数据的分数阈值。分类将会识别出95％新的实例（只要我们的新数据与训练数据相似）。
计算一致预测阈值
现在我们将预测校准集的分类概率。这将用于为新数据设置分类阈值。
# 获取校准集的预测y_pred = classifier.predict(X_Cal)y_pred_proba = classifier.predict_proba(X_Cal)# 显示前5个实例y_pred_proba[0:5]
array([[4.65677826e-04, 1.29602253e-03, 9.98238300e-01], [1.73428257e-03, 1.20718182e-02, 9.86193899e-01], [2.51649788e-01, 7.48331668e-01, 1.85434981e-05], [5.97545130e-04, 3.51642214e-04, 9.99050813e-01], [4.54193815e-06, 9.99983628e-01, 1.18300819e-05]])计算非一致性得分：
我们只根据与观察类别相关的概率来计算s_i分数。对于每个实例，我们会获取该实例类别的预测概率。s_i分数（非一致性）为1减去概率。s_i分数越高，该实例与其他类别相比，越不符合该类别。
si_scores = []# 循环遍历所有校准实例for i, true_class in enumerate(y_cal): # 获取观察/真实类别的预测概率 predicted_prob = y_pred_proba[i][true_class] si_scores.append(1 - predicted_prob) # 转换为 NumPy 数组si_scores = np.array(si_scores)# 显示前5个实例si_scores[0:5]
array([1.76170035e-03, 1.38061008e-02, 2.51668332e-01, 9.49187344e-04, 1.63720201e-05])获取95th百分位阈值：
阈值确定了我们的分类的覆盖范围。覆盖率指的是实际包含真实结果的预测比例。
阈值是与1-𝛼相对应的百分位数。要获得95％的覆盖率，我们设置𝛼为0.05。
在实际应用中，基于𝛼的分位水平需要进行有限样本修正来计算相应的分位数𝑞。我们将0.95乘以$(n+1)/n$，这意味着当n = 1000时，𝑞𝑙𝑒𝑣𝑒𝑙将为0.951。
number_of_samples = len(X_Cal)alpha = 0.05qlevel = (1 - alpha) * ((number_of_samples + 1) / number_of_samples)threshold = np.percentile(si_scores, qlevel*100)print(f'阈值：{threshold:0.3f}')
阈值：0.598显示s_i值的图表，带有截断阈值。
x = np.arange(len(si_scores)) + 1sorted_si_scores = np.sort(si_scores)index_of_95th_percentile = int(len(si_scores) * 0.95)# 通过截断阈值上色conform = 'g' * index_of_95th_percentilenonconform = 'r' * (len(si_scores) - index_of_95th_percentile)color = list(conform + nonconform)fig = plt.figure(figsize=((6,4)))ax = fig.add_subplot()# 添加条形图ax.bar(x, sorted_si_scores, width=1.0, color = color)# 添加95th百分位线ax.plot([0, index_of_95th_percentile],[threshold, threshold], c='k', linestyle='--')ax.plot([index_of_95th_percentile, index_of_95th_percentile], [threshold, 0], c='k', linestyle='--')# 添加文本txt = '95th百分位一致性阈值'ax.text(5, threshold + 0.04, txt)# 添加坐标轴标签ax.set_xlabel('样本实例（按$s_i$排序）')ax.set_ylabel('$S_i$（非一致性）')plt.show()
获取被分类为阳性的测试集样本/类别
现在可以找出所有模型输出小于阈值的样本/类别。
一个个体示例可能没有值，或者大于一个值小于阈值。
prediction_sets = (1 - classifier.predict_proba(X_test) <= threshold)# 显示前十个实例prediction_sets[0:10]
array([[ True, False, False], [False, False, True], [ True, False, False], [False, False, True], [False, True, False], [False, True, False], [False, True, False], [ True, True, False], [False, True, False], [False, True, False]])获取预测集标签，并与标准分类进行比较。
# 获取标准预测结果y_pred = classifier.predict(X_test)# 函数用于获取预测集标签def get_prediction_set_labels(prediction_set, class_labels): # 为预测集中的每个实例获取类别标签集合 prediction_set_labels = [ set([class_labels[i] for i, x in enumerate(prediction_set) if x]) for prediction_set in prediction_sets] return prediction_set_labels# 整理预测结果results_sets = pd.DataFrame()results_sets['观测'] = [class_labels[i] for i in y_test]results_sets['标签'] = get_prediction_set_labels(prediction_sets, class_labels)results_sets['分类'] = [class_labels[i] for i in y_pred]results_sets.head(10)
观测 标签 分类0 blue {blue} blue1 green {green} green2 blue {blue} blue3 green {green} green4 orange {orange} orange5 orange {orange} orange6 orange {orange} orange7 orange {blue, orange} blue8 orange {orange} orange9 orange {orange} orange注意实例7实际上是橙色类别，但简单分类器将其分类为蓝色。遵循性预测将其分类为橙色和蓝色的集合。
绘制数据，显示被预测为可能属于2个类别的实例7：
# 绘制数据fig = plt.subplots(figsize=(5, 5))ax = plt.subplot(111)for i in range(n_classes): ax.scatter(X_test[y_test == i, 0], X_test[y_test == i, 1], label=class_labels[i], alpha=0.5, s=10)# 添加实例7set_label = results_sets['标签'].iloc[7]ax.scatter(X_test[7, 0], X_test[7, 1], color='k', s=100, marker='*', label=f'实例7')legend = ax.legend()legend.set_title("类别")ax.set_xlabel("特征1")ax.set_ylabel("特征2")txt = f"实例7的预测集：{set_label}"ax.text(-20, 18, txt)plt.show()
显示覆盖率和平均集合大小
覆盖率是实际包含真实结果的预测集比例。
平均集合大小是每个实例中预测的类别的平均数。
我们将定义一些函数来计算结果。
# 获取每个类别的计数def get_class_counts(y_test): class_counts = [] for i in range(n_classes): class_counts.append(np.sum(y_test == i)) return class_counts# 获取每个类别的覆盖率def get_coverage_by_class(prediction_sets, y_test): coverage = [] for i in range(n_classes): coverage.append(np.mean(prediction_sets[y_test == i, i])) return coverage# 获取每个类别的平均集合大小def get_average_set_size(prediction_sets, y_test): average_set_size = [] for i in range(n_classes): average_set_size.append( np.mean(np.sum(prediction_sets[y_test == i], axis=1))) return average_set_size # 获取加权覆盖率（按类别大小加权）def get_weighted_coverage(coverage, class_counts): total_counts = np.sum(class_counts) weighted_coverage = np.sum((coverage * class_counts) / total_counts) weighted_coverage = round(weighted_coverage, 3) return weighted_coverage# 获取加权集合大小（按类别大小加权）def get_weighted_set_size(set_size, class_counts): total_counts = np.sum(class_counts) weighted_set_size = np.sum((set_size * class_counts) / total_counts) weighted_set_size = round(weighted_set_size, 3) return weighted_set_size显示每个类别的结果。
results = pd.DataFrame(index=class_labels)results['类别计数'] = get_class_counts(y_test)results['覆盖率'] = get_coverage_by_class(prediction_sets, y_test)results['平均集合大小'] = get_average_set_size(prediction_sets, y_test)results
类别计数 覆盖率 平均集合大小blue 241 0.817427 1.087137orange 848 0.954009 1.037736green 828 0.977053 1.016908显示整体结果。
weighted_coverage = get_weighted_coverage( results['覆盖率'], results['类别计数'])weighted_set_size = get_weighted_set_size( results['平均集合大小'], results['类别计数'])print (f'整体覆盖率：{weighted_coverage}')print (f'平均集合大小：{weighted_set_size}')
整体覆盖率：0.947平均集合大小：1.035注意：虽然我们的整体覆盖率符合要求，非常接近95％，但不同类别的覆盖率有所差异，并且对于最小的类别，覆盖率最低（83％）。如果个别类别的覆盖率很重要，我们可以为类别单独设定阈值，这就是我们现在要做的。
具有相等覆盖率的一致分类
当我们想要确保对所有类别的覆盖率时，我们可以为每个类别单独设置阈值。
注意：如果使用共享阈值导致问题，我们也可以对数据的子组（例如，确保对不同种族群体的诊断具有相等覆盖率）进行此操作。
为每个类别单独获取阈值
# 设置 alpha（1 - 覆盖率）alpha = 0.05thresholds = []# 获取校准集的预测概率y_cal_prob = classifier.predict_proba(X_Cal)# 获取每个类别的第95个百分位数s-scoresfor class_label in range(n_classes): mask = y_cal == class_label y_cal_prob_class = y_cal_prob[mask][:, class_label] s_scores = 1 - y_cal_prob_class q = (1 - alpha) * 100 class_size = mask.sum() correction = (class_size + 1) / class_size q *= correction threshold = np.percentile(s_scores, q) thresholds.append(threshold)print(thresholds)
[0.9030202125697161, 0.6317149025299887, 0.26033562285411]对每个类别的分类应用特定的阈值
# 获取测试集的 Si 分数predicted_proba = classifier.predict_proba(X_test)si_scores = 1 - predicted_proba# 对于每个类别，检查每个实例是否低于阈值prediction_sets = []for i in range(n_classes): prediction_sets.append(si_scores[:, i] <= thresholds[i])prediction_sets = np.array(prediction_sets).T# 获取预测集标签并显示前10个prediction_set_labels = get_prediction_set_labels(prediction_sets, class_labels)# 获取标准预测y_pred = classifier.predict(X_test)# 收集预测结果results_sets = pd.DataFrame()results_sets['观测值'] = [class_labels[i] for i in y_test]results_sets['标签'] = get_prediction_set_labels(prediction_sets, class_labels)results_sets['分类'] = [class_labels[i] for i in y_pred]# 显示前10个结果results_sets.head(10)
观测值 标签 分类0 蓝色 {蓝色} 蓝色1 绿色 {绿色} 绿色2 蓝色 {蓝色} 蓝色3 绿色 {绿色} 绿色4 橙色 {橙色} 橙色5 橙色 {橙色} 橙色6 橙色 {橙色} 橙色7 橙色 {蓝色, 橙色} 蓝色8 橙色 {橙色} 橙色9 橙色 {橙色} 橙色检查覆盖率并在类别间设置大小
我们现在在所有类别中有大约95%的覆盖率。与标准的分类方法相比，适合性预测方法给我们的少数类别提供了更好的覆盖率。
results = pd.DataFrame(index=class_labels)results['Class counts'] = get_class_counts(y_test)results['Coverage'] = get_coverage_by_class(prediction_sets, y_test)results['Average set size'] = get_average_set_size(prediction_sets, y_test)results
类别数量 覆盖率 平均集合大小
蓝色 241 0.954357 1.228216
橙色 848 0.956368 1.139151
绿色 828 0.942029 1.006039
加权覆盖率 = get_weighted_coverage( results['Coverage'], results['Class counts'])加权集合大小 = get_weighted_set_size( results['Average set size'], results['Class counts'])print (f'整体覆盖率：{weighted_coverage}')print (f'平均集合大小：{weighted_set_size}')
整体覆盖率：0.95平均集合大小：1.093总结
适合性预测被用来将实例分类为集合而不是单个预测结果。位于两个类别之间的实例被标记为两个类别，而不是选择概率最高的类别。
当重要的是所有类别都以相同的覆盖率被检测到时，可以单独设置实例分类的阈值（这种方法也可以用于数据的子群，例如为不同族群保证相同的覆盖率）。
适合性预测不会改变模型的预测结果。它只是以不同的方式使用这些结果进行传统分类之外的分类。它可以与传统方法并用。
（所有图片由作者提供）