Python-可解释的机器学习-三-

Python 可解释的机器学习（三）

原文：annas-archive.org/md5/9b99a159e8340372894ac9bde8bbd5d9

译者：飞龙

协议：CC BY-NC-SA 4.0

第七章：可视化卷积神经网络

到目前为止，我们只处理了表格数据，以及在第五章，局部模型无关解释方法中简要提到的文本数据。本章将专门探讨适用于图像的解释方法，特别是训练图像分类器的卷积神经网络（CNN）模型。通常，深度学习模型被视为黑盒模型的典范。然而，CNN 的一个优点是它很容易进行可视化，因此我们不仅可以可视化结果，还可以通过激活来可视化学习过程中的每一步。在所谓的黑盒模型中，解释这些步骤的可能性是罕见的。一旦我们掌握了 CNN 的学习方式，我们将研究如何使用最先进的基于梯度的属性方法，如显著性图和Grad-CAM来调试类别属性。最后，我们将通过基于扰动的属性方法，如遮挡敏感性和KernelSHAP来扩展我们的属性调试知识。

这些是我们将要讨论的主要主题：

• 使用传统解释方法评估 CNN 分类器

• 使用基于激活的方法可视化学习过程

• 使用基于梯度的属性方法评估误分类

• 使用基于扰动的属性方法理解分类

技术要求

本章的示例使用了mldatasets、pandas、numpy、sklearn、tqdm、torch、torchvision、pytorch-lightning、efficientnet-pytorch、torchinfo、matplotlib、seaborn和captum库。如何安装所有这些库的说明在前言中。

本章的代码位于此处：packt.link/qzUvD。

任务

全球每年产生超过 20 亿吨垃圾，预计到 2050 年将增长到超过 35 亿吨。近年来，全球垃圾产量急剧上升和有效废物管理系统需求日益迫切。在高收入国家，超过一半的家庭垃圾是可回收的，在低收入国家为 20%，并且还在上升。目前，大多数垃圾最终都堆放在垃圾填埋场或焚烧，导致环境污染和气候变化。考虑到全球范围内，很大一部分可回收材料没有得到回收，这是可以避免的。

假设可回收垃圾被收集，但仍可能很难且成本高昂地进行分类。以前，废物分类技术包括：

• 通过旋转圆柱形筛网（“摇筛”）按尺寸分离材料

• 通过磁力和磁场分离铁和非铁金属（“涡流分离器”）

• 通过空气按重量分离

• 通过水按密度分离（“沉浮分离”）

• 由人工执行的手动分类

即使对于大型、富裕的城市市政府，有效地实施所有这些技术也可能具有挑战性。为了应对这一挑战，智能回收系统应运而生，利用计算机视觉和人工智能高效、准确地分类废物。

智能回收系统的发展可以追溯到 2010 年代初，当时研究人员和革新者开始探索计算机视觉和人工智能改善废物管理流程的潜力。他们首先开发了基本的图像识别算法，利用颜色、形状和纹理等特征来识别废物材料。这些系统主要用于研究环境，商业应用有限。随着机器学习和人工智能的进步，智能回收系统经历了显著的改进。卷积神经网络（CNN）和其他深度学习技术使这些系统能够从大量数据中学习并提高其废物分类的准确性。此外，人工智能驱动的机器人集成使得废物材料的自动化分拣和处理成为可能，从而提高了回收工厂的效率。

摄像头、机器人和用于低延迟、高容量场景运行深度学习模型的芯片等成本与十年前相比显著降低，这使得最先进的智能回收系统对甚至更小、更贫穷的城市废物管理部门也变得可负担。巴西的一个城市正在考虑翻新他们 20 年前建成的一个由各种机器拼凑而成的回收厂，这些机器的集体分拣准确率仅为 70%。人工分拣只能部分弥补这一差距，导致不可避免的污染和污染问题。该巴西市政府希望用一条单条传送带替换当前系统，这条传送带由一系列机器人高效地将 12 种不同类别的废物分拣到垃圾桶中。

他们购买了传送带、工业机器人和摄像头。然后，他们支付了一家人工智能咨询公司开发一个用于分类可回收物的模型。然而，他们想要不同大小的模型，因为他们不确定这些模型在他们的硬件上运行的速度有多快。

如请求，咨询公司带回了 4 到 6400 万参数之间各种大小的模型。最大的模型（b7）比最小的模型（b0）慢六倍以上。然而，最大的模型在验证 F1 分数上显著更高，达到 96%（F1 val），而最小的模型大约为 90%：

图 7.1：由人工智能咨询公司提供的模型 F1 分数

市政领导对结果感到非常满意，但也感到惊讶，因为顾问们要求不要提供任何领域知识或用于训练模型的数据，这使得他们非常怀疑。他们要求回收厂的工人用一批可回收物测试这些模型。他们用这一批次的模型得到了 25%的错误分类率。

为了寻求第二意见和模型的诚实评估，市政厅联系了另一家 AI 咨询公司——你的公司！

第一项任务是组装一个更符合回收工厂工人在误分类中发现的边缘情况的测试数据集。你的同事使用测试数据集获得了 62%到 66%的 F1 分数（F1 测试）。接下来，他们要求你理解导致这些误分类的原因。

方法

没有一种解释方法完美无缺，即使是最好的情况也只能告诉你故事的一部分。因此，你决定首先使用传统的解释方法来评估模型的预测性能，包括以下方法：

• ROC 曲线和 ROC-AUC

• 混淆矩阵以及由此派生的一些指标，如准确率、精确率、召回率和 F1

然后，你将使用基于激活的方法检查模型：

• 中间激活

这之后是使用三种基于梯度的方法评估决策：

• 显著性图

• Grad-CAM

• 集成梯度

以及一个基于反向传播的方法：

• DeepLIFT

这之后是三种基于扰动的算法：

• 遮蔽敏感性

• 特征消除

• Shapley 值采样

我希望你在这一过程中理解为什么模型的表现不符合预期，以及如何修复它。你还可以利用你将生成的许多图表和可视化来向市政厅的行政人员传达这个故事。

准备工作

你会发现这个示例的大部分代码都在这里：github.com/PacktPublishing/Interpretable-Machine-Learning-with-Python-2E/blob/main/07/GarbageClassifier.ipynb

加载库

要运行这个示例，你需要安装以下库：

• 使用torchvision加载数据集

• 使用mldatasets、pandas、numpy和sklearn（scikit-learn）来操作数据集

• 使用torch、pytorch-lightning、efficientnet-pytorch和torchinfo模型进行预测并显示模型信息

• 使用matplotlib、seaborn、cv2、tqdm和captum来制作和可视化解释

你首先应该加载所有这些库：

import math
import os, gc
import random
import mldatasets
import numpy as np
import pandas as pd
from sklearn.preprocessing import OneHotEncoder
import torchvision
import torch
import pytorch_lightning as pl
import efficientnet_pytorch
from torchinfo import summary
import matplotlib.pyplot as plt
from matplotlib.cm import ScalarMappable 
from matplotlib.colors import LinearSegmentedColormap
import seaborn as sns
import cv2
from tqdm.notebook import tqdm
from captum import attr 

接下来，我们将加载和准备数据。

理解和准备数据

训练模型所使用的数据在 Kaggle 上公开可用（www.kaggle.com/datasets/mostafaabla/garbage-classification）。它被称为“垃圾分类”，是几个不同在线资源的汇编，包括网络爬取。它已经被分割成训练集和测试集，还附带了一个额外的较小的测试数据集，这是你的同事用来测试模型的。这些测试图像的分辨率也略高。

我们像这样从 ZIP 文件中下载数据：

dataset_file = "garbage_dataset_sample"
dataset_url = f"https://github.com/PacktPublishing/Interpretable-Machine-Learning-with-Python-2E/raw/main/datasets/{dataset_file}.zip"
torchvision.datasets.utils.download_url(dataset_url, ".")
torchvision.datasets.utils.extract_archive(f"{dataset_file}.zip",\
                                           remove_finished=True) 

它还会将 ZIP 文件提取到四个文件夹中，分别对应三个数据集和更高分辨率的测试数据集。请注意，garbage_dataset_sample只包含训练和验证数据集的一小部分。如果你想下载完整的数据集，请使用dataset_file = "garbage_dataset"。无论哪种方式，都不会影响测试数据集的大小。接下来，我们可以这样初始化数据集的转换和加载：

X_train, norm_mean = (0.485, 0.456, 0.406)
norm_std  = (0.229, 0.224, 0.225)
transform = torchvision.transforms.Compose(
    [
        torchvision.transforms.ToTensor(),
        torchvision.transforms.Normalize(norm_mean, norm_std),
    ]
)
train_data = torchvision.datasets.ImageFolder(
    f"{dataset_file}/train", transform
)
val_data = torchvision.datasets.ImageFolder(
    f"{dataset_file}/validation", transform
)
test_data = torchvision.datasets.ImageFolder(
    f"{dataset_file}/test", transform
)
test_400_data = torchvision.datasets.ImageFolder(
    f"{dataset_file}/test_400", transform
) 

上述代码所做的就是组合一系列标准转换，如归一化和将图像转换为张量。然后，它实例化与每个文件夹对应的 PyTorch 数据集——即一个用于训练、验证和测试数据集，以及更高分辨率的测试数据集（test_400_data）。这些数据集也包括转换。这样，每次从数据集中加载图像时，它都会自动进行转换。我们可以使用以下代码来验证数据集的形状是否符合我们的预期：

print(f"# Training Samples:    \t{len(train_data)}")
print(f"# Validation Samples:  \t{len(val_data)}")
print(f"# Test Samples:        \t{len(test_data)}")
print(f"Sample Dimension:      \t{test_data[0][0].shape}")
print("="*50)
print(f"# Test 400 Samples:    \t{len(test_400_data)}")
print(f"# 400 Sample Dimension:\t{test_400_data[0][0].shape}") 

上述代码输出了每个数据集中的图像数量和图像的维度。你可以看出，有超过 3,700 张训练图像，900 张验证图像和 120 张测试图像，它们的维度为 3 x 224 x 224。第一个数字对应于通道（红色、绿色和蓝色），接下来的两个数字对应于像素的宽度和高度，这是模型用于推理的。Test 400 数据集与 Test 数据集相同，只是图像的高度和宽度更大。我们不需要 Test 400 数据集进行推理，所以它不符合模型的维度要求也是可以的：

# Training Samples:        3724
# Validation Samples:      931
# Test Samples:            120
Sample Dimension:          torch.Size([3, 224, 224])
==================================================
# Test 400 Samples:        120 
# 400 Sample Dimension:    torch.Size([3, 400, 400]) 

数据准备

如果你打印(test_data[0])，你会注意到它首先会输出一个包含图像的张量，然后是一个单独的整数，我们称之为标量。这个整数是一个介于 0 到 11 之间的数字，对应于使用的标签。为了快速参考，以下是 12 个标签：

labels_l = ['battery', 'biological', 'brown-glass', 'cardboard',\
            'clothes', 'green-glass', 'metal', 'paper', 'plastic',\
            'shoes', 'trash', 'white-glass'] 

解释通常涉及从数据集中提取单个样本，以便稍后使用模型进行推理。为此，熟悉从数据集中提取任何图像，比如测试数据集的第一个样本是很重要的：

tensor, label = test_400_data[0]
img = mldatasets.tensor_to_img(tensor, norm_std, norm_mean)
plt.figure(figsize=(5,5))
plt.title(labels_l[label], fontsize=16)
plt.imshow(img)
plt.show() 
0) from the higher resolution version of the test dataset (test_400_data) and extracting the tensor and label portion from it. Then, we are using the convenience function tensor_to_img to convert the PyTorch tensor to a numpy array but also reversing the standardization that had been previously performed on the tensor. Then, we plot the image with matplotlib's imshow and use the labels_l list to convert the label into a string, which we print in the title. The result can be seen in *Figure 7.2*:

图 7.2：一个可回收碱性电池的测试样本

我们还需要执行的一个预处理步骤是对y标签进行独热编码（OHE），因为我们需要 OHE 形式来评估模型的预测性能。一旦我们初始化了OneHotEncoder，我们需要将其fit到测试标签（y_test）的数组格式中。但首先，我们需要将测试标签放入一个列表（y_test）。我们也可以用同样的方法处理验证标签，因为这些标签也便于评估：

y_test = np.array([l for _, l in test_data])
y_val = np.array([l for _, l in val_data])
ohe = OneHotEncoder(sparse=False).\
              fit(np.array(y_test).reshape(-1, 1)) 

此外，为了确保可重复性，始终这样初始化你的随机种子：

rand = 42
os.environ['PYTHONHASHSEED']=str(rand)
np.random.seed(rand)
random.seed(rand)
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
if device == 'cuda':
    torch.cuda.manual_seed(rand)
else:
    torch.manual_seed(rand) 

深度学习中确定性的实现非常困难，并且通常依赖于会话、平台和架构。如果你使用NVIDIA GPU，你可以尝试使用 PyTorch 通过命令torch.use_deterministic_algorithms(True)来避免非确定性算法。这并不保证，但如果尝试的操作无法以确定性完成，它将引发错误。如果成功，它将运行得慢得多。只有在你需要使模型结果一致时才值得这样做——例如，用于科学研究或合规性。有关可重复性和 PyTorch 的更多详细信息，请查看此处：pytorch.org/docs/stable/notes/randomness.html。

检查数据

现在，让我们看看我们的数据集中有哪些图像。我们知道训练和验证数据集非常相似，所以我们从验证数据集开始。我们可以迭代labels_l中的每个类别，并使用np.random.choice从验证数据集中随机选择一个。我们将每个图像放置在一个 4×3 的网格中，类别标签位于其上方：

plt.subplots(figsize=(14,10))
for c, category in enumerate(labels_l):
    plt.subplot(3, 4, c+1)
    plt.title(labels_l[c], fontsize=12)
    idx = np.random.choice(np.where(y_test==c)[0], 1)[0]
    im = mldatasets.tensor_to_img(test_data[idx][0], norm_std,\
                                  norm_mean)
    plt.imshow(im, interpolation='spline16')
    plt.axis("off")
plt.show() 

上一段代码生成了图 7.3。你可以看出，物品的边缘存在明显的像素化；有些物品比其他物品暗得多，而且有些图片是从奇怪的角度拍摄的：

图 7.3：验证数据集的随机样本

现在我们对测试数据集做同样的处理，以便与验证/训练数据集进行比较。我们可以使用之前的相同代码，只需将y_val替换为y_test，将val_data替换为test_data。生成的代码生成了图 7.4。你可以看出，测试集的像素化较少，物品的照明更一致，主要是从正面和侧面角度拍摄的：

图 7.4：测试数据集的随机样本

在本章中，我们不需要训练 CNN。幸运的是，客户已经为我们提供了它。

CNN 模型

其他咨询公司训练的模型是微调后的 EfficientNet 模型。换句话说，AI 咨询公司使用 EfficientNet 架构的先前训练模型，并使用垃圾分类数据集进一步训练它。这种技术被称为迁移学习，因为它允许模型利用从大型数据集（在这种情况下，来自 ImageNet 数据库的百万张图片）中学习到的先前知识，并将其应用于具有较小数据集的新任务。其优势是显著减少了训练时间和计算资源，同时保持高性能，因为它已经学会了从图像中提取有用的特征，这可以成为新任务的宝贵起点，并且只需要适应手头的特定任务。

选择 EfficientNet 是有道理的。毕竟，EfficientNet 是由 Google AI 研究人员在 2019 年引入的一组 CNN。EfficientNet 的关键创新是其复合缩放方法，这使得模型能够比其他 CNN 实现更高的准确性和效率。此外，它基于这样的观察：模型的各个维度，如宽度、深度和分辨率，以平衡的方式对整体性能做出贡献。EfficientNet 架构建立在称为 EfficientNet-B0 的基线模型之上。采用复合缩放方法创建基线模型更大、更强大的版本，同时提高网络的宽度、深度和分辨率。这产生了一系列模型，从 EfficientNet-B1 到 EfficientNet-B7，容量和性能逐渐提高。最大的模型 EfficientNet-B7 在多个基准测试中实现了最先进的性能，例如 ImageNet。

加载 CNN 模型

在我们能够加载模型之前，我们必须定义 EfficientLite 类的类——一个继承自 PyTorch Lightning 的 pl.LightningModule 的类。这个类旨在创建基于 EfficientNet 架构的定制模型，对其进行训练并执行推理。我们只需要它来执行后者，这就是为什么我们还将其修改为包含一个 predict() 函数——类似于 scikit-learn 模型，以便能够使用类似的评估函数：

class **EfficientLite**(pl.LightningModule):
    def __init__(self, lr: float, num_class: int,\
                 pretrained="efficientnet-b0", *args, **kwargs):
        super().__init__()
        self.save_hyperparameters()
        self.model = efficientnet_pytorch.EfficientNet.\
                                      from_pretrained(pretrained)
        in_features = self.model._fc.in_features
        self.model._fc = torch.nn.Linear(in_features, num_class)
    def forward(self, x):
        return self.model(x)
    def predict(self, dataset):
        self.model.eval()
        device = torch.device("cuda" if torch.cuda.is_available()\
                              else "cpu")
        with torch.no_grad():
            if isinstance(dataset, np.ndarray):
                if len(dataset.shape) == 3:
                    dataset = np.expand_dims(dataset, axis=0)
                dataset = [(x,0) for x in dataset]
            loader = torch.utils.data.DataLoader(dataset,\
                                                 batch_size=32)
            probs = None
        for X_batch, _ in tqdm(loader):
            X_batch = X_batch.to(device, dtype=torch.float32)
            logits_batch =  self.model(X_batch)
            probs_batch = torch.nn.functional.softmax(logits_batch,\
                                       dim=1).cpu().detach().numpy()
            if probs is not None:
                probs = np.concatenate((probs, probs_batch))
            else:
                probs = probs_batch
            clear_gpu_cache()
        return probs 

你会注意到这个类有三个函数：

• __init__: 这是 EfficientLite 类的构造函数。它通过使用 efficientnet_pytorch.EfficientNet.from_pretrained() 方法加载预训练的 EfficientNet 模型来初始化模型。然后，它将最后一个全连接层 (_fc) 替换为一个新创建的 torch.nn.Linear 层，该层具有相同数量的输入特征，但输出特征的数量不同，等于类别的数量 (num_class)。

• forward: 此方法定义了模型的正向传播。它接收一个输入张量 x 并将其通过模型传递，返回输出。

• predict: 此方法接收一个数据集并使用训练好的模型进行推理。它首先将模型设置为评估模式 (self.model.eval())。输入数据集被转换为具有 32 个批次的 DataLoader 对象。该方法遍历 DataLoader，处理每个数据批次，并使用 softmax 函数计算概率。在每个迭代之后调用 clear_gpu_cache() 函数以释放未使用的 GPU 内存。最后，该方法返回计算出的概率作为 numpy 数组。

如果你正在使用支持 CUDA 的 GPU，有一个名为clear_gpu_cache()的实用函数，每次进行 GPU 密集型操作时都会运行。根据你的 GPU 性能如何，你可能需要更频繁地运行它。你可以自由地使用另一个便利函数print_gpu_mem_used()来检查在任何给定时刻 GPU 内存的使用情况，或者使用print(torch.cuda.memory_summary())来打印整个摘要。接下来的代码下载预训练的 EfficientNet 模型，将模型权重加载到 EfficientLite 中，并准备模型进行推理。最后，它打印了一个摘要：

model_weights_file = **"garbage-finetuned-efficientnet-b4"**
model_url = f"https://github.com/PacktPublishing/Interpretable-Machine-Learning-with-Python-2E/raw/main/models/{model_weights_file}.ckpt"
torchvision.datasets.utils.download_url(model_url, ".")
garbage_mdl = EfficientLite.load_from_checkpoint(
    f"{model_weights_file}.ckpt"
)
garbage_mdl = garbage_mdl.to(device).eval()
print(summary(garbage_mdl)) 

代码相当直接，但重要的是要注意，我们在这个章节选择了 b4 模型，它在大小、速度和准确性方面介于 b0 和 b7 之间。你可以根据你的硬件能力更改最后一位数字，但这可能会改变本章代码的一些结果。前面的代码片段输出了以下摘要：

=======================================================================
Layer (type:depth-idx) Param # =======================================================================
EfficientLite -- 
├─EfficientNet: 1-1 – 
│ └─Conv2dStaticSamePadding: 2-1 1,296 
│ │ └─ZeroPad2d: 3-1 – 
│ └─BatchNorm2d: 2-2 96 
│ └─ModuleList: 2-3 – 
│ │ └─MBConvBlock: 3-2 2,940 
│ │ └─MBConvBlock: 3-3 1,206 
│ │ └─MBConvBlock: 3-4 11,878 
│ │ └─MBConvBlock: 3-5 18,120 
│ │ └─MBConvBlock: 3-6 18,120 
│ │ └─MBConvBlock: 3-7 18,120 
│ │ └─MBConvBlock: 3-8 25,848 
│ │ └─MBConvBlock: 3-9 57,246 
│ │ └─MBConvBlock: 3-10 57,246 
│ │ └─MBConvBlock: 3-11 57,246 
│ │ └─MBConvBlock: 3-12 70,798 
│ │ └─MBConvBlock: 3-13 197,820 
│ │ └─MBConvBlock: 3-14 197,820 
│ │ └─MBConvBlock: 3-15 197,820 
│ │ └─MBConvBlock: 3-16 197,820 
│ │ └─MBConvBlock: 3-17 197,820 
│ │ └─MBConvBlock: 3-18 240,924 
│ │ └─MBConvBlock: 3-19 413,160 
│ │ └─MBConvBlock: 3-20 413,160 
│ │ └─MBConvBlock: 3-21 413,160 
│ │ └─MBConvBlock: 3-22 413,160 
│ │ └─MBConvBlock: 3-23 413,160 
│ │ └─MBConvBlock: 3-24 520,904 
│ │ └─MBConvBlock: 3-25 1,159,332 
│ │ └─MBConvBlock: 3-26 1,159,332 
│ │ └─MBConvBlock: 3-27 1,159,332 
│ │ └─MBConvBlock: 3-28 1,159,332 
│ │ └─MBConvBlock: 3-29 1,159,332 
│ │ └─MBConvBlock: 3-30 1,159,332 
│ │ └─MBConvBlock: 3-31 1,159,332 
│ │ └─MBConvBlock: 3-32 1,420,804 
│ │ └─MBConvBlock: 3-33 3,049,200 
│ └─Conv2dStaticSamePadding: 2-4 802,816 
│ │ └─Identity: 3-34 – 
│ └─BatchNorm2d: 2-5 3,584 
│ └─AdaptiveAvgPool2d: 2-6 – 
│ └─Dropout: 2-7 – 
│ └─Linear: 2-8 21,516 
│ └─MemoryEfficientSwish: 2-9 
======================================================================
Total params: 17,570,132 
Trainable params: 17,570,132 
Non-trainable params: 0 ====================================================================== 

它几乎包含了我们需要的关于模型的所有信息。它有两个自定义卷积层（Conv2dStaticSamePadding），每个卷积层后面跟着一个批归一化层（BatchNorm2d）和 32 个MBConvBlock模块。

网络还有一个 Swish 激活函数的内存高效实现（MemoryEfficientSwish），就像所有激活函数一样，它将非线性引入模型。它是平滑且非单调的，有助于它更快地收敛，同时学习更复杂和细微的模式。它还有一个全局平均池化操作（AdaptiveAvgPool2d），它减少了特征图的空间维度。然后有一个用于正则化的第一个Dropout层，后面跟着一个将节点数从 1792 减少到 12 的完全连接层（Linear）。Dropout 通过在每个更新周期中使一部分神经元不活跃来防止过拟合。如果你想知道每个层之间的输出形状是如何减少的，可以将input_size输入到摘要中——例如summary(garbage_mdl, input_size=(64, 3, 224, 224))——因为网络是针对 64 个批次的尺寸设计的。如果你对这些术语不熟悉，不要担心。我们稍后会重新讨论它们。

使用传统解释方法评估 CNN 分类器

我们将首先使用evaluate_multiclass_mdl函数和验证数据集来评估模型。参数包括模型（garbage_mdl）、我们的验证数据（val_data）、类别名称（labels_l）以及编码器（ohe）。最后，我们不会绘制 ROC 曲线（plot_roc=False）。此函数返回预测标签和概率，我们可以将它们存储在变量中以供以后使用：

y_val_pred, y_val_prob = mldatasets.evaluate_multiclass_mdl(
    garbage_mdl, val_data,\
    class_l=labels_l, ohe=ohe, plot_roc=False
) 

前面的代码生成了带有混淆矩阵的图 7.5和每个类别的性能指标的图 7.6：

图 7.5：验证数据集的混淆矩阵

尽管图 7.5中的混淆矩阵似乎表明分类完美，但一旦你看到图 7.6中的精确率和召回率分解，你就可以知道模型在金属、塑料和白色玻璃方面存在问题：

图 7.6：验证数据集的分类报告

如果你使用最优的超参数对模型进行足够的轮次训练，你可以期望模型总是达到100%的训练准确率。接近完美的验证准确率更难实现，这取决于这两个值之间的差异。我们知道验证数据集只是来自同一集合的图像样本，所以达到 94.7%并不特别令人惊讶。

plot_roc=True) but only the averages, and not on a class-by-class basis (plot_roc_class=False) because there are only four pictures per class. Given the small number of samples, we can display the numbers in the confusion matrix rather than percentages (pct_matrix=False):

y_test_pred, y_test_prob = mldatasets.evaluate_multiclass_mdl(
    garbage_mdl, test_data,\
    class_l=labels_l, ohe=ohe,\
    plot_roc=True, plot_roc_class=False, pct_matrix=False
) 
generated the ROC curve in *Figure 7.7*, the confusion matrix in *Figure 7.8*, and the classification report in *Figure 7.9*:

图 7.7：测试数据集的 ROC 曲线

测试 ROC 图(图 7.7)显示了宏平均和微平均的 ROC 曲线。这两者的区别在于它们的计算方式。宏度量是独立地对每个类别进行计算然后平均，对待每个类别不同，而微平均则考虑了每个类别的贡献或代表性；一般来说，微平均更可靠。

图 7.8：测试数据集的混淆矩阵

如果我们查看图 7.8中的混淆矩阵，我们可以看出只有生物、绿色玻璃和鞋子得到了 10/10 的分类。然而，很多物品被错误地分类为生物和鞋子。另一方面，很多物品经常被错误分类，比如金属、纸张和塑料。许多物品在形状或颜色上相似，所以你可以理解为什么会这样，但金属怎么会和白色玻璃混淆，或者纸张会和电池混淆呢？

图 7.9：测试数据集的预测性能指标

在商业环境中讨论分类模型时，利益相关者通常只对一个数字感兴趣：准确率。很容易让这个数字驱动讨论，但其中有很多细微差别。例如，令人失望的测试准确率（68.3%）可能意味着很多事情。这可能意味着六个类别得到了完美的分类，而其他所有类别都没有，或者 12 个类别只有一半被错误分类。可能发生的事情有很多。

在任何情况下，处理多类分类问题时，即使准确率低于 50%也可能不像看起来那么糟糕。考虑到无信息率代表了在数据集中总是预测最频繁类别的朴素模型所能达到的准确率。它作为一个基准，确保开发出的模型提供了超越这种简单方法的见解。并且，如果数据集被平均分成 12 类，那么无信息率可能大约是 8.33%（100%/12 类），所以 68%仍然比这高得多。实际上，距离 100%的差距还要小！对于一个机器学习从业者来说，这意味着如果我们仅仅根据测试准确率结果来判断，模型仍在学习一些有价值的东西，这些是可以进一步改进的。

在任何情况下，测试数据集在图 7.9中的预测性能指标与我们在混淆矩阵中看到的一致。生物类别召回率高但精确率低，而金属、纸张、塑料和垃圾的召回率都很低。

确定要关注的错误分类

我们已经注意到一些有趣的错误分类，我们可以集中关注：

• 金属的假阳性：测试数据集中有 120 个样本中的 16 个被错误分类为金属。这是所有错误分类的 42%！模型为什么如此容易将金属与其他垃圾混淆，这是怎么回事？

• 塑料的假阴性：70%的所有真实塑料样本都被错误分类。因此，塑料在所有材料中除了垃圾之外，召回率最低。很容易理解为什么垃圾分类如此困难，因为它极其多样，但不是塑料。

我们还应该检查一些真实阳性，以对比这些错误分类。特别是电池，因为它们作为金属和塑料有很多假阳性，以及白色玻璃，因为它 30%的时间作为金属有假阴性。由于金属的假阳性很多，我们应该将它们缩小到仅仅是电池的。

为了可视化前面的任务，我们可以创建一个 DataFrame（preds_df），其中包含一个列的真实标签（y_true）和另一个列的预测标签。为了了解模型对这些预测的确定性，我们可以创建另一个包含概率的 DataFrame（probs_df）。我们可以为这些概率生成列总计，以便根据模型在所有样本中最确定哪个类别来排序列。然后，我们可以将我们的预测 DataFrame 与概率 DataFrame 的前 12 列连接起来：

preds_df = pd.DataFrame({'y_true':[labels_l[o] for o in y_test],\
                         'y_pred':y_test_pred})
probs_df = pd.DataFrame(y_test_prob*100).round(1)
probs_df.loc['Total']= probs_df.sum().round(1)
probs_df.columns = labels_l
probs_df = probs_df.sort_values('Total', axis=1, ascending=False)
probs_df.drop(['Total'], axis=0, inplace=True)
probs_final_df = probs_df.iloc[:,0:12]
preds_probs_df = pd.concat([preds_df, probs_final_df], axis=1) 

现在我们输出感兴趣的预测实例的 DataFrame，并对其进行颜色编码。一方面，我们有金属的假阳性，另一方面，我们有塑料的假阴性。但我们还有电池和白色玻璃的真实阳性。最后，我们将所有超过 50%的概率加粗，并将所有 0%的概率隐藏起来，这样更容易发现任何高概率的预测：

num_cols_l = list(preds_probs_df.columns[2:])
num_fmt_dict = dict(zip(num_cols_l, ["{:,.1f}%"]*len(num_cols_l)))
preds_probs_df[
    (preds_probs_df.y_true!=preds_probs_df.y_pred)
    | (preds_probs_df.y_true.isin(['battery', 'white-glass']))
].style.format(num_fmt_dict).apply(
    lambda x: ['background: lightgreen' if (x[0] == x[1])\
                else '' for i in x], axis=1
).apply(
    lambda x: ['background: orange' if (x[0] != x[1] and\
                x[1] == 'metal' and x[0] == 'battery')\
                else '' for i in x], axis=1
).apply(
    lambda x: ['background: yellow' if (x[0] != x[1] and\
                                        x[0] == 'plastic')\
                else '' for i in x], axis=1
).apply(
    lambda x: ['font-weight: bold' if isinstance(i, float)\
                                                 and i >= 50\
                else '' for i in x], axis=1
).apply(
    lambda x: ['color:transparent' if i == 0.0\
                else '' for i in x], axis=1) 
Figure 7.10. We can tell by the highlights which are the metal false positives and the plastic false negatives, as well as which would be the true positives: #0-6 for battery, and #110-113 and #117-119 for white glass:

图 7.10：测试数据集中所有 38 个错误分类、选定的真实正例及其真实和预测标签，以及它们的预测概率的表格

我们可以使用以下代码轻松地将这些实例的索引存储在列表中。这样，为了未来的参考，我们可以遍历这些列表来评估单个预测，或者用它们来对整个组执行解释任务。正如你所看到的，我们有所有四个组的列表：

plastic_FN_idxs = preds_df[
    (preds_df['y_true'] !=preds_df['y_pred'])
    & (preds_df['y_true'] == 'plastic')
].index.to_list()
metal_FP_idxs = preds_df[
    (preds_df['y_true'] != preds_df['y_pred'])
    & (preds_df['y_pred'] == 'metal')
    & (preds_df['y_true'] == 'battery')
].index.to_list()
battery_TP_idxs = preds_df[
    (preds_df['y_true'] ==preds_df['y_pred'])
    & (preds_df['y_true'] == 'battery')
].index.to_list()
wglass_TP_idxs = preds_df[
    (preds_df['y_true'] == preds_df['y_pred'])
    & (preds_df['y_true'] == 'white-glass')
].index.to_list() 

现在我们已经预处理了所有数据，模型已完全加载并列出要调试的预测组。现在我们可以继续前进。让我们开始解释！

使用基于激活的方法可视化学习过程

在我们开始讨论激活、层、过滤器、神经元、梯度、卷积、核以及构成卷积神经网络（CNN）的所有神奇元素之前，让我们首先简要回顾一下 CNN 的机制，特别是其中一个机制。

卷积层是 CNN 的基本构建块，它是一个顺序神经网络。它通过可学习的过滤器对输入进行卷积，这些过滤器相对较小，但会在特定的距离或步长上应用于整个宽度、高度和深度。每个过滤器产生一个二维的激活图（也称为特征图）。之所以称为激活图，是因为它表示图像中激活的位置——换句话说，特定“特征”所在的位置。在这个上下文中，特征是一个抽象的空间表示，在处理过程的下游，它反映在完全连接（线性）层的所学权重中。例如，在垃圾 CNN 案例中，第一个卷积层有 48 个过滤器，3 × 3 的核，2 × 2 的步长和静态填充，这确保输出图保持与输入相同的大小。过滤器是模板匹配的，因为当在输入图像中找到某些模式时，它们最终会在激活图中激活区域。

但在我们到达完全连接层之前，我们必须减小过滤器的尺寸，直到它们达到可工作的尺寸。例如，如果我们展平第一个卷积的输出（48 × 112 × 112），我们就会有超过 602,000 个特征。我想我们都可以同意，这会太多以至于无法输入到完全连接层中。即使我们使用了足够的神经元来处理这项工作负载，我们可能也没有捕捉到足够的空间表示，以便神经网络能够理解图像。因此，卷积层通常与池化层配对，池化层对输入进行下采样——换句话说，它们减少了数据的维度。在这种情况下，有一个自适应平均池化层（AdaptiveAvgPool2d），它在所有通道上执行平均，以及许多在Mobile Inverted Bottleneck Convolution Blocks（MBConvBlock）内的池化层。

顺便提一下，MBConvBlock、Conv2dStaticSamePadding和BatchNorm2d是 EfficientNet 架构的构建块。这些组件共同工作，创建了一个高度高效且准确的卷积神经网络：

• MBConvBlock: 形成 EfficientNet 架构核心的移动倒置瓶颈卷积块。在传统的卷积层中，过滤器同时应用于所有输入通道，导致计算量很大，但MBConvBlocks将这个过程分为两个步骤：首先，它们应用深度卷积，分别处理每个输入通道，然后使用点卷积（1 x 1）来结合来自不同通道的信息。因此，在 B0 的MBConvBlock模块中，有三个卷积层：一个深度卷积，一个点卷积（称为项目卷积），以及在某些块中的另一个点卷积（称为扩展卷积）。然而，第一个块只包含两个卷积层（深度卷积和项目卷积），因为它没有扩展卷积。对于 B4，架构类似，但每个块中堆叠的卷积更多，MBConvBlocks的数量也翻倍。自然地，B7 有更多的块和卷积层。对于 B4，总共有 158 次卷积操作分布在 32 个MBConvBlocks之间。

• Conv2dStaticSamePadding: 与传统的卷积层（如Conv2d）不同，这些层不会减少维度。它确保输入和输出特征图具有相同的空间维度。

• BatchNorm2d: 批标准化层，通过归一化输入特征来帮助稳定和加速训练，这有助于在训练过程中保持输入特征的分布一致性。

一旦执行了超过 230 次的卷积和池化操作，我们得到一个更易于处理的扁平化输出：1,792 个特征，全连接层将这些特征转换为 12 个，利用softmax激活函数，为每个类别输出介于 0 和 1 之间的概率。在垃圾 CNN 中，有一个dropout层用于帮助正则化训练。我们可以完全忽略这一点，因为在推理过程中，它们是被忽略的。

如果这还不够清晰，不要担心！接下来的部分将通过激活、梯度和扰动直观地展示网络可能学习或未学习的图像表示方式。

中间激活

对于推理，图像通过网络的输入，预测通过输出穿过每个层。然而，具有顺序和分层架构的一个优点是我们可以提取任何层的输出，而不仅仅是最终层。中间激活是任何卷积或池化层的输出。它们是激活图，因为激活函数应用后，亮度较高的点映射到图像的特征。在这种情况下，模型在所有卷积层上使用了 ReLU，这就是激活点的原理。我们只对卷积层的中间激活感兴趣，因为池化层只是这些层的下采样版本。为什么不去看更高分辨率的版本呢？

随着滤波器宽度和高度的减小，学习到的表示将会更大。换句话说，第一个卷积层可能关于细节，如纹理，下一个关于边缘，最后一个关于形状。然后我们必须将卷积层的输出展平，以便将其输入到从那时起接管的多层感知器。

我们现在要做的是提取一些卷积层的激活。在 B4 中，有 158 个，所以我们不能全部做！为此，我们将使用model.children()获取第一层的层，并遍历它们。我们将从这个顶层将两个Conv2dStaticSamePadding层添加到conv_layers列表中。但我们会更深入，将ModuleList层中的前六个MBConvBlock层的第一个卷积层也添加进去。最后，我们应该有八个卷积层——中间的六个属于 Mobile Inverted Bottleneck Convolution 块：

conv_layers = []
model_children = list(garbage_mdl.model.children())
for model_child in model_children:
    if (type(model_child) ==\
                efficientnet_pytorch.utils.Conv2dStaticSamePadding):
        conv_layers.append(model_child)
    elif (type(model_child) == torch.nn.modules.container.ModuleList):
        module_children = list(model_child.children())
        module_convs = []
        for module_child in module_children:
            module_convs.append(list(module_child.children())[0])
        conv_layers.extend(module_convs[:6])
print(conv_layers) 

在我们遍历所有它们，为每个卷积层生成激活图之前，让我们为单个滤波器和层做一下：

idx = battery_TP_idxs[0]
tensor = test_data[idx][0][None, :].to(device)
label = y_test[idx]
method = attr.LayerActivation(garbage_mdl, conv_layers[layer]
attribution = method.attribute(tensor).detach().cpu().numpy()
print(attribution.shape) 
tensor for the first battery true positive (battery_TP_idxs[0]). Then, it initializes the LayerActivation attribution method with the model (garbage_mdl) and the first convolutional layer (conv_layers[0]). Using the attribute function, it creates an attribution with this method. For the shape of the attribution, we should get (1, 48, 112, 112). The tensor was for a single image, so it makes sense that the first number is a one. The next number corresponds to the number of filters, followed by the width and height dimensions of each filter. Regardless of the kind of attribution, the numbers inside each attribution relate to how a pixel in the input is seen by the model. Interpretation varies according to the method. However, generally, it is interpreted that higher numbers mean more of an impact on the outcome, but attributions may also have negative numbers, which mean the opposite.

让我们可视化第一个滤波器，但在我们这样做之前，我们必须决定使用什么颜色图。颜色图将决定将不同数字分配给哪些颜色作为渐变。例如，以下颜色图将白色分配给0（十六进制中的#ffffff），中等灰色分配给0.25，黑色（十六进制中的#000000）分配给1，这些颜色之间有一个渐变：

cbinary_cmap = LinearSegmentedColormap.from_list('custom binary',
                                                 [(0, '#ffffff'),
                                                  (0.25, '#777777'),
                                                  (1, '#000000')]) 

你也可以使用matplotlib.org/stable/tutorials/colors/colormaps.html上的任何命名颜色图，而不是使用你自己的。接下来，让我们像这样绘制第一个滤波器的属性图：

filter = 0
filter_attr = attribution[0,filter]
filter_attr = mldatasets.apply_cmap(filter_attr, cbinary_cmap, 'positive')
y_true = labels_l[label]
y_pred = y_test_pred[idx]
fig, ax = plt.subplots(1, 1, figsize=(8, 6))
fig.suptitle(f"Actual label: {y_true}, Predicted: {y_pred}", fontsize=16)
ax.set_title(
    f"({method.get_name()} Attribution for Filter #{filter+1} for\
        Convolutional Layer #{layer+1})",
    fontsize=12
)
ax.imshow(filter_attr)
ax.grid(False)
fig.colorbar(
    ScalarMappable(norm='linear', cmap=cbinary_cmap),
    ax=ax,
    orientation="vertical"
)
plt.show() 
Figure 7.11:

图 7.11：第一个真实正样本电池样本的第一个卷积层的第一个滤波器的中间激活图

如你在图 7.11中可以看到，第一个滤波器的中间激活似乎在寻找电池的边缘和最突出的文本。

接下来，我们将遍历所有计算层和每个电池，并可视化每个的归因。现在，一些这些归因操作可能计算成本很高，因此在这些操作之间清除 GPU 缓存（clear_gpu_cache()）是很重要的：

for l, layer in enumerate(conv_layers):
    layer = conv_layers[l]
    method = attr.LayerActivation(garbage_mdl, layer)
    for idx in battery_TP_idxs:
        orig_img = mldatasets.tensor_to_img(test_400_data[idx][0],\
                                            norm_std, norm_mean,\
                                            to_numpy=True)
        tensor = test_data[idx][0][None, :].to(device)
        label = int(y_test[idx])
        attribution = method.attribute(tensor).detach().cpu().numpy()
        viz_img =  mldatasets.**create_attribution_grid**(attribution,\
                            cmap='copper', cmap_norm='positive')
        y_true = labels_l[label]
        y_pred = y_test_pred[idx]
        probs_s = probs_df.loc[idx]
        name = method.get_name()
        title = f'CNN Layer #{l+1} {name} Attributions for Sample #{idx}'
        mldatasets.**compare_img_pred_viz**(orig_img, viz_img, y_true,\
                                        y_pred, probs_s, title=title)
    clear_gpu_cache() 
look fairly familiar. Where it’s different is that it’s placing every attribution map for every filter in a grid (viz_img) with create_attribition_grid. It could just then display it with plt.imshow as before, but instead, we will leverage a utility function called compare_img_pred_viz to visualize the attribution(s) side by side with the original image (orig_img). It also takes the sample’s actual label (y_true) and predicted label (y_pred). Optionally, we can provide a pandas series with the probabilities for this prediction (probs_s) and a title. It generates 56 images in total, including *Figures 7.12*, *7.13*, and *7.14*.

如您从图 7.12中可以看出，第一层卷积似乎在捕捉电池的字母以及其轮廓：

图 7.12：电池#4 的第一卷积层的中间激活

然而，图 7.13显示了网络如何通过第四层卷积更好地理解电池的轮廓：

图 7.13：电池#4 的第四卷积层的中间激活

在图 7.14中，最后一层卷积层难以解释，因为这里有 1,792 个 7 像素宽和高的过滤器，但请放心，那些微小的图中编码了一些非常高级的特征：

图 7.14：电池#4 的最后一层卷积层的中间激活

提取中间激活可以为你提供基于样本的某些洞察。换句话说，它是一种局部模型解释方法。这绝对不是唯一的逐层归因方法。Captum 有超过十种层归因方法：github.com/pytorch/captum#about-captum。

使用基于梯度的归因方法评估误分类

基于梯度的方法通过 CNN 的前向和反向传递计算每个分类的归因图。正如其名所示，这些方法利用反向传递中的梯度来计算归因图。所有这些方法都是局部解释方法，因为它们只为每个样本推导出一个解释。顺便提一下，在这个上下文中，归因意味着我们将预测标签归因于图像的某些区域。在学术文献中，它们也常被称为敏感性图。

要开始，我们首先需要创建一个数组，包含测试数据集（test_data）中所有我们的误分类样本（X_misclass），使用所有我们感兴趣的误分类的合并索引（misclass_idxs）。由于误分类并不多，我们正在加载它们的一个批次（next）：

misclass_idxs = metal_FP_idxs + plastic_FN_idxs[-4:]
misclass_data = torch.utils.data.Subset(test_data, misclass_idxs)
misclass_loader = torch.utils.data.DataLoader(misclass_data,\
                                              batch_size = 32)
X_misclass, y_misclass = next(iter(misclass_loader))
X_misclass, y_misclass = X_misclass.to(device), y_misclass.to(device) 

下一步是创建一个我们可以重用的实用函数来获取任何方法的归因图。可选地，我们可以使用名为NoiseTunnel的方法（github.com/pytorch/captum#getting-started）来平滑地图。我们将在稍后更详细地介绍这种方法：

def get_attribution_maps(**method**, model, device,X,y=None,\
                         init_args={}, nt_type=None, nt_samples=10,\
                         stdevs=0.2, **kwargs):
    attr_maps_size = tuple([0] + list(X.shape[1:]))
    attr_maps = torch.empty(attr_maps_size).to(device)
    **attr_method** = **method**(model, **init_args)
    if nt_type is not None:
        noise_tunnel = attr.NoiseTunnel(attr_method)
        nt_attr_maps = torch.empty(attr_maps_size).to(device)
    for i in tqdm(range(len(X))):
        X_i = X[i].unsqueeze(0).requires_grad_()
        model.zero_grad()
        extra_args = {**kwargs}
        if y is not None:
            y_i = y[i].squeeze_()
            extra_args.update({"target":y_i})

        attr_map = **attr_method.attribute**(X_i, **extra_args)
        attr_maps = torch.cat([attr_maps, attr_map])
        if nt_type is not None:
            model.zero_grad()
            nt_attr_map = noise_tunnel.attribute(
                X_i, nt_type=nt_type, nt_samples=nt_samples,\
                stdevs=stdevs, nt_samples_batch_size=1, **extra_args)
            nt_attr_maps = torch.cat([nt_attr_maps, nt_attr_map])
        clear_gpu_cache()
    if nt_type is not None:
        return attr_maps, nt_attr_maps
    return attr_maps 

上述代码可以为给定模型和设备的任何 Captum 方法创建归因图。为此，它需要图像的张量X及其相应的标签y。标签是可选的，只有在归因方法是针对特定目标时才需要 - 大多数方法都是。大多数归因方法（attr_method）仅使用模型初始化，但一些需要一些额外的参数（init_args）。它们通常在用attribute函数生成归因时具有最多的参数，这就是为什么我们在get_attribution_maps函数中收集额外的参数（**kwargs），并将它们放在这个调用中。

需要注意的一个重要事项是，在这个函数中，我们遍历X张量中的所有样本，并为每个样本独立创建属性图。这通常是不必要的，因为属性方法都配备了同时处理一批数据的能力。然而，存在硬件无法处理整个批次的风险，在撰写本文时，非常少的方法带有internal_batch_size参数，这可能会限制一次可以处理的样本数量。我们在这里所做的是本质上等同于每次都将这个数字设置为1，以努力确保我们不会遇到内存问题。然而，如果你有强大的硬件，你可以重写函数以直接处理X和y张量。

接下来，我们将执行我们的第一个基于梯度的归因方法。

显著性图

显著性图依赖于梯度的绝对值。直觉上，它会找到图像中可以扰动最少且输出变化最大的像素。它不执行扰动，因此不验证假设，而绝对值的使用阻止它找到相反的证据。

这种首次提出的显著性图方法在当时具有开创性，并激发了许多不同的方法。它通常被昵称为“vanilla”，以区别于其他显著性图。

使用我们的get_attribution_maps函数为所有误分类的样本生成显著性图相对简单。你所需要的是 Captum 归因方法（attr.Saliency）、模型（garbage_mdl）、设备以及误分类样本的张量（X_misclass和y_misclass）：

saliency_maps = get_attribution_maps(attr.Saliency, garbage_mdl,\
                                     device, X_misclass, y_misclass) 

我们可以绘制其中一个显著性图的输出，第五个，与样本图像并排显示以提供上下文。Matplotlib 可以通过subplots网格轻松完成此操作。我们将创建一个 1 × 3 的网格，并将样本图像放在第一个位置，其显著性热图放在第二个位置，第三个位置是叠加在一起的。就像我们之前对归因图所做的那样，我们可以使用tensor_to_img将图像转换为numpy数组，同时应用归因的调色板。它默认使用 jet 调色板（cmap='jet'）使显著的区域看起来更加突出：

pos = 4
orig_img = mldatasets.tensor_to_img(X_misclass[pos], norm_std,\
                                    norm_mean, to_numpy=True)
attr_map = mldatasets.tensor_to_img(
    saliency_maps[pos], to_numpy=True,\
    cmap_norm='positive'
)
fig, axs = plt.subplots(1, 3, figsize=(15,5))
axs[0].imshow(orig_img)
axs[0].grid(None)
axs[0].set_title("Original Image")
axs[1].imshow(attr_map)
axs[1].grid(None)
axs[1].set_title("Saliency Heatmap")
axs[2].imshow(np.mean(orig_img, axis=2), cmap="gray")
axs[2].imshow(attr_map, alpha=0.6)
axs[2].grid(None)
axs[2].set_title("Saliency Overlayed")
idx = misclass_idxs[pos]
y_true = labels_l[int(y_test[idx])]
y_pred = y_test_pred[idx]
plt.suptitle(f"Actual label: {y_true}, Predicted: {y_pred}")
plt.show() 

上述代码生成了图 7.15中的图表：

图 7.15：将塑料误分类为生物废物的显著性图

图 7.15中的样本图像看起来像是被撕碎的塑料，但预测结果是生物废物。标准的显著性图将这个预测主要归因于塑料上较平滑、较暗的区域。看起来是缺乏镜面高光让模型产生了偏差，但通常，较旧的破损塑料会失去光泽。

镜面高光是在物体表面反射光线时出现的明亮光点。它们通常是光源的直接反射，并且在光滑或光亮的表面上更为明显，例如金属、玻璃或水。

引导 Grad-CAM

要讨论引导 Grad-CAM，我们首先应该讨论CAM，它代表类别激活图。CAM 的工作方式是移除除了最后一层全连接层之外的所有层，并用全局平均池化（GAP）层替换最后一个最大池化层。GAP 层计算每个特征图的平均值，将其减少到每个图的单个值，而最大池化层通过从图的一个局部区域中的值集中选择最大值来减小特征图的大小。例如，在这个案例中：

1. 最后一个卷积层输出一个1792 × 7 × 7的张量。

2. GAP 通过仅平均这个张量的最后两个维度来减少维度，产生一个1792 × 1 × 1的张量。

3. 然后，它将这个结果输入到一个有 12 个神经元的全连接层中，每个神经元对应一个类别。

4. 一旦重新训练了一个 CAM 模型并通过样本图像通过 CAM 模型，它将从最后一层（一个1792 × 12的张量）中提取与预测类别相对应的值（一个1792 × 1的张量）。

5. 然后，你计算最后一个卷积层输出（1792 × 7 × 7）与权重张量（1792 x 1）的点积。

6. 这个加权的总和将结束于一个1 × 7 × 7的张量。

7. 通过双线性插值将其拉伸到1 × 224 × 224，这变成了一个上采样后的激活图。当你上采样数据时，你增加了其维度。

CAM 背后的直觉是，CNN 在卷积层中本质上保留了空间细节，但遗憾的是，这些细节在全连接层中丢失了。实际上，最后一个卷积层中的每个滤波器代表不同空间位置上的视觉模式。一旦加权，它们就代表了整个图像中最显著的区域。然而，要应用 CAM，你必须彻底修改模型并重新训练它，而且有些模型并不容易适应这种修改。

如其名称所示，Grad-CAM 是一个类似的概念，但避免了修改和重新训练的麻烦，并使用梯度代替——具体来说，是关于卷积层激活图的类别分数（在 softmax 之前）的梯度。对这些梯度执行 GAP 操作以获得神经元重要性权重。然后，我们使用这些权重计算激活图的加权线性组合，随后是 ReLU。ReLU 非常重要，因为它确保定位只对结果产生正面影响的特征。像 CAM 一样，它通过双线性插值上采样以匹配图像的尺寸。

Grad-CAM 也有一些缺点，例如无法识别多个发生或由预测类别表示的物体的全部。像 CAM 一样，激活图的分辨率可能受到最终卷积层维度的限制，因此需要上采样。

因此，我们使用引导 Grad-CAM。引导 Grad-CAM 是 Grad-CAM 和引导反向传播的结合。引导反向传播是另一种可视化方法，它计算目标类别相对于输入图像的梯度，但它修改了反向传播过程，只传播正激活的正梯度。这导致了一个更高分辨率、更详细的可视化。这是通过将 Grad-CAM 热图（上采样到输入图像分辨率）与引导反向传播结果进行逐元素乘法来实现的。输出是一个可视化，强调给定类别在图像中最相关的特征，比单独的 Grad-CAM 具有更高的空间细节。

使用我们的get_attribution_maps函数为所有误分类样本生成 Grad-CAM 归因图。你需要的是 Captum 归因方法（attr.GuidedGradCam）、模型（garbage_mdl）、设备以及误分类样本的张量（X_misclass和y_misclass），并在方法初始化参数中，一个用于计算 Grad-CAM 归因的层：

gradcam_maps = get_attribution_maps(
    attr.GuidedGradCam, garbage_mdl, device, X_misclass,\
    y_misclass, init_args={'layer':conv_layers[3]}
) 

注意，我们并没有使用最后一层（可以用7或-1索引）而是第四层（3）。这样做只是为了保持事情有趣，但我们也可以更改它。接下来，让我们像之前一样绘制归因图。代码几乎相同，只是将saliency_maps替换为gradcam_maps。输出结果如图7.16所示。

图 7.16：将塑料误分类为生物废物的引导 Grad-CAM 热图

如你在图 7.16中观察到的，与显著性归因图一样，类似的平滑哑光区域被突出显示，除了引导 Grad-CAM 产生一些亮区和边缘。

对所有这些内容都要持保留态度。在 CNN 解释领域，仍然存在许多持续的争论。研究人员仍在提出新的和更好的方法，甚至对于大多数用例几乎完美的技术仍然存在缺陷。关于类似 CAM 的方法，有许多新的方法，例如Score-CAM、Ablation-CAM和Eigen-CAM，它们提供了类似的功能，但不需要依赖梯度，而梯度可能是不稳定的，因此有时是不可靠的。我们在这里不会讨论它们，因为当然，它们不是基于梯度的！但是，尝试不同的方法以查看哪些适用于您的用例是有益的。

集成梯度

集成梯度（IG），也称为路径积分梯度，是一种不限于 CNN 的技术。您可以将它应用于任何神经网络架构，因为它计算了输出相对于输入的梯度，这些梯度是在从基线到实际输入之间的路径上平均计算的。它对卷积层的存在不敏感。然而，它需要定义一个基线，这个基线应该传达信号缺失的概念，比如一个均匀着色的图像。在实践中，特别是对于 CNN 来说，这表示零基线，对于每个像素来说，通常意味着一个完全黑色的图像。尽管名称暗示了使用路径积分，但积分并不是计算的，而是用足够小的区间内的求和来近似，对于一定数量的步骤。对于 CNN 来说，这意味着它使输入图像的变体逐渐变暗或变亮，直到它成为对应于预定义步骤数的基线。然后它将这些变体输入 CNN，为每个变体计算梯度，并取平均值。IG 是图像与梯度平均值之间的点积。

与 Shapley 值一样，IG 建立在坚实的数学理论基础上。在这种情况下，它是线积分的基本定理。IG 方法的数学证明确保了所有特征的归因之和等于模型在输入数据上的预测与在基线输入上的预测之间的差异。除了他们称之为完备性的这种属性之外，还有线性保持、对称保持和敏感性。我们在这里不会描述这些属性中的每一个。然而，重要的是要注意，一些解释方法满足显著的数学属性，而其他方法则从实际应用中证明了它们的有效性。

除了 IG 之外，我们还将利用NoiseTunnel对样本图像进行小的随机扰动——换句话说，就是添加噪声。它多次创建相同样本图像的不同噪声版本，然后计算每个版本的归因方法。然后它对这些归因进行平均，这可能是使归因图更加平滑的原因，这就是为什么这种方法被称为SmoothGrad。

但等等，你可能要问：那它不应该是一种基于扰动的算法吗？！在这本书中，我们之前已经处理了几种基于扰动的算法，从 SHAP 到锚点，它们共有的特点是它们扰动输入以测量对输出的影响。SmoothGrad 并不测量对输出的影响。它只帮助生成一个更鲁棒的归因图，因为扰动输入的平均归因应该会生成更可靠的归因图。我们进行交叉验证来评估机器学习模型也是出于同样的原因：在不同分布的测试数据集上执行的平均指标会生成更好的指标。

对于 IG，我们将使用与 Saliency 相同的非常相似的代码，除了我们将添加几个与NoiseTunnel相关的参数，例如噪声隧道的类型（nt_type='smoothgrad'）、用于生成的样本变化（nt_samples=20）以及添加到每个样本中的随机噪声的量（以标准差计stdevs=0.2）。我们会发现，生成的置换样本越多，效果越好，但达到一定程度后，效果就不会有太大变化。然而，噪声过多也是一种情况，如果你使用得太少，就不会有任何效果：

ig_maps, smooth_ig_maps = get_attribution_maps(
    attr.IntegratedGradients, garbage_mdl, device, X_misclass,\
    y_misclass, nt_type='smoothgrad', nt_samples=20, stdevs=0.2
) 

我们还可以选择性地定义 IG 的步数（n_steps）。默认设置为50，我们还可以修改基线，默认情况下是一个全零的张量。正如我们使用 Grad-CAM 所做的那样，我们可以将第一个样本图像与 IG 图并排显示，但这次，我们将修改代码以在第三个位置绘制 SmoothGrad 集成梯度（smooth_ig_maps），如下所示：

nt_attr_map = mldatasets.tensor_to_img(
    smooth_ig_maps[pos], to_numpy=True, cmap_norm='positive'
)
axs[2].imshow(nt_attr_map)
axs[2].grid(None)
axs[2].set_title("SmoothGrad Integrated Gradients") 
Figure 7.17:

图 7.17：将塑料误分类为生物垃圾的集成梯度热图

在图 7.17 的 IG 热图中的区域与显著性图和引导 Grad-CAM 图检测到的许多区域相吻合。然而，在明亮的黄色区域以及棕色的阴影区域中，有更多的强归因簇，这与某些食物被丢弃时的外观（如香蕉皮和腐烂的叶状蔬菜）一致。另一方面，明亮的橙色和绿色区域则不是这样。

至于 SmoothGrad IG 热图，与不平滑的 IG 热图相比，这张图非常不同。这并不总是如此；通常，它只是更平滑的版本。可能发生的情况是0.2噪声对归因的影响过大，或者 20 个扰动样本不够。然而，很难说，因为也有可能 SmoothGrad 更准确地描绘了真实的故事。

我们现在不会做这件事，但你可以直观地“调整”stdevs和nt_samples参数。你可以尝试使用更少的噪声和更多的样本，使用一系列组合，例如0.1和80，以及0.15和40，试图找出它们之间是否存在共性。你所选择的那个最能清楚地描绘出这个一致的故事。SmoothGrad 的一个缺点是必须定义最优参数。顺便提一下，IG 在定义基线和步数（n_steps）方面也存在相同的问题。默认的基线在输入图像太大或太小时将不起作用，因此必须更改，IG 论文的作者建议 20-300 步将使积分在 5%以内。

奖励方法：DeepLIFT

IG 有一些批评者，他们已经创建了避免使用梯度的类似方法，例如DeepLIFT。IG 对零值梯度和梯度的不连续性可能很敏感，这可能导致误导性的归因。但这些指向的是所有基于梯度的方法共有的缺点。因此，我们引入了深度学习重要特征算法（DeepLIFT）。它既不是基于梯度的，也不是基于扰动的。它是一种基于反向传播的方法！

在本节中，我们将将其与 IG 进行对比。像 IG 和 Shapley 值一样，DeepLIFT 是为了完整性而设计的，因此符合显著的数学性质。除此之外，像 IG 一样，DeepLIFT 也可以应用于各种深度学习架构，包括 CNN 和循环神经网络（RNN），使其适用于不同的用例。

DeepLIFT 通过使用“参考差异”的概念，将模型的输出预测分解为每个输入特征的贡献。它通过网络层反向传播这些贡献，为每个输入特征分配一个重要性分数。

更具体地说，像 IG 一样，它使用一个基线，该基线代表关于任何类别的信息。然而，它随后计算输入和基线之间每个神经元的激活差异，并通过网络反向传播这些差异，计算每个神经元对输出预测的贡献。然后我们为每个输入特征求和其贡献，以获得其重要性分数（归因）。

它相对于 IG 的优势如下：

• 基于参考的：与 IG 等基于梯度的方法不同，DeepLIFT 明确地将输入与参考输入进行比较，这使得归因更加可解释和有意义。

• 非线性交互：DeepLIFT 在计算归因时考虑了神经元之间的非线性交互。它通过考虑神经网络每一层的乘数（由于输入的变化而导致的输出的变化）来捕捉这些交互。

• 稳定性：DeepLIFT 比基于梯度的方法更稳定，因为它对输入的小变化不太敏感，提供了更一致的归因。因此，在 DeepLIFT 归因上使用 SmoothGrad 是不必要的，尽管对于基于梯度的方法来说强烈推荐。

总体而言，DeepLIFT 提供了一种更可解释、更稳定和更全面的归因方法，使其成为理解和解释深度学习模型的有价值工具。

接下来，我们将以类似的方式创建 DeepLIFT 归因图：

deeplift_maps = get_attribution_maps(attr.DeepLift, garbage_mdl,\
                                     device, X_misclass, y_misclass) 

要绘制一个归因图，使用的代码几乎与 Grad-CAM 相同，只是将gradcam_maps替换为deeplift_maps。输出在图 7.18中展示。

图 7.18：将塑料误分类为生物废物的 DeepLIFT 热图

图 7.18的归因不如 IG 那样嘈杂。但它们似乎也聚集在阴影中的一些单调的黄色和深色区域；它还指向右上角附近的一些单调的绿色。

将所有这些结合起来

现在，我们将运用我们所学到的关于基于梯度的归因方法的一切知识，来理解所有选择的错误分类（塑料的假阴性金属的假阳性）的原因。正如我们处理中间激活图一样，我们可以利用compare_img_pred_viz函数将高分辨率的样本图像与四个归因图并排显示：显著性、Grad-CAM、SmoothGrad IG 和 DeepLift。为此，我们首先必须迭代所有错误分类的位置和索引，并提取所有图。请注意，我们正在使用tensor_to_img函数中的overlay_bg来生成一个新的图像，每个图像都叠加了原始图像和热图。最后，我们将四个归因输出连接成一个单独的图像（viz_img）。正如我们之前所做的那样，我们提取实际的标签（y_true）、预测标签（y_pred）和带有概率的pandas系列（probs_s），以便为我们将生成的图表添加一些上下文。for循环将生成六个图表，但为了简洁起见，我们只将讨论其中的三个：

for pos, idx in enumerate(misclass_idxs):
    orig_img = mldatasets.tensor_to_img(test_400_data[idx][0],\
                                   norm_std, norm_mean, to_numpy=True)
    bg_img = mldatasets.tensor_to_img(test_data[idx][0],\
                                   norm_std, norm_mean, to_numpy=True)
    map1 = mldatasets.tensor_to_img(
        saliency_maps[pos], to_numpy=True,\
        cmap_norm='positive', overlay_bg=bg_img
    )
    map2 = mldatasets.tensor_to_img(
        smooth_ig_maps[pos],to_numpy=True,\
        cmap_norm='positive', overlay_bg=bg_img
    )
    map3 = mldatasets.tensor_to_img(
        gradcam_maps[pos], to_numpy=True,\
        cmap_norm='positive', overlay_bg=bg_img
    )
    map4 = mldatasets.tensor_to_img(
        deeplift_maps[pos], to_numpy=True,\
        cmap_norm='positive', overlay_bg=bg_img
    )
    viz_img = cv2.vconcat([
        cv2.hconcat([map1, map2]),
        cv2.hconcat([map3, map4])
    ])
    label = int(y_test[idx])
    y_true = labels_l[label]
    y_pred = y_test_pred[idx]
    probs_s = probs_df.loc[idx]
    title = 'Gradient-Based Attr for Misclassification Sample #{}'.\
                                                           format(idx)
    mldatasets.compare_img_pred_viz(orig_img, viz_img, y_true,\
                                    y_pred, probs_s, title=title) 

之前的代码生成了图 7.19到图 7.21。重要的是要注意，在所有生成的图表中，我们都可以观察到左上角的显著性归因、右上角的 SmoothGrad IG、左下角的引导 Grad-CAM 和右下角的 DeepLIFT：

图 7.19：金属误分类为电池的基于梯度的归因 #8

在图 7.19中，所有四种归因方法之间缺乏一致性。显著性归因图显示，所有电池的中心部分都被视为金属表面，除了纸箱的白色部分。另一方面，SmoothGrad IG 主要聚焦于白色纸箱，而 Grad-CAM 几乎完全聚焦于蓝色纸箱。最后，DeepLIFT 的归因非常稀疏，仅指向白色纸箱的一些部分。

在图 7.20中，归因比图 7.19中的一致性要好得多。哑光白色区域明显让模型感到困惑。考虑到训练数据中的塑料主要是空塑料容器的单个部件——包括白色牛奶壶——这是有道理的。然而，人们确实回收玩具、塑料工具如勺子和其他塑料物品。有趣的是，尽管所有归因方法都在白色和浅黄色表面上都很显著，SmoothGrad IG 还突出了某些边缘，如一只鸭子的帽子和另一只的领子：

图 7.20：塑料误分类的基于梯度的归因 #86

继续探讨回收玩具的主题，乐高积木是如何被错误分类为电池的？参见图 7.21以获取解释：

图 7.21：塑料误分类的基于梯度的归因 #89

图 7.21展示了在所有归因方法中，主要是黄色和绿色的积木（以及较少的浅蓝色）是误分类的罪魁祸首，因为这些颜色在电池制造商中很受欢迎，正如训练数据所证明的那样。此外， studs 之间的平面表面获得了最多的归因，因为这些表面与电池的接触相似，尤其是 9 伏方形电池。与其他示例一样，显著性是最嘈杂的方法。然而，这次，引导 Grad-CAM 是最不嘈杂的。它也比其他方法在边缘上的显著性更强，而不是在表面上。

我们接下来将尝试通过在真实正例上执行的基于扰动的归因方法，来发现模型关于电池（除了白色玻璃之外）学到了什么。

通过扰动归因方法理解分类

到目前为止，这本书已经对基于扰动的方 法进行了大量的介绍。因此，我们介绍的大多数方法，包括 SHAP、LIME、锚点，甚至排列特征重要性，都采用了基于扰动的策略。这些策略背后的直觉是，如果你从你的输入数据中删除、更改或屏蔽特征，然后使用它们进行预测，你将能够将新预测与原始预测之间的差异归因于你在输入中做出的更改。这些策略可以在全局和局部解释方法中加以利用。

我们现在将像对错误分类样本所做的那样做，但针对选定的真阳性，并在单个张量（X_correctcls）中收集每个类别的四个样本：

correctcls_idxs = wglass_TP_idxs[:4] + battery_TP_idxs[:4] 
correctcls_data = torch.utils.data.Subset(test_data, correctcls_idxs)
correctcls_loader = torch.utils.data.DataLoader(correctcls_data,\
                                                batch_size = 32)
X_correctcls, y_correctcls = next(iter(correctcls_loader))
X_correctcls, y_correctcls = X_correctcls.to(device),\
                             y_correctcls.to(device) 

在图像上执行排列方法的一个更复杂方面是，不仅有几十个特征，而是有成千上万个特征需要排列。想象一下：224 x 224 等于 50,176 像素，如果我们想测量每个像素独立变化对结果的影响，我们至少需要为每个像素制作 20 个排列样本。所以，超过一百万！因此，几个排列方法接受掩码来确定一次要排列哪些像素块。如果我们将它们分成 32 x 32 像素的块，这意味着我们总共只有 49 个块需要排列。然而，尽管这会加快归因方法的速度，但如果我们块越大，就会错过对较小像素集的影响。

我们可以使用许多方法来创建掩码，例如使用分割算法根据表面和边缘将图像分割成直观的块。分割是按图像进行的，因此段的数量和位置将在图像之间变化。scikit-learn 的图像分割库（skimage.segmentation）有许多方法：scikit-image.org/docs/stable/api/skimage.segmentation.html。然而，我们将保持简单，并使用以下代码为所有 224 x 224 图像创建一个掩码：

feature_mask = torch.zeros(3, 224, 224).int().to(device)
counter = 0
strides = 16
for row in range(0, 224, strides):
    for col in range(0, 224, strides):
        feature_mask[:, row:row+strides, col:col+strides] = counter
        counter += 1 

前面的代码所做的初始化一个与模型输入大小相同的零张量。将这个张量概念化为一个空图像会更简单。然后它沿着 16 像素宽和高的步长移动，从图像的左上角到右下角。在移动过程中，它使用counter设置连续数字的值。最终你得到一个所有值都填充了 0 到 195 之间数字的张量，如果你将其可视化为一幅图像，它将是一个从左上角的黑色到右下角浅灰色的对角渐变。重要的是要注意，具有相同值的每个块都被归因方法视为相同的像素。

在我们继续前进之前，让我们讨论一下基线。在 Captum 归因方法中，正如其他库的情况一样，默认基线是一个全零张量，当图像由介于 0 和 1 之间的浮点数组成时，这通常等同于一个黑色图像。然而，在我们的情况下，我们正在标准化我们的输入张量，这样模型就不会看到最小值为 0 但平均值为 0 的张量！因此，对于我们的垃圾模型，全零张量对应于中等灰色图像，而不是黑色图像。对于基于梯度的方法，灰色图像基线本身并没有固有的错误，因为很可能存在许多步骤介于它和输入图像之间。然而，基于扰动的方 法可能对基线过于接近输入图像特别敏感，因为如果你用基线替换输入图像的部分，模型将无法区分出来！

对于我们的垃圾模型的情况，一个黑色图像由张量-2.1179组成，因为我们对输入张量执行标准化操作之一是(x-0.485)/0.229，当x=0时，这恰好等于大约-2.1179。你还可以计算当x=1时的张量；它转换为白色图像的2.64。话虽如此，假设在我们的真实阳性样本中，至少有一个像素具有最低值，另一个具有最高值，是没有害处的，因此我们将只使用max()和min()来创建亮暗基线：

baseline_light = float(X_correctcls.max().detach().cpu())
baseline_dark = float(X_correctcls.min().detach().cpu()) 

我们将只对除了一个扰动方法之外的所有方法使用一个基线，但请随意切换它们。现在，让我们继续为每种方法创建归因图！

特征消除

特征消除是一种相对简单的方法。它所做的是通过用基线替换它来遮挡样本输入图像的一部分，默认情况下，基线为零。目标是通过对改变它的效果进行观察，了解每个输入特征（或特征组）在做出预测中的重要性。

这就是特征消除是如何工作的：

1. 获取原始预测：首先，获取模型对原始输入的预测。这作为比较扰动输入特征效果的基准。

2. 扰动输入特征：接下来，对于每个输入特征（或由特征掩码设置的特征组），它被替换为基线值。这创建了一个“消除”版本的输入。

3. 获取扰动输入的预测：计算消除输入的模型预测。

4. 计算归因：计算原始输入和消除输入之间模型预测的差异。这个差异归因于改变的特征，表明它在预测中的重要性。

特征消除是一种简单直观的方法，用于理解模型预测中输入特征的重要性。然而，它也有一些局限性。它假设特征是独立的，可能无法准确捕捉特征之间交互的影响。此外，对于具有大量输入特征或复杂输入结构的模型，它可能计算成本高昂。尽管存在这些局限性，特征消除仍然是理解和解释模型行为的一个有价值的工具。

要生成归因图，我们将使用之前使用的get_attribution_maps函数，并输入额外的feature_mask和baselines参数：

ablation_maps = get_attribution_maps(
    attr.FeatureAblation,garbage_mdl,\
    device,X_correctcls,y_correctcls,\
    feature_mask=feature_mask,\
    baselines=baseline_dark
) 

要绘制归因图的示例，你可以复制我们用于显著性的相同代码，只是将saliency_maps替换为ablation_maps，并且我们使用occlusion_maps数组中的第二个图像，如下所示：

pos = 2
orig_img = mldatasets.tensor_to_img(X_correctcls[pos], norm_std,\
                                    norm_mean, to_numpy=True)
attr_map = mldatasets.tensor_to_img(occlusion_maps[pos],to_numpy=True,\
                         cmap_norm='positive') 
Figure 7.22:

图 7.22：测试数据集中白色玻璃真阳性的特征消除图

在图 7.22中，酒杯底部的特征组似乎是最重要的，因为它们的缺失对结果的影响最大，但酒杯的其他部分也有一定程度的显著性，除了酒杯的茎。这是有道理的，因为没有茎的酒杯仍然是一个类似玻璃的容器。

接下来，我们将讨论一种类似的方法，它将能够以更详细的方式展示归因。

遮挡敏感性

遮挡敏感性与特征消除非常相似，因为它也用基线替换了图像的部分。然而，与特征消除不同，它不使用特征掩码来分组像素。相反，它使用滑动窗口和步长自动将连续特征分组，在这个过程中，它创建了多个重叠区域。当这种情况发生时，它会对输出差异进行平均，以计算每个像素的归因。

在这个场景中，除了重叠区域及其对应平均值之外，遮挡敏感性和特征消除是相同的。事实上，如果我们使用滑动窗口和 3 x 16 x 16 的步长，就不会有任何重叠区域，特征分组将与由 16 x 16 块组成的feature_mask定义的特征分组相同。

那么，你可能想知道，熟悉这两种方法有什么意义？意义在于遮挡敏感性仅在固定分组连续特征很重要时才适用，比如图像和可能的其他空间数据。由于其使用步长，它可以捕捉特征之间的局部依赖性和空间关系。然而，尽管我们使用了连续的特征块，特征消融不必如此，因为feature_mask可以以任何对输入进行分段最有意义的方式排列。这个细节使其非常适用于其他数据类型。因此，特征消融是一种更通用的方法，可以处理各种输入类型和模型架构，而遮挡敏感性则是专门针对图像数据和卷积神经网络定制的，重点关注特征之间的空间关系。

要生成遮挡的归因图，我们将像以前一样操作，并输入额外的参数baselines、sliding_window_shapes和strides：

occlusion_maps = get_attribution_maps(
    attr.Occlusion, garbage_mdl,\
    device,X_correctcls,y_correctcls,\
    baselines=baseline_dark,\
    sliding_window_shapes=(3,16,16),\
    strides=(3,8,8)
) 

请注意，我们通过将步长设置为仅 8 像素，而滑动窗口为 16 像素，创建了充足的重叠区域。要绘制归因图，你可以复制我们用于特征消融的相同代码，只是将ablation_maps替换为occlusion_maps。输出如图图 7.23所示：

图 7.23：测试数据集中白色玻璃真阳性的遮挡敏感性图

通过图 7.23，我们可以看出遮挡的归因与消融的归因惊人地相似，只是分辨率更高。考虑到前者的特征掩码与后者的滑动窗口对齐，这种相似性并不令人惊讶。

无论我们使用 16 x 16 像素的非重叠块还是 8 x 8 像素的重叠块，它们的缺失影响都是独立测量的，以创建归因。因此，消融和遮挡方法都没有装备来测量非连续特征组之间的交互。当两个非连续特征组的缺失导致分类发生变化时，这可能会成为一个问题。例如，没有把手或底座的酒杯还能被认为是酒杯吗？当然可以被认为是玻璃，人们希望如此，但也许模型学到了错误的关系。

说到关系，接下来，我们将回顾一个老朋友：Shapley！

Shapley 值采样

如果你还记得第四章，全局模型无关解释方法，Shapley 提供了一种非常擅长衡量和归因特征联盟对结果影响的方法。Shapley 通过一次对整个特征联盟进行排列，而不是像前两种方法那样一次排列一个特征，来实现这一点。这样，它可以揭示多个特征或特征组如何相互作用。

创建归因图的代码现在应该非常熟悉了。这种方法使用feature_mask和baselines，但也测试了特征排列的数量（n_samples）。这个最后的属性对方法的保真度有巨大影响。然而，它可能会使计算成本变得非常昂贵，所以我们不会使用默认的每个排列 25 个样本来运行它。相反，我们将使用 5 个样本来使事情更易于管理。然而，如果你的硬件能够处理，请随意调整它：

svs_maps = get_attribution_maps(
    attr.ShapleyValueSampling,garbage_mdl,\
    device, X_correctcls, y_correctcls,\
    baselines=baseline_dark,\
    n_samples=5, feature_mask=feature_mask
) 
occlusion_maps is replaced by svs_maps. The output is shown in *Figure 7.24*:

图 7.24：测试数据集中白色玻璃真阳性的 Shapley 值采样图

图 7.24显示了一些一致的归因，例如最显著的区域位于酒杯碗的左下角。此外，底部似乎比遮挡和消融方法更重要。

然而，这些归因比之前的要嘈杂得多。这部分的理由是因为我们没有使用足够数量的样本来覆盖所有特征和交互的组合，部分原因是因为交互的混乱性质。对于单个独立特征的归因集中在几个区域是有意义的，例如酒杯的碗。然而，交互可能依赖于图像的几个部分，例如酒杯的底部和边缘。它们可能只有在它们一起出现时才变得重要。更有趣的是背景的影响。例如，如果你移除背景的一部分，酒杯是否不再像酒杯？也许背景比你想象的更重要，尤其是在处理半透明材料时。

KernelSHAP

既然我们谈论到了 Shapley 值，那么让我们尝试一下来自第四章，全局模型无关解释方法中的KernelSHAP。它利用 LIME 来更高效地计算 Shapley 值。Captum 的实现与 SHAP 类似，但它使用的是线性回归而不是 Lasso，并且计算核的方式也不同。此外，对于 LIME 图像解释器，最好使用有意义的特征组（称为超像素）而不是我们在特征掩码中使用的连续块。同样的建议也适用于KernelSHAP。然而，为了这个练习的简单性，我们也将保持一致性，并与其他三种基于排列的方法进行比较。

我们现在将创建归因图，但这次我们将使用浅色基线和深色基线各做一个。因为KernelSHAP是对 Shapley 采样值的近似，并且计算成本不是很高，所以我们可以将n_samples设置为 300。然而，这并不一定能保证高保真度，因为KernelSHAP需要大量的样本来近似相对较少的样本可以用 Shapley 彻底做到的事情：

kshap_light_maps = get_attribution_maps(attr.KernelShap, garbage_mdl,\
                                  device, X_correctcls, y_correctcls,\
                                  baselines=baseline_light,\
                                  n_samples=300,\
                                  feature_mask=feature_mask)
kshap_dark_maps = get_attribution_maps(attr.KernelShap, garbage_mdl,\
                                  device, X_correctcls, y_correctcls,\
                                  baselines=baseline_dark,\
                                  n_samples=300,\
                                  feature_mask=feature_mask) 
svs_maps is replaced by kshap_light_maps, and we modify the code to plot the attributions with the dark baselines in the third position, like this:

axs[2].imshow(attr_dark_map)
axs[2].grid(None)
axs[2].set_title("Kernel Shap Dark Baseline Heatmap") 
Figure 7.25:

图 7.25：测试数据集中白色玻璃真正阳性的 KernelSHAP 图

图 7.25中的两个归因图在大多数情况下并不一致，更重要的是，与之前的归因不一致。有时，某些方法比其他方法更难，或者需要一些调整才能按预期工作。

将所有这些结合起来

现在，我们将利用关于基于扰动归因方法的所有知识，来理解所有选择的真正阳性分类（无论是白色玻璃还是电池）的原因。正如我们之前所做的那样，我们可以利用compare_img_pred_viz函数将高分辨率样本图像与四个归因图并排放置：特征消融、遮挡敏感性、Shapley 和KernelSHAP。首先，我们必须迭代所有分类的位置和索引，并提取所有图。请注意，我们正在使用overlay_bg来生成一个新的图像，该图像将每个归因的热图叠加到原始图像上，就像我们在基于梯度的部分所做的那样。最后，我们将四个归因输出连接成一个单独的图像（viz_img）。正如我们之前所做的那样，我们提取实际的标签（y_true）、预测标签（y_pred）和包含概率的pandas系列（probs_s），以便为我们将生成的图表添加一些上下文。for循环将生成六个图表，但我们只会讨论其中的两个：

for pos, idx in enumerate(correctcls_idxs):
    orig_img = mldatasets.tensor_to_img(test_400_data[idx][0],\
                                        norm_std, norm_mean,\
                                        to_numpy=True)
    bg_img = mldatasets.tensor_to_img(test_data[idx][0],\
                                      norm_std, norm_mean,\
                                      to_numpy=True)
    map1 = mldatasets.tensor_to_img(ablation_maps[pos],\
                                    to_numpy=True,\
                                    cmap_norm='positive',\
                                    overlay_bg=bg_img)
    map2 = mldatasets.tensor_to_img(svs_maps[pos], to_numpy=True,\
                                    cmap_norm='positive',\
                                    overlay_bg=bg_img)
    map3 = mldatasets.tensor_to_img(occlusion_maps[pos],\
                                    to_numpy=True,\
                                    cmap_norm='positive',\
                                    overlay_bg=bg_img)
    map4 = mldatasets.tensor_to_img(kshap_dark_maps[pos],\
                                    to_numpy=True,\
                                    cmap_norm='positive',\
                                    overlay_bg=bg_img)
    viz_img = cv2.vconcat([
            cv2.hconcat([map1, map2]),
            cv2.hconcat([map3, map4])
        ])
    label = int(y_test[idx])
    y_true = labels_l[label]
    y_pred = y_test_pred[idx]
    probs_s = probs_df.loc[idx]
    title = 'Pertubation-Based Attr for Correct classification #{}'.\
                format(idx)
    mldatasets.compare_img_pred_viz(orig_img, viz_img, y_true,\
                                    y_pred, probs_s, title=title) 
Figures 7.26 to *7.28*. For your reference, ablation is in the top-left corner and occlusion is at the bottom left. Then, Shapley is at the top right and KernelSHAP is at the bottom right.

总体来说，你可以看出消融和遮挡非常一致，而 Shapley 和KernelSHAP则不太一致。然而，Shapley 和KernelSHAP的共同之处在于归因更加分散。

在图 7.26中，所有归因方法都突出了文本，以及至少电池的左侧接触。这与图 7.28相似，那里的文本也被大量突出显示，以及顶部接触。这表明，对于电池，模型已经学会了文本和接触都很重要。至于白色玻璃，则不太清楚。图 7.27中的所有归因方法都指向破碎花瓶的一些边缘，但并不总是相同的边缘（除了消融和遮挡，它们是一致的）：

图 7.26：电池分类的第 1 个基于扰动的归因

白色玻璃是三种玻璃中最难分类的，原因也不难理解：

图 7.27：基于扰动的白色玻璃分类#113 的归因

如图 7.27和其他测试示例中所述，模型很难区分白色玻璃和浅色背景。它设法用这些例子正确分类。然而，这并不意味着它在其他例子中也能很好地泛化，例如玻璃碎片和照明不足的情况。只要归因显示背景有显著的影响，就很难相信它仅凭镜面高光、纹理和边缘就能识别玻璃。

图 7.28：基于扰动的电池分类#2 的归因

对于图 7.28，在所有归因图中背景也被显著突出。但也许这是因为基线是暗的，整个物体也是如此。如果你用黑色方块替换电池边缘外的区域，模型会感到困惑是有道理的。因此，在使用基于排列的方法时，选择一个合适的基线非常重要。

任务完成

任务是提供一个对市政回收厂垃圾分类模型的客观评估。在样本外验证图像上的预测性能非常糟糕！你本可以就此停止，但那样你就不知道如何制作一个更好的模型。

然而，预测性能评估对于推导出特定的误分类以及正确的分类，以评估使用其他解释方法至关重要。为此，你运行了一系列的解释方法，包括激活、梯度、扰动和基于反向传播的方法。所有方法的一致意见是模型存在以下问题：

• 区分背景和物体

• 理解不同物体共享相似的颜色色调

• 混乱的照明条件，例如像酒杯那样的特定材料特性产生的镜面高光

• 无法区分每个物体的独特特征，例如乐高砖块中的塑料凸起与电池接触

• 被由多种材料组成的物体所困惑，例如塑料包装和纸盒包装中的电池

为了解决这些问题，模型需要用更多样化的数据集进行训练——希望是一个反映回收厂真实世界条件的数据集；例如，预期的背景（在输送带上）、不同的照明条件，甚至被手、手套、袋子等部分遮挡的物体。此外，他们应该为由多种材料组成的杂项物体添加一个类别。

一旦这个数据集被编译，利用数据增强使模型对各种变化（角度、亮度、对比度、饱和度和色调变化）更加鲁棒是至关重要的。他们甚至不需要从头开始重新训练模型！他们甚至可以微调 EfficientNet！

摘要

阅读本章后，你应该了解如何利用传统的解释方法来更全面地评估 CNN 分类器的预测性能，并使用基于激活的方法可视化 CNN 的学习过程。你还应该了解如何使用基于梯度和扰动的方法比较和对比误分类和真实正例。在下一章中，我们将研究 NLP 变换器的解释方法。

进一步阅读

• Smilkov, D., Thorat, N., Kim, B., Viégas, F., and Wattenberg, M., 2017, SmoothGrad: 通过添加噪声去除噪声. ArXiv, abs/1706.03825: arxiv.org/abs/1706.03825

• Sundararajan, M., Taly, A., and Yan, Q., 2017, 深度网络的公理化归因. 机器学习研究会议论文集，第 3319–3328 页，国际会议中心，悉尼，澳大利亚: arxiv.org/abs/1703.01365

• Zeiler, M.D., and Fergus, R., 2014, 视觉化和理解卷积网络. 在欧洲计算机视觉会议，第 818–833 页: arxiv.org/abs/1311.2901

• Shrikumar, A., Greenside, P., and Kundaje, A., 2017, 通过传播激活差异学习重要特征: arxiv.org/abs/1704.02685

在 Discord 上了解更多

要加入本书的 Discord 社区——在那里你可以分享反馈、向作者提问，并了解新版本——请扫描下面的二维码：

packt.link/inml

第八章：解释 NLP Transformer

在上一章中，我们学习了如何将解释方法应用于特定类型的深度学习模型架构，即卷积神经网络。在本章中，我们将提供一些工具来对 Transformer 模型架构执行相同的操作。Transformer 模型越来越受欢迎，它们最常见的使用案例是自然语言处理（NLP）。我们在第五章，局部模型无关解释方法中提到了 NLP。在本章中，我们也将这样做，但使用 Transformer 特定的方法和工具。首先，我们将讨论如何可视化注意力机制，然后是解释集成梯度属性，最后是探索瑞士军刀般的学习可解释性工具（LIT）。

我们将涵盖的主要主题包括：

• 使用 BertViz 可视化注意力

• 使用集成梯度解释标记属性

• LIME、反事实和其他 LIT 的可能性

技术要求

本章的示例使用了mldatasets、pandas、numpy、torch、transformers、bertviz、captum和lit-nlp库。如何安装所有这些库的说明在序言中。

本章的代码位于此处：packt.link/Yzf2L

使命

你是一名在纽约市一家即将推出的初创公司工作的数据科学家。这家初创公司旨在成为寻找城市中最佳、最新和最激动人心的美食目的地的首选之地！

目标是超越关于餐厅的典型结构化数据，深入挖掘网络上丰富的文本数据，从社交媒体网站到目录网站。初创公司认为，虽然评分可能提供对体验的简单量化，但评论包含更丰富的细节，可以提供多维度见解，了解什么使餐厅变得特别。

评论表达详细的情感，捕捉多样化的用户体验，与提供单一、非比较视角的评分不同。通过利用评论中的粒度，初创公司可以更精确地定制其推荐，以满足各种受众群体。

你们团队一直在讨论如何利用评论中的情感分析来确定如何最好地寻找体现用户在推荐系统中寻求体验感受的情感。二元情感分析（正面/负面）无法提供区分通常和独特体验，或针对特定群体（如旅行者、家庭或情侣）所需的细微差别。此外，初创公司的创始人认为用餐体验是多方面的。一种可能被视为“正面”的体验可能从“温馨而怀旧”到“刺激而冒险”。区分这些细微差别将使推荐系统能够更加个性化且有效。

你的经理遇到了一个情感分类模型，该模型使用名为 GoEmotions 的数据集进行训练，该数据集由谷歌发布。GoEmotions 提供了一种更详细的情感分类，比二元分类模型更有效地捕捉人类情感的丰富性。然而，首席策略师决定分类太多，决定将它们分组到另一个名为埃克曼的不同情绪分类法中（参见图 8.1）：

图 8.1：情绪分类法

埃克曼的情绪分类法是由心理学家保罗·埃克曼开发的一种分类系统，它识别出六种基本情绪，他认为这些情绪在所有人类文化中都是普遍体验到的。这六种情绪是快乐、悲伤、恐惧、愤怒、厌恶和惊讶。埃克曼提出，这些是基本情绪，它们被硬编码在我们的大脑中，并且世界各地的人们都以相同的方式表达，无论他们的文化如何。这些情绪可以通过特定的面部表情来识别，理解它们可以帮助心理学、沟通和社会学等领域。埃克曼的分类法提供了一套更简洁的情绪类别，同时仍然保留了细微差别。这使得对情绪的解释对开发团队和其他利益相关者来说更加可管理和可操作。然而，我们不得不保持中立，这不能归类到任何埃克曼类别中。

现在，下一步是使用 Tripadvisor 评论数据集来解释 GoEmotions Ekman 分类器模型，以了解模型学到了什么，并揭示可能对推荐系统开发有用的模式。这是一个开放性的任务。一般目标是理解模型在评论中识别出的模式以及这些模式如何与埃克曼的分类相对应。然而，这条路径可能导致许多发现或陷入死胡同。领导层强调，像你这样的数据科学家必须运用他们的判断力，在数据探索和模型解释中寻找机会。

通过揭示这些模式，初创公司可以微调其算法，寻找与这些情绪产生共鸣的评论。这些见解还可以指导餐厅合作、营销策略和平台功能增强。

方法

你决定采取三管齐下的方法：

1. 你将深入到转换器模型的内部，使用 BertViz 可视化注意力权重，以寻找这些机制中的相关模式。

2. 然后，你将生成显著性图，其中每个感兴趣的评论中的每个标记的归因都使用集成梯度法进行着色编码。

3. 最后，你将使用 LIT 来检验反事实情况。

你希望这些步骤能够向领导团队提供一些可操作的见解。

准备工作

你可以在这里找到这个示例的代码：github.com/PacktPublishing/Interpretable-Machine-Learning-with-Python-2E/tree/master/08/ReviewSentiment.ipynb

加载库

要运行此示例，您需要安装以下库：

• mldatasets 用于加载数据集

• pandas 和 numpy 用于操作

• torch（PyTorch）和 transformers 用于加载和配置模型

• bertviz、captum 和 lit-nlp 用于生成和可视化模型解释

您应该首先加载所有这些库：

import math
import os, random, re, gc
import warnings
warnings.filterwarnings("ignore")
import mldatasets
import numpy as np
import pandas as pd
import torch
from transformers import AutoTokenizer,\
                         AutoModelForSequenceClassification, pipeline
from bertviz import head_view, model_view
from captum.attr import LayerIntegratedGradients,\
                        TokenReferenceBase, visualization
from lit_nlp import notebook
from lit_nlp.api import dataset as lit_dataset
from lit_nlp.api import model as lit_model
from lit_nlp.api import types as lit_types 

接下来，我们进行数据处理和理解。

理解和准备数据

我们将数据这样加载到我们称为 reviews_df 的 DataFrame 中：

reviews_df = mldatasets.load("nyc-reviews", prepare=True) 

应该有超过 380,000 条记录和 12 列。我们可以使用 info() 来验证这一点：

reviews_df.info() 

输出检查无误。没有缺失值。然而，只有三个数值特征，一个日期，其余的都是对象数据类型，因为它们大多是文本。鉴于本章重点介绍 NLP，这并不令人惊讶。让我们检查数据字典，以了解我们将从该 DataFrame 中使用什么。

数据字典

这些是 DataFrame 中的 12 列，其中大部分是为了参考：

• review_id: ID – 评论的唯一标识符（仅作参考）

• author_id: ID – 作者的唯一标识符（仅作参考）

• restaurant_name: 文本 – 餐厅的名称（仅作参考）

• url_restaurant: URL – 统一资源标识符，用于定位包含餐厅评论的网页（仅作参考）

• review_date: 日期 – 评论制作的日期（仅作参考）

• review_title: 文本 – 作者为评论写的标题

• review_preview: 文本 – 为审查生成的预览

• review_full: 文本 – 作者撰写的完整评论

• rating: 序数 – 作者对场所给出的评级（1–5 级别）

• positive_sentiment: 二进制 – 根据二元情感模型，评论是否有积极情感（积极/消极）

• label: 分类 – GoEmotions 分类器预测的情绪（根据 Ekman 七类分类：快乐、中性、悲伤、厌恶、恐惧、愤怒和惊讶）

• score: 连续 – 预测评论属于预测类别的概率

这是一个多分类模型，因此它为每个类别预测了分数。然而，我们只存储了最可能类别的 score（标签）。因此，最后两列代表模型的输出。至于输入，让我们检查前三行来说明：

reviews_df[["review_title","review_full","label","score"]].head(3) 
Figure 8.2:

图 8.2：数据集的前三个评论

在图 8.2的前两列中，review_title和review_full代表模型的输入。它将它们作为一段单独的文本，因此当我们讨论评论时，我们指的是将这两个字符串通过冒号和空格连接起来的字符串，如下所示：令人失望：食物最多只是平庸。羊排是他们网站上首页展示的图片。

但这些并不是分析中唯一可能重要的列。我们当然可以按作者、餐厅、日期等分析评论，甚至将餐厅与地图上的特定坐标连接起来，以了解情感在地理上的变化。这一切都非常有趣。然而，我们不会在这里深入探讨，因为尽管这可能对情感分析的一般任务相关，但它会偏离本章的技术主题，即解释 Transformer 模型。

尽管如此，我们将探索一些与模型结果高度相关的特征，即作者提供的评分和二元情感分析模型的正面情感结果。你肯定期望这些会匹配，因为评论通常与评分一致——也就是说，正面评论的评分会比负面评论高。同样，一些情绪比负面情绪更积极。

为了更好地理解这些相关性，让我们将评论汇总起来，为每种情绪计算平均评分和正面情感，如下所示：

sum_cols_l = ["score","positive_sentiment","rating"]
summary_df = reviews_df.groupby("label")[sum_cols_l].agg(
    {"score":["count","mean"], "positive_sentiment":"mean", "rating":"mean"}
)
summary_df.columns = ["count", "avg. score", "% positive", "avg. rating"]
summary_df.sort_values(by="avg. rating", ascending=False).style.format(
    {
        "count":"{:,}",
        "avg. score":"{:.1%}",
        "% positive":"{:.1%}" ,
        "avg. rating":"{:.2f}"
    }
).bar(subset=["avg. score", "% positive", "avg. rating"],\
              color="#4EF", width=60) 

上述代码将在图 8.3中生成输出：

图 8.3：预测的评论数据集情绪的汇总表

正如你在图 8.3中可以看到的，380,000 条评论中的大多数都被放在了“喜悦”标签或类别中。喜悦是一种积极的情绪，因此我们的二元情感分类器将超过 90%的它们分类为正面，并且喜悦评论的平均评分几乎为 4.5。DataFrame 按平均评分排序，因为它不是模型（可能出错）的结果，所以它可能给我们提供最清晰的预测情绪的指示，即最终用户认为最积极的情绪。随着你向下查看列表，首先是积极的情绪，然后是中性的，最后是负面的。请注意，二元分类器认为正面评论的百分比与平均评分提供的相同顺序大致一致。

另一方面，每个标签的平均分数告诉我们平均来说预测有多自信它属于该标签。快乐、悲伤和惊讶是最自信的。由于多类预测是七个加起来等于 1 的数字，愤怒的平均分数为 56.6%表明，在愤怒是最可能情绪的许多预测中，其他情绪也可能有很大的概率——甚至可能是看似不兼容的情绪。我们将对此进行标记，因为这将会很有趣，以后可以探索这个问题。

你可以对图 8.3做出另一个有趣的解释，即尽管惊讶通常被认为是一种积极的情绪，但其中很大一部分是负面的。此外，平均评分低于 4，可能有很多负面评分拖累了它们。我们不会在本章中探索数据，但确实有很多表达惊讶情绪的负面评论。鉴于这一发现，为了适应任务，假设你向老板展示了这一点，他们决定关注惊讶是有道理的，因为市场研究显示人们喜欢发现和被“隐藏的宝藏”所惊喜。因此，一个推荐引擎能够帮助挖掘任何令人惊喜的餐厅，同时抑制那些持续令人失望的餐厅，这是至关重要的。

加载模型

之后，我们将从我们的数据集中随机选择，为了保持一致性，最好设置一个随机种子。在所有相关的库中初始化种子总是好的做法，尽管在这种情况下，它对 PyTorch 推理操作没有影响：

rand = 42
os.environ["PYTHONHASHSEED"]=str(rand)
random.seed(rand)
np.random.seed(rand)
torch.manual_seed(rand) 

接下来，让我们定义一个device变量，因为如果你有一个 CUDA 支持的 GPU，模型推理将执行得更快。然后，我们将使用from_pretrained函数从 Hugging Face 加载分词器（goemotions_tok）和模型（goemotions_mdl）。最后，我们将使用model.to(device)函数将所有权重和偏差移动到你的设备上，并使用model.eval()将模型设置为评估模式：

device = torch.device("cuda:0" if torch.cuda.is_available() else "cpu")
goemotions_mdl_path = "monologg/bert-base-cased-goemotions-ekman"
goemotions_tok = AutoTokenizer.from_pretrained(goemotions_mdl_path)
goemotions_mdl = AutoModelForSequenceClassification.from_pretrained(
    goemotions_mdl_path, output_attentions=True
)
goemotions_mdl.to(device)
goemotions_mdl.eval() 

模型加载后，我们总是使用print(goemotions_mdl)检查其架构。然而，在解释 Transformer 模型时，从广义上讲，最重要的是它们有多少层和注意力头。我们可以使用以下代码片段轻松检查：

num_layers = goemotions_mdl.config.num_hidden_layers
num_attention_heads = goemotions_mdl.config.num_attention_heads
print(f"The model has {num_layers} layers.")
print(f"Each layer has {num_attention_heads} attention heads.") 

应该说明的是，有 12 层和 12 个注意力头。在下一节中，我们将深入了解注意力机制的工作原理，以及如何使用 BertViz 可视化层和头：

使用 BertViz 可视化注意力

注意力是什么？让我们想象你正在阅读一本书，遇到一个提到你之前读过的角色的句子，但你忘记了关于他们的某些细节。你不太可能从头开始重新阅读，而是可能会快速浏览前面的页面，专注于具体讨论这个角色的部分。你的大脑会给予相关信息“注意力”，同时过滤掉不那么相关的部分。

变换器模型（如 BERT）中的注意力机制以类似的方式工作。在处理信息时，它不会平等地对待所有数据片段。相反，它“关注”最相关的部分，在当前任务的上下文中给予它们更多的重视。这种选择性地关注特定部分的能力有助于模型理解数据中的复杂模式和关系。

变换器由两个主要组件组成：编码器和解码器。每个组件都利用注意力机制，但它们以不同的方式使用：

• 编码器: 编码器的任务是理解输入数据。它通过使用注意力机制来确定输入的每个部分（如句子中的一个单词）如何与其他所有部分相关联。这使得编码器能够创建一个丰富的输入表示，捕捉其中的关系和上下文。这就像阅读一个句子并理解在这个句子中每个单词的含义。

• 解码器: 一旦编码器创建了这种表示，解码器就会使用它来生成输出。解码器还使用注意力机制，但以两种方式使用。首先，它关注编码器的表示来理解输入是什么。其次，它关注自己的先前输出，以确保当前输出与其迄今为止产生的输出一致。这就像根据你所读的和已经写下的内容来写一个有意义的句子。

然而，并非所有变换器模型都具有这两个组件。本质上，用例决定了变换器需要哪些部分：

• 编码器模型（如 BERT）：这些模型只使用变换器的编码器组件。它们通常用于涉及理解输入数据的任务，如情感分析（确定文本是积极的还是消极的）、命名实体识别（在文本中识别人、组织和地点）或其他分类任务。这是因为编码器的任务是创建一个表示输入数据的表示，捕捉其中的关系和上下文。

• 解码器模型（如 GPT、LLaMa 等）：这些用于涉及生成新数据的任务，如文本生成。变换器的解码器部分确保生成的输出与迄今为止产生的输出一致。

• 编码器-解码器模型（如 FLAN）：这些用于涉及将一块数据转换为另一块数据的任务，如翻译。编码器理解输入，解码器生成输出。

既然我们已经涵盖了注意力模型，那么让我们深入探讨 BERT。

BERT，代表来自转换器的双向编码器表示，是由 Google 开发的一种转换器模型。它用于理解和分析各种语言的文本数据。BERT 是一种读取文本双向以更好地理解单词上下文的转换器模型。它只使用转换器的编码部分，因为它的任务是理解文本，而不是生成文本。这使得 BERT 对于涉及理解文本的广泛任务非常有效。

因此，我们的 BERT 转换器模型有 12 层和 12 个注意力头。但这些是什么，它们是如何工作的呢？

• BERT 层：12 个隐藏层是 BERT 层，这些层是构成这个编码器转换器的其他层的堆叠。与我们在第七章中检查的 CNN 中的卷积层类似，可视化卷积神经网络，BERT 层代表抽象层。随着输入数据通过层传递，模型学习数据的越来越抽象的表示。在文本的上下文中，底层可能捕获基本句法信息，如单词在句子中的作用。随着你向上移动到层，它们倾向于捕获更高级别的语义，如整体句子意义或主题。层的数量，称为模型的深度，通常与其理解上下文和表示复杂关系的能力相关。然而，更多的层也需要更多的计算资源，并且如果没有足够的数据，可能会更容易过拟合。

• 注意力头：自注意力机制是任何转换器模型的核心。注意力头包含几个并行工作的自注意力机制。在多头自注意力机制内部，有多个独立的注意力头并行工作。每个注意力头都学会关注输入数据的不同部分（如不同标记之间的关系，这些通常是指单词）。拥有多个注意力头允许模型同时捕捉各种类型的关系。例如，一个头可能关注形容词和名词之间的关系，而另一个头可能捕捉动词-主语关系。在每个头计算其自己的注意力加权值表示之后，所有头的输出被连接起来并线性转换，以产生下一层的最终值表示。

让我们使用一些真实的评论来检查 GoEmotions 模型的内部工作原理。为此，我们将取四个样本评论，并使用以下代码打印它们的详细信息。在此过程中，我们将样本评论保存到字典（sample_reviews_dict）中，以便以后参考：

surprise_sample_reviews_l = [174067, 284154, 480395, 47659]
line_pattern = r"(?<=[.!?])\s+"
sample_reviews_dict = {}
for i, review_idx in enumerate(surprise_sample_reviews_l):
    review_s = reviews_df.loc[review_idx, :]
    sentiment = "Positive" if review_s["positive_sentiment"]\
                            else "Negative"
    review_lines_l = re.split(
        line_pattern, review_s["review_full"], maxsplit=1
    )
    review_txt = "\r\n\t\t".join(review_lines_l)
    print(f"{review_s["restaurant_name"]}") 
    print(f"\tSentiment:\t\t{sentiment}")
    print(f"\tRating:\t\t\t{review_s["rating"]}")
    print(f"\tGoEmotions Label:\t{review_s["label"]}")
    print(f"\tGoEmotions Score:\t{review_s["score"]:.1%}")
    print(f"\tTitle:\t{review_s["review_title"]}")
    print(f"\tReview:\t {review_txt}")
    sample_reviews_dict[i] = review_lines_l 
Figure 8.4:

图 8.4：数据集中的一些样本惊喜评论

正如你在图 8.4中看到的，所有评论样本的共同点是惊喜，无论是积极的还是消极的。

接下来，我们将利用 BertViz，尽管名字叫这个名字，但它可以可视化仅编码器 transformer 模型（如 BERT 及其所有变体）、仅解码器 transformer（如 GPT 及其所有变体）以及编码器-解码器 transformer（如 T5）的注意力。它非常灵活，但重要的是要注意，它是一个交互式工具，所以本节中用图表示的打印屏幕并不能完全体现其功能。

接下来，我们将创建一个函数，该函数使用分词器、模型和句子元组，可以创建两种不同的 BertViz 可视化：

def view_attention(tokenizer, model, sentences, view="model"):
    sentence_a, sentence_b = sentences
    # Encode sentences with tokenizer
    inputs = tokenizer.encode_plus(
        sentence_a, sentence_b, return_tensors="pt"
    )
    # Extract components from inputs
    input_ids = inputs["input_ids"]
    token_type_ids = inputs["token_type_ids"]
    # Get attention weights from model given the inputs
    attention = model(input_ids, token_type_ids=token_type_ids)[-1]
    # Get 2nd sentence start and tokens
    sentence_b_start = token_type_ids[0].tolist().index(1)
    input_id_list = input_ids[0].tolist()
    tokens = tokenizer.convert_ids_to_tokens(input_id_list)
    # BertViz visualizers
    if view=="head":
        head_view(attention, tokens, sentence_b_start)
    elif view=="model":
        model_view(attention, tokens, sentence_b_start) 

为了可视化注意力，我们需要取一对输入句子并使用我们的分词器（inputs）对它们进行编码。然后，我们提取这些输入的标记 ID（input_ids）和值，这些值表示每个标记属于哪个句子（token_type_ids）——换句话说，0代表第一句话，1代表第二句话。然后，我们将输入（input_ids和token_type_ids）传递给模型并提取attention权重。最后，有两个 BertViz 可视化器，head_view和model_view，为了使它们工作，我们只需要我们的输入产生的attention权重，将标记 ID 转换为tokens，以及第二句话开始的position（sentence_b_start）。

接下来，我们将使用model view在整个模型中可视化注意力。

使用模型视图绘制所有注意力

 sentences in the 1st sample review:

view_attention(
    goemotions_tok, goemotions_mdl, sample_reviews_dict[0], view="model"
) 

上述代码创建了一个大图，就像图 8.5中描绘的那样：

图 8.5：2nd Avenue Deli 的第一个样本评论的模型视图

图 8.5是一个 12 x 12 的网格，其中包含了 BERT 模型中的每个注意力头。我们可以点击任何注意力头来查看 BERT 输入中的两个句子，它们之间有线条相连，代表从左边的标记（一个句子）到右边的标记（另一个句子）的注意力权重。我们可以选择“Sentence A -> Sentence A”、“Sentence A -> Sentence B”以及所有介于两者之间的组合，以查看所有注意力权重的一个子集。权重接近一的线条看起来非常不透明，而接近零的权重则显示为透明，以至于几乎看不见。

一眼就能看出，有些注意力头有更多的线条、更粗的线条，或者看起来比另一个方向更明显的线条。我们可以点击单个注意力头来单独检查它们——例如：

• 层 1 头部 11 主要关注，从同一句子中的一个标记移动到下一个标记。这是一个非常常见的模式，并且完全符合逻辑，因为我们从左到右阅读英语，并且主要按照这个顺序理解它，尽管当然还有其他无疑在注意力头部中的模式。我们还看到了另一个常见模式的证据，即一个或多个标记对 [CLS] 标记具有注意力权重。[CLS] 标记是一个特殊的标记，在为分类任务使用类似 BERT 的模型时，会添加到每个输入序列的开头。它通常用于获取整个序列的聚合表示，以便进行分类。这意味着对于这个特定的注意力头部，标记在分类决策中起着作用。当注意力从分隔符标记 [SEP] 转移到分类标记 [CLS] 时，这可以被视为模型识别上下文或句子的结束，并反映其语义结论，从而可能影响分类决策。

• 层 9 头部 似乎执行一个更复杂的任务，即关联单词与“伟大”和“惊讶”，甚至跨越句子。这是另一个常见模式，其中连接的单词预测一个单词。

注意注意力头部中的模式，并注意一些其他模式，比如句子中的注意力向后移动，或者连接同义词。然后，将 sample_reviews_dict[0] 中的零改为一、二或三，看看注意力头部是否显示相同的模式。样本相当不同，但如果注意力头部没有做相同的事情，那么它们可能正在做非常相似的事情。然而，对于更大的图景，最好眯起眼睛，看看不同层中明显的模式。

接下来，我们将通过头部视图使这个任务变得更简单。

通过头部视图深入层注意力

我们可以通过以下代码开始选择第一个样本：

view_attention(
    goemotions_tok, goemotions_mdl, sample_reviews_dict[0], view="head"
) 

我们可以选择 12 层中的任何一层（0-11）。线条透明度意味着与模型视图相同的意思。然而，区别在于我们可以通过点击它们来隔离单个标记的注意力权重，因此当我们选择左侧的任何标记时，我们会看到线条与右侧的标记相连。此外，右侧将出现彩色编码的方框，表示每个 12 个注意力头部中每个的注意力权重。

如果我们恰好选择右侧的任何标记，我们将获得来自左侧标记的所有指向它的注意力。在 图 8.6 中，我们可以看到几个样本句子如何出现：

图 8.6：样本 1、2 和 4 的头部视图

在图 8.6中，看看我们如何像处理模型视图一样深入挖掘模式。例如，我们可以检查不同的注意力模式，从相邻的相关的单词（“伟大的三明治”和“尽管相当昂贵”）到仅在句子上下文中存在关系的那些（“三明治”和“昂贵”），甚至跨句子（“坏”和“被忽视的”，“厨师”和“菜单”）。

通过这样做，我们可以意识到可视化注意力头不仅仅是出于好奇，还可以帮助我们理解模型如何将单词之间的点连接起来，从而完成下游任务，如分类。也许这可以帮助我们了解下一步该做什么，无论是进一步微调模型以包含代表性不足的单词和情况，还是以不同的方式准备数据，以防止模型被特定的单词或一组单词所困惑。然而，考虑到注意力头的复杂性，寻找模型问题就像在 haystack 中找针一样。我们之所以在本章中从层和注意力头开始，是因为它提供了一种直观的方式来理解变压器如何编码标记之间的关系。

一个更好的开始方式是使用归因。归因方法是一种将计算输入的一部分对模型预测的贡献程度的方法。在第七章中，关于可视化卷积神经网络的图像中，我们计算归因的输入部分是像素。对于文本，等效的部分是标记，在这种情况下，由（主要是）单词组成，所以接下来，我们将生成标记归因。

使用集成梯度解释标记归因

集成梯度是一种流行的方法，在第七章中，我们解释并利用它为图像中的每个像素生成归因。该方法具有相同的步骤：

1. 选择一个基线输入：基线代表没有信息。对于图像，通常是一个纯黑色图像。对于文本，这可能是一个所有单词都被占位符如[PAD]替换的句子，或者只是一个空句子。

2. 逐步将这个基线输入变为实际的输入句子（例如，评论），一步一步地。在每一步中，你将基线稍微改变一点，使其接近实际输入。

3. 计算输出变化：对于每一步，计算模型预测的变化量。

4. 总结句子中每个单词的所有变化。这将为每个单词提供一个分数，表示它对模型最终预测的贡献程度。

然而，在我们能够使用集成梯度之前，最好定义一个将输入分词并在一步中执行模型推理的变压器管道：

goemotions = pipeline(
    model=goemotions_mdl,
    tokenizer=goemotions_tok,
    task="text-classification",
    function_to_apply="softmax",
    device=device,
    top_k=None
) 

你可以这样测试goemotions管道：

goemotions(["this restaurant was unexpectedly disgusting!",
            "this restaurant was shockingly amazing!"]) 

应该输出以下列表的列表字典：

[[{"label": "disgust", "score": 0.961812436580658},
  {"label": "surprise", "score": 0.022211072966456413},
  {"label": "sadness", "score": 0.004870257806032896},
  {"label": "anger", "score": 0.0034139526542276144},
  {"label": "joy", "score": 0.003016095608472824},
  {"label": "fear", "score": 0.0027414397336542606},
  {"label": "neutral", "score": 0.0019347501220181584}],
 [{"label": "joy", "score": 0.6631762385368347},
  {"label": "surprise", "score": 0.3326122760772705},
  {"label": "neutral", "score": 0.001732577453367412},
  {"label": "anger", "score": 0.0011324150254949927},
  {"label": "sadness", "score": 0.0010195496724918485},
  {"label": "fear", "score": 0.00021178492170292884},
  {"label": "disgust", "score": 0.00011514205834828317}]] 

如您所见，在第一个列表中，有两个预测（每个文本一个），每个预测都有一个包含七个字典的列表，其中一个字典包含每个类别的分数。由于字典是按最高分数到最低分数排序的，您可以判断第一个餐厅评论主要被预测为厌恶，第二个评论为喜悦，占 66%，但也有相当数量的惊讶。

接下来，我们将创建一个函数，该函数可以接受任何包含我们的评论和我们的转换器的 DataFrame 行，并为每个超过 10%概率的预测生成和输出归因：

def visualize_ig_review(interpret_s:pd.Series,
                          pline:pipeline,
                          max_prob_thresh:float=0.1,
                          max_classes=np.PINF,
                          concat_title=True,
                          summary_df=None
) -> pd.DataFrame:
    print(f"{interpret_s.name}: {interpret_s['restaurant_name']}")
    # Init some variables
    if concat_title:
        text = interpret_s["review_title"] + ":" + interpret_s["review_full"]
    else:
        text = interpret_s["review_full"]
    true_label = "Positive" if interpret_s["positive_sentiment"]\
                            else "Negative"
    rating = interpret_s["rating"]
    # Get predictions
    prediction = pline(text)[0]
    prediction_df = pd.DataFrame(prediction)
    if summary_df is not None:
        prediction_df["label_avg_rating"] = prediction_df.label.\
            replace(summary_df["avg. rating"].to_dict())
        prediction_df = prediction_df.sort_values("label_avg_rating",\
           ascending=False).reset_index(drop=True)
    # Process predictions
    prediction_tuples = [(p["label"], p["score"]) for p in prediction]
    sorted_prediction_tuples = sorted(prediction_tuples,\
        key=lambda x: x[1], reverse=True)
    pred_class, pred_prob = sorted_prediction_tuples[0]
    # Initialize Integrated Gradients
    forward_func = lambda inputs, position=0: pline.model(
        inputs, attention_mask=torch.ones_like(inputs)
    )[position]
    layer = getattr(pline.model, "bert").embeddings
    lig = LayerIntegratedGradients(forward_func, layer)
    # Prepare tokens and baseline
    device = torch.device("cuda:0" if torch.cuda.is_available()\
                            else "cpu")
    inputs = torch.tensor(pline.tokenizer.encode(text,\
        add_special_tokens=False), device = device).unsqueeze(0)
    tokens = pline.tokenizer.convert_ids_to_tokens(
        inputs.detach().numpy()[0]
    )
    sequence_len = inputs.shape[1]
    baseline = torch.tensor(
        [pline.tokenizer.cls_token_id]\
        + [pline.tokenizer.pad_token_id] * (sequence_len - 2)\
        + [pline.tokenizer.sep_token_id],\
        device=device
    ).unsqueeze(0)
    # Iterate over every prediction
    vis_record_l = []
    for i, (attr_class, attr_score) in\ 
        enumerate(sorted_prediction_tuples):
        if (attr_score > max_prob_thresh) and (i < max_classes):
            # Sets the target class
            target = pline.model.config.label2id[attr_class]
            # Get attributions
            with torch.no_grad():
                attributes, delta = lig.attribute(
                    inputs=inputs,
                    baselines=baseline,
                    target=target,
                    return_convergence_delta = True
                )
            # Post-processing attributions
            attr = attributes.sum(dim=2).squeeze(0)
            attr = attr / torch.norm(attr)
            attr = attr.cpu().detach().numpy()
            # Generate & Append Visualization Data Record
            vis_record = visualization.VisualizationDataRecord(
                    word_attributions=attr,
                    pred_prob=pred_prob,
                    pred_class=pred_class,
                    true_class=f"{true_label} ({rating})",
                    attr_class=attr_class,
                    attr_score=attr_score,
                    raw_input_ids=tokens,
                    convergence_score=delta
            )
            vis_record_l.append(vis_record)
    # Display list of visualization data records
    _ = visualization.visualize_text(vis_record_l)
    return prediction_df 

虽然代码量看起来很复杂，但单独解释时，有很多步骤相对简单。我们将从模型推理开始，逐步进行：

1. 获取预测结果：这是一个非常直接的步骤。它只需将text输入到管道（pline）中。它只取返回的第一个项目（[0]），因为它只预期输入并由此通过管道返回一个预测。接下来的几行显示了如果函数接收到sample_df，模型会做什么，它实际上只需要按平均最佳评分的顺序对预测进行排序。

2. 处理预测：在这里，代码确保预测被排序并放入元组中，以便在后续的for循环中更容易迭代。

3. 初始化集成梯度：定义了一个正向函数，它接受输入并返回给定位置的模型输出，以及一个用于计算归因的层，在这个例子中是嵌入层。然后，使用正向函数和指定的层初始化一个LayerIntegratedGradients（lig）实例。

4. 准备标记和基线：首先，对文本进行标记并转换为张量，然后将其移动到指定的设备。然后，将标记 ID 转换回tokens，以供潜在的视觉化或分析使用。为集成梯度方法创建一个baseline。它由[CLS][ token at the start, ]{custom-style="P - Code"}[SEP][ token at the end, and ]{custom-style="P - Code"}[PAD][ tokens in the middle matching the length of the input ]{custom-style="P - Code"}text组成。

5. 遍历每个预测：这是一个for循环，它遍历每个预测，只要概率超过由max_prob_threshold定义的 10%，我们将在循环中进行以下操作：

   1. 设置目标类别：集成梯度是一种有向归因方法，因此我们需要知道为哪个target类别生成归因；因此，我们需要模型内部用于预测类别的 ID。

   2. 获取归因：使用与第七章中使用的相同的 Captum attribute方法，我们为我们的文本（inputs）、基线、目标和决定是否返回 IG 方法的增量（一个近似误差的度量）生成 IG 归因。

   3. 后处理归因：IG 方法返回的归因形状为 (num_inputs, sequence_length, embedding_dim)，其中对于此模型，embedding_dim=768，sequence_length对应输入中的标记数量，而num_inputs=1，因为我们一次只执行一个归因。因此，每个标记的嵌入都有一个归因分数，但我们需要的每个标记一个归因。因此，这些分数在嵌入维度上求和，以得到序列中每个标记的单个归因值。然后，对归因进行归一化，确保归因的幅度在 0 到 1 之间，并且处于可比较的尺度。最后，将归因从计算图中分离出来，移动到 CPU，并转换为numpy数组以进行进一步处理或可视化。

   4. 生成并附加可视化数据记录：Captum 有一个名为VisualizationDataRecord的方法，用于创建每个归因的记录以供可视化，因此在这一步中，我们使用这些归因、增量、标记和与预测相关的元数据创建这些记录。然后，将此数据记录附加到列表中。

6. 显示可视化数据记录列表：利用visualize_text显示记录列表。

现在，让我们创建一些样本以执行集成梯度归因：

neg_surprise_df = reviews_df[
    (reviews_df["label"]=="surprise")
    & (reviews_df["score"]>0.9)
    & (reviews_df["positive_sentiment"]==0)
    & (reviews_df["rating"]<3)
] #43
neg_surprise_samp_df = neg_surprise_df.sample(
    n=10, random_state=rand
) 

在上述代码片段中，我们选取所有概率超过 90%的惊喜评论，但为了确保它们是负面的，我们将选择一个负面情感和低于三的评分。然后，我们将从这些评论中随机抽取 10 条。

接下来，我们将遍历列表中的每个评论并生成一些可视化。其他一些可视化显示在图 8.7的屏幕截图上：

for i in range(10):
    sample_to_interpret = neg_surprise_samp_df.iloc[i]
    _ = visualize_ig_review(
        sample_to_interpret, goemotions, concat_title=True, summary_df=summary_df
) 
Figure 8.7:

图 8.7：IG 可视化对于负面惊喜

如图 8.7所示，惊喜预测归因于“震惊”、“意识到”和“神秘”等单词，以及“不确定如何”等短语。这些都很有意义，因为它们表明某件事是未知的。自然地，也有一些情况，单词“surprise”或“surprised”就足以得到惊喜预测。然而，有时事情并不那么简单。在最后一个例子中，并不是一个单词似乎表明了惊喜，而是许多单词，大意是某件事不协调。更具体地说，这些来自伦敦的游客对纽约市的一家熟食店如此昂贵感到非常惊讶。请注意，“负面”和“正面”的颜色编码并不意味着一个单词是负面的或正面的，而是它对（负面）或有利于（正面）归因标签的权重。

接下来，我们将重复运行与生成样本惊喜负面评论的 IG 解释相似的代码，但这次是为了正面评论。为了确保它们是正面的，我们将使用 4 分以上的评分。这次，我们将确保从样本中移除包含“惊喜”一词的任何评论，以使事情更有趣：

pos_surprise_df = reviews_df[
    (reviews_df["label"]=="surprise")
    & (reviews_df["score"]>0.97)
    & (reviews_df["positive_sentiment"]==1)
    & (reviews_df["rating"]>4)
]
pos_surprise_samp_df = pos_surprise_df[
    ~pos_surprise_df["review_full"].str.contains("surprise")
]
for i in range(10):
    sample_to_interpret = pos_surprise_samp_df.iloc[i]
    _ = visualize_ig_review(
        sample_to_interpret, goemotions,\
        concat_title=False, summary_df=summary_df
) 

上一段代码将在图 8.8中生成可视化。

图 8.8：积极惊喜的 IG 可视化

图 8.8展示了“困惑”和“难以置信”等单词，以及“难以置信的是”这样的短语如何表明惊喜。还有一些情况，标记对惊喜预测有负面影响。例如，对于最后一家餐厅，“适合每个人”并不使其非常令人惊讶。此外，您会注意到中间的日本餐厅被预测为同时体现惊喜和快乐情绪。有趣的是，有些单词与一种情绪相关，但与另一种情绪不太相关，有时它们甚至表示相反的意思，比如“很难找到地方”中的“hard”表示惊喜但不表示快乐。

找到像日本餐厅这样的混合情感评论可能有助于解释为什么有些评论难以完全用单一情感进行分类。因此，现在我们将生成一些正负混合评论样本。我们可以通过确保预测标签的分数永远不会超过 50%来轻松做到这一点：

pos_mixed_samp_df = reviews_df[
    (~reviews_df["label"].isin(["neutral","joy"]))
    & (reviews_df["score"] < 0.5)
    & (reviews_df["positive_sentiment"]==1)
    & (reviews_df["rating"]< 5)
].sample(n=10, random_state=rand)
neg_mixed_samp_df = reviews_df[
    (~reviews_df["label"].isin(["neutral","joy"]))
    & (reviews_df["score"] < 0.5)
    & (reviews_df["positive_sentiment"]==0)
    & (reviews_df["rating"]>2)
].sample(n=10, random_state=rand) 
mldatasets.plot_polar, which plots a polar line chart for the predictions with plotly. You’ll need both plotly and kaleido to make this work:

for i in range(10):
    sample_to_interpret = pos_mixed_samp_df.iloc[i]
    prediction_df = visualize_ig_review(
        sample_to_interpret,\
        goemotions, concat_title=False,\
        summary_df=summary_df
    )
    rest_name = sample_to_interpret["restaurant_name"]
    mldatasets.plot_polar(
    prediction_df, "score", "label", name=rest_name
) 

上一段代码将在图 8.9中生成 IG 可视化和极线图：

图 8.9：混合情感评论的 IG 可视化

图 8.9展示了这家三明治店的评论似乎表现出快乐、恐惧和中性情绪。单词“terrific”和“friendly”与快乐相连，但与中性不相干。然而，奇怪的是，“terrific”这个词也与恐惧相关。也许这与对这个词进行的 WordPiece 分词有关。注意，“terrific”以三个子词标记出现，分别是 te、##rri 和##fic。这很可能是因为用于训练模型的原语料库（Reddit 评论）中“terrific”这个词的频率不够高，无法将其作为独立单词包含，但这些子词却包含了。这种技术的缺点是，可能“te”和“rri”标记经常用于像“terrifying”这样的单词，以及像“horrific”和“mortific”这样的其他可怕单词中的“fic”。另一方面，“fic”出现在“magnificent”和“beneficial”中。所以，尽管有上下文嵌入，子词标记可能会引起一些歧义。

我们现在可以运行与之前相同的代码，但针对neg_mixed_samp_df，以检查其他示例并得出自己的结论。接下来，我们可以使用 LIT 扩展我们的 XAI NLP 特定工具集。

LIT、反事实以及其他可能性

LIT 是一个由People+AI Research（PAIR）倡议开发的开源平台，用于可视化和理解 NLP 模型。PAIR 开发了第六章中特色展示的“如果工具”（What-If Tool）。

LIT 提供了一个交互式和可视化的界面，以便深入探究 NLP 模型的行为。使用 LIT，用户可以：

• 识别模型表现不佳的示例类型。

• 确定特定模型预测背后的原因。

• 测试模型在文本变化（如风格、动词时态或代词性别）下的一致性。

LIT 提供了各种内置功能，包括显著性图、注意力可视化、指标计算和反事实生成。然而，它也支持定制，允许添加专门的解释性技术、可视化等。

尽管 LIT 的主要关注点是文本语言数据，但它也支持在图像和表格数据上操作的模式。它与包括 TensorFlow 和 PyTorch 在内的多种机器学习框架兼容。该工具可以作为独立服务器运行，也可以在 Colab、Jupyter 和 Google Cloud Vertex AI 笔记本等笔记本环境中运行。

为了与任何自定义数据集一起工作，LIT 提供了一个Dataset子类来创建一个兼容 LIT 的数据集加载器。您必须包含一个__init__，它加载数据集，以及一个 spec 函数，它指定数据集中返回的数据类型，而lit_nlp.api.types提供了一种确保 LIT 识别您数据集中的每个特征的方法。在这种情况下，我们提供了评论（TextSegment），标签（CategoryLabel）以及七个标签，以及两个额外的类别，这些类别可用于切片和分箱：

class GEDataset(lit_dataset.Dataset):
    GE_LABELS = ["anger", "disgust", "fear", "joy",\
                 "neutral", "sadness", "surprise"]
    def __init__(self, df: pd.DataFrame):
        self._examples = [{
            "review": row["review_title"] + ":" + row["review_full"],
            "label": row["label"],
            "rating": row["rating"],
            "positive": row["positive_sentiment"]
        } for _, row in df.iterrows()]
    def spec(self):
        return {
            "review": lit_types.TextSegment(),
            "label": lit_types.CategoryLabel(vocab=self.GE_LABELS),
            "rating": lit_types.CategoryLabel(),
            "positive": lit_types.CategoryLabel()
        } 

为了使 LIT 能够适应任何模型，有一个Model子类来创建一个 LIT 兼容的模型加载器。它还需要__init__函数来初始化模型，以及一个predict_minibatch函数来使用它进行预测。为此，我们还需要为predict函数的输入（input_spec）和输出（output_spec）创建规范。在这种情况下，我们输入一个评论（类型为TextSegment），并返回类型为MulticlassPreds的概率。请注意，模型的输出并不总是一致的，因为每个预测都是从最高分到最低分排列的。注意，为了使predict_minibatch的输出符合MulticlassPreds，我们必须将概率作为与标签（GE_LABELS）相对应的列表排列，与提供给vocab的顺序相同：

class GEModel(lit_model.Model):
    GE_LABELS = ["anger", "disgust", "fear", "joy",\
                 "neutral", "sadness", "surprise"]
    def __init__(self, model, tokenizer, **kw):
        self._model = pipeline(
            model=model,
            tokenizer=tokenizer,
            task="text-classification",
            function_to_apply="softmax",
            device=device,
            top_k=None
        )
    def input_spec(self):
        return {
            "review": lit_types.TextSegment()
        }
    def output_spec(self):
        return {
            "probas": lit_types.MulticlassPreds(vocab=self.GE_LABELS,\
                                                parent="label")
        }
    def predict_minibatch(self, inputs):
        examples = [d["review"] for d in inputs]
        with torch.no_grad():
            preds = self._model(examples)
        preds = [{p["label"]:p["score"] for p in pred_dicts}\
                for pred_dicts in preds]
        preds = [dict(sorted(pred_dict.items())) for pred_dict in preds]
        preds = [{"probas": list(pred_dict.values())} for pred_dict in preds]
        return preds 

好的，现在我们拥有了 LIT 运行所需的两个类。GoEmotions 模型初始化器（GEModel）接收模型（goemotions_mdl）和分词器（goemotions_tok）。我们将这些放入字典中，因为 LIT 可以接受多个模型和多个数据集进行比较。对于数据集，为了使其快速加载，我们将使用 100 个样本（samples100_df），由我们迄今为止创建的四个 10 样本 DataFrame 组成，再加上从整个评论数据集中随机抽取的 60 个额外样本。然后，我们将我们的 100 样本 DataFrame 输入到 GoEmotions 数据集初始化器（GEDataset）中，并将其放置到我们的数据集字典中，命名为NYCRestaurants。最后，我们通过输入我们的模型和数据集字典以及render它来创建小部件（notebook.LitWidget）。请注意，如果您想在笔记本环境之外运行此代码，可以使用Server命令使其在 LIT 服务器上运行：

models = {"GoEmotion":GEModel(goemotions_mdl, goemotions_tok)}
samples100_df = pd.concat(
    [
        neg_surprise_samp_df,
        pos_surprise_samp_df,
        neg_mixed_samp_df,
        pos_mixed_samp_df,
        reviews_df.sample(n=60, random_state=rand)
    ]
)
datasets = {"NYCRestaurants":GEDataset(samples100_df)}
widget = notebook.LitWidget(models, datasets)
widget.render(height=600)
# litserver = lit_nlp.dev_server.Server(models, datasets, port=4321)
# litserver.serve() 
Figure 8.10:

图 8.10：打开预测标签的笔记本视图

正如您在图 8.10中可以看到的，LIT 具有：

• 一个顶部栏，有一个下拉菜单用于选择模型和数据集（但在这里您不能选择，因为每种只有一个）和三种不同的视图（简单、默认和笔记本）。笔记本是默认选项。

• 一个选择栏，用于选择数据点和查看哪些被固定。

• 一个标签栏，有三个标签（预测、解释和分析）。默认情况下，预测被选中，此标签页左侧有数据表，您可以在其中选择和固定单个数据点，右侧有分类结果面板。

尽管笔记本视图比简单视图有更多功能，但它缺少默认视图中许多可用的功能。从现在开始，我们将检查图 8.11中描述的默认视图：

图 8.11：打开预测标签的默认视图

如你在图 8.11中看到的，默认视图有两个永久窗格，顶窗格中有数据表和数据点编辑器，底窗格中有标签页。这不是一个小笔记本单元格的好布局，但它可以让你轻松固定、选择和编辑数据点，同时在下面的标签页上执行任务。注意，标签页不止三个。我们将简要解释每个标签页：

• 预测: 这让你可以看到所选和固定的数据点的分类结果。注意，它用“P”表示预测标签，用“T”表示真实标签。然而，由于我们没有用这个数据集训练模型，提供的标签与预测的标签没有区别，但这可以证明在检查错误分类时非常有用。在分类结果的右侧，我们有标量，这允许我们比较固定和所选数据点的分数与数据集中所有其他数据点的分数。

• 解释: 在这里，我们可以对我们的数据点使用多种解释/归因方法，例如 LIME 和集成梯度。

• 显著性聚类: 我们对许多数据点进行归因，并将结果聚类以了解标记是如何聚类的。鉴于我们只使用 100 个数据集，我们不会在这里详细介绍。

• 指标: 如果我们使用带有真实标签的训练数据集，这个标签页将非常有用，因为它可以以多种方式切片和分组性能指标。

• 反事实: 与第六章类似，这里的反事实概念是相同的，即找出你可以改变的特征（在这种情况下是一个标记），这样就可以修改模型结果（预测标签）。这里提供了几种反事实发现方法。

因此，我们将按顺序处理这个列表，排除显著性聚类和预测（我们已经在图 8.11中解释过），所以接下来，我们将查看解释，如图8.12所示：

图 8.12：在“解释”标签页中，比较了固定和选择的评论的 LIME 解释

图 8.12展示了 LIME 解释在固定数据和选定点之间的差异。LIME 之前在第五章，局部模型无关解释方法中，以及在 NLP 的背景下都有涉及。这里也是一样。顺便提一下，尽管至少有四种方法可用，包括积分梯度，但只有 LIME 能与这个模型一起工作。这是因为 LIME 是一种模型无关的基于排列的方法，它不需要访问模型的所有内在参数，而其他方法不是模型无关的。如果你还记得，我们的GEModel没有暴露任何内在参数。如果我们想在 LIT 中使用基于梯度的方法如 IG，我们就需要不使用管道，然后以这种方式指定输入和输出，即暴露标记嵌入。LIT 网站上有一些示例可以帮助你完成这项工作。

接下来，我们将查看指标选项卡，如图 8.13 所示：

图 8.13：指标选项卡中的混淆矩阵

在图 8.13中的指标面板中，通常会有关于整个数据集、你所做的选择以及你可能定义的任何附加维度的信息性指标。然而，如果你展开选项卡，你总是会看到这个数据集的 100%准确率，因为没有真实标签。也许右侧的混淆矩阵面板在这种情况下更有信息性，因为我们可以看到标签和评分，或标签和正面的交叉表，因为我们定义了评分和正为CategoryLabel。请注意，从技术上讲，这并不是一个混淆矩阵，因为它没有将预测情感标签与相应的真实标签进行比较，但你可以看到预测标签和评分之间有多少一致性。

最后，让我们检查反事实选项卡，如图 8.14 所示：

图 8.14：在反事实选项卡中生成消融翻转反事实

图 8.14中的反事实选项卡提供了几种方法来改变输入，从而修改预测标签。由于一些反事实方法是模型无关的，而其他方法需要内在参数，因此并非所有方法都适用于GEModel。在这里，我们使用模型无关的方法消融翻转，从输入中移除标记。消融翻转简单地尝试从输入中删除标记，以找出哪个标记改变了预测。正如你所看到的，第一个移除的是“Not”从“review”中，第二个移除的是“impressed”。

通过反事实，你可以测试在 te ##rri ##fic（如图 8.9 所示）中添加和删除子词标记是否会在许多不同的上下文中引起歧义。例如，你可以从被认为具有中性或负面情绪的评论中逐个移除 te ##rri ##fic 中的一个标记，看看预测是否向积极方向改变。你也可以用“magnificent”这样的同义词替换所有三个标记。

任务完成

从摘要表（图 8.3）中很明显，并且通过集成梯度和一定程度上的注意力可视化练习得到证实，许多 Ekman 情绪在评论中难以辨别，恐惧、厌恶、愤怒和悲伤产生了许多混合情绪的评论。而这些难以辨别的都是负面情绪。

此外，在 GoEmotions 和 Ekman 分类法中，许多情绪在推荐引擎的上下文中并不重要，因此考虑合并一些负面情绪是有意义的。鉴于惊喜类别有时是积极的，有时是消极的，因此将它们拆分以包括好奇心和困惑是有道理的。

另一个重要的发现是，鉴于你发现的许多一致模式，惊喜并不难分类，但却是预测的关键情绪。然而，有好的惊喜和坏的惊喜。并且，在适当的训练数据下，模型很可能以高精度区分两者。我们可以确保标记化永远不会将传达情感的单词分开。

摘要

阅读本章后，你应该了解如何利用 BertViz 可视化变压器模型、层和注意力头，以及如何使用 Captum 的归因方法，特别是集成梯度，以及可视化数据记录来查看哪些标记负责预测的标签。最后，你应该对如何开始使用 LIT 有一个扎实的掌握。在下一章中，我们将探讨解释多元时间序列模型。

进一步阅读

• Vig, J., 2019, Transformer 模型中的注意力多尺度可视化。ArXiv：arxiv.org/abs/1906.05714

• Kokhlikyan, N., Miglani, V., Martin, M., Wang, E., Alsallakh, B., Reynolds, J., Melnikov, A., Kliushkina, N., Araya, C., Yan, S., & Reblitz-Richardson, O., 2020, Captum：PyTorch 的一个统一和通用的模型可解释性库。ArXiv：arxiv.org/abs/2009.07896

• Tenney, I., Wexler, J., Bastings, J., Bolukbasi, T., Coenen, A., Gehrmann, S., Jiang, E., Pushkarna, M., Radebaugh, C., Reif, E., & Yuan, A., 2020, *《语言可解释性工具：NLP 模型的可扩展、交互式可视化和分析工具》。实证自然语言处理方法会议：arxiv.org/abs/2008.05122

在 Discord 上了解更多信息

要加入本书的 Discord 社区——在那里您可以分享反馈、向作者提问，并了解新书发布——请扫描下面的二维码：

packt.link/inml

第九章：多元预测和敏感性分析的解释方法

在整本书中，我们学习了我们可以用来解释监督学习模型的各种方法。它们在评估模型的同时，也能揭示其最有影响力的预测因子及其隐藏的相互作用。但是，正如“监督学习”这个术语所暗示的，这些方法只能利用已知的样本以及基于这些已知样本分布的排列。然而，当这些样本代表过去时，事情可能会变得复杂！正如诺贝尔物理学奖获得者尼尔斯·玻尔著名地打趣说，“预测是非常困难的，尤其是如果它关乎未来。”

事实上，当你看到时间序列中的数据点波动时，它们可能看起来像是在按照可预测的模式有节奏地跳舞——至少在最佳情况下是这样。就像一个舞者随着节奏移动，每一次重复的动作（或频率）都可以归因于季节性模式，而音量（或振幅）的逐渐变化则可以归因于同样可预测的趋势。这种舞蹈不可避免地具有误导性，因为总会有一些缺失的拼图碎片稍微改变数据点，比如供应商供应链中的延迟导致今天销售数据的意外下降。更糟糕的是，还有不可预见的、十年一遇、一代人一遇或甚至一次性的灾难性事件，这些事件可以彻底改变人们对时间序列运动的一些理解，就像一个舞厅舞者在抽搐。例如，在 2020 年，由于 COVID-19，无论好坏，各地的销售预测都变得毫无用处！

我们可以将这种情况称为极端异常事件，但我们必须认识到，模型并不是为了预测这些重大事件而构建的，因为它们几乎完全基于可能发生的情况进行训练。未能预测这些不太可能但后果严重的意外事件，这就是我们一开始就不应该过度依赖预测模型的原因，尤其是在没有讨论确定性或置信区间的情况下。

本章将探讨一个使用长短期记忆（LSTM）模型的多元预测问题。我们将首先使用传统的解释方法评估模型，然后使用我们在第七章“可视化卷积神经网络”中学习的集成梯度方法来生成我们模型的局部属性。

但更重要的是，我们将更好地理解 LSTM 的学习过程和局限性。然后，我们将采用预测近似方法以及 SHAP 的KernelExplainer进行全局和局部解释。最后，预测和不确定性是内在相关的，敏感性分析是一系列旨在衡量模型输出不确定性相对于其输入的方法，因此在预测场景中非常有用。我们还将研究两种这样的方法：Morris用于因素优先级排序和Sobol用于因素固定，这涉及到成本敏感性。

下面是我们将要讨论的主要主题：

• 使用传统解释方法评估时间序列模型

• 使用积分梯度生成 LSTM 属性

• 使用 SHAP 的KernelExplainer计算全局和局部属性

• 使用因素优先级识别有影响力的特征

• 使用因素固定量化不确定性和成本敏感性

让我们开始吧！

技术要求

本章的示例使用了mldatasets、pandas、numpy、sklearn、tensorflow、matplotlib、seaborn、alibi、distython、shap和SALib库。如何安装所有这些库的说明可以在本书的序言中找到。

本章的代码位于此处：packt.link/b6118。

任务

高速公路交通拥堵是一个影响世界各地的城市的问题。随着发展中国家每千人车辆数量的稳步增加，而道路和停车基础设施不足以跟上这一增长，拥堵水平已经达到了令人担忧的程度。在美国，每千人车辆统计数字是世界上最高的之一（2019 年为每千人 838 辆）。因此，美国城市占全球 381 个城市中至少有 15%拥堵水平的 62 个城市。

明尼阿波利斯就是这样一座城市（参见图 9.1），那里的阈值最近已经超过并持续上升。为了将这个大都市地区置于适当的背景中，拥堵水平在 50%以上时极为严重，但中等程度的拥堵（15-25%）已经是未来可能出现严重拥堵的预警信号。一旦拥堵达到 25%，就很难逆转，因为任何基础设施的改善都将非常昂贵，而且还会进一步扰乱交通。最严重的拥堵点之一是在明尼阿波利斯和圣保罗这对双城之间的 94 号州际公路（I-94），当通勤者试图缩短旅行时间时，这会拥堵替代路线。了解这一点后，这两座城市的市长们已经获得了一些联邦资金来扩建这条公路：

图 9.1：TomTom 的 2019 年明尼阿波利斯交通指数

市长们希望能够吹嘘一个完成的扩建项目作为共同成就，以便在第二任期内再次当选。然而，他们很清楚，一个嘈杂、脏乱和阻碍交通的扩建项目可能会给通勤者带来很大的麻烦，因此，如果扩建项目不是几乎看不见，它可能会在政治上适得其反。因此，他们规定建筑公司尽可能在其他地方预制，并在低流量时段进行组装。这些时段的每小时交通量少于 1,500 辆车。他们一次只能在一个方向的高速公路上工作，并且在他们工作时只能阻挡不超过一半的车道。为了确保遵守这些规定，如果他们在任何交通量超过这个阈值时阻挡超过四分之三的高速公路，他们将对公司处以每辆车 15 美元的罚款。

除了这些，如果施工队伍在每小时交通量超过 1,500 辆车时在现场阻挡一半的高速公路，他们每天将花费 5,000 美元。为了更直观地了解这一点，在典型的交通高峰时段阻挡可能会使建筑公司每小时损失 67,000 美元，再加上每天 5,000 美元的费用！当地当局将使用沿路线的自动交通记录器（ATR）站点来监控交通流量，以及当地交通警察来记录施工时车道被阻挡的情况。

该项目已被规划为一个为期 2 年的建设项目；第一年将在 I-94 路线的西行车道进行扩建，而第二年将扩建东行车道。施工现场的建设仅从 5 月到 10 月进行，因为在这几个月里下雪不太可能延误施工。在整个余下的年份，他们将专注于预制。他们将尝试只在工作日工作，因为工人联盟为周末谈判了慷慨的加班费。因此，只有在有重大延误的情况下，周末才会进行施工。然而，工会同意在 5 月至 10 月期间以相同的费率在节假日工作。

建筑公司不想承担任何风险！因此，他们需要一个模型来预测 I-94 路线的交通流量，更重要的是，要了解哪些因素会创造不确定性并可能增加成本。他们已经聘请了一位机器学习专家来完成这项工作：就是你！

建筑公司提供的 ATR 数据包括截至 2018 年 9 月的每小时交通量，以及同一时间尺度的天气数据。它只包括西行车道，因为那部分扩建将首先进行。

方法

您已经使用近四年的数据（2012 年 10 月 – 2016 年 9 月）训练了一个有状态的双向 LSTM模型。您保留了最后一年用于测试（2017 年 9 月–2018 年）和之前一年用于验证（2016 年 9 月 –2017 年）。这样做是有道理的，因为测试和验证数据集与高速公路扩建项目预期的条件（3 月 – 11 月）相吻合。您曾考虑使用仅利用这些条件数据的其他分割方案，但您不想如此大幅度地减少训练数据，也许它们最终还是可能需要用于冬季预测。回望窗口定义了时间序列模型可以访问多少过去数据。您选择了 168 小时（1 周）作为回望窗口大小。鉴于模型的这种有状态性质，随着模型在训练数据中的前进，它可以学习每日和每周的季节性，以及一些只能在几周内观察到的趋势和模式。您还训练了另外两个模型。您概述了以下步骤以满足客户期望：

1. 使用 RMSE、回归图、混淆矩阵 等等，您将访问模型的预测性能，更重要的是，了解误差的分布情况。

2. 使用 集成梯度，您将了解是否采取了最佳建模策略，因为它可以帮助您可视化模型到达决策的每一条路径，并帮助您根据这一点选择模型。

3. 使用 SHAP 的 KernelExplainer 和预测近似方法，您将推导出对所选模型有重要意义的特征的全局和局部理解。

4. 使用 Morris 敏感性分析，您将识别 因子优先级，它根据它们可以驱动输出变异性的程度对因素（换句话说，特征）进行排序。

5. 使用 Sobol 敏感性分析，您将计算 因子固定，这有助于确定哪些因素不具有影响力。它是通过量化输入因素对输出变异性的贡献和相互作用来做到这一点的。有了这个，您可以了解哪些因素可能对潜在的罚款和成本影响最大，从而产生基于变异性的成本敏感性分析。

准备工作

您可以在此处找到此示例的代码：github.com/PacktPublishing/Interpretable-Machine-Learning-with-Python-2E/blob/main/09/Traffic_compact1.ipynb。

加载库

要运行此示例，您需要安装以下库：

• mldatasets 用于加载数据集

• pandas 和 numpy 用于操作数据集

• tensorflow 用于加载模型

• scikit-learn、matplotlib、seaborn、alibi、distython、shap 和 SALib 用于创建和可视化解释

您应该首先加载所有这些内容：

import math
import os
import mldatasets
import pandas as pd
import numpy as np
from sklearn.preprocessing import MinMaxScaler
from sklearn import metrics
import tensorflow as tf
from tensorflow import keras
from tensorflow.keras.preprocessing.sequence import \ TimeseriesGenerator
from keras.utils import get_file
import matplotlib.pyplot as plt
from matplotlib.colors import TwoSlopeNorm
import seaborn as sns
from alibi.explainers import IntegratedGradients
from distython import HEOM
import shap
from SALib.sample import morris as ms
from SALib.analyze import morris as ma
from SALib.plotting import morris as mp
from SALib.sample.saltelli import sample as ss
from SALib.analyze.sobol import analyze as sa
from SALib.plotting.bar import plot as barplot 

让我们通过使用print(tf.__version__)命令来检查 TensorFlow 是否加载了正确的版本。它应该是 2.0 或更高版本。

理解和准备数据

are loading the data into a DataFrame called traffic_df. Please note that the prepare=True parameter is important because it performs necessary tasks such as subsetting the DataFrame to the required timeframe, since October 2015, some interpolation, correcting holidays, and performing one-hot encoding:

traffic_df = mldatasets.load("traffic-volume-v2", prepare=True) 

应该有超过 52,000 条记录和 16 列。我们可以使用traffic_df.info()来验证这一点。输出应该符合预期。所有特征都是数值型的，没有缺失值，并且分类特征已经为我们进行了独热编码。

数据字典

由于分类编码，只有九个特征，但它们变成了 16 列：

• dow: 序数型；以星期一开始的星期几（介于 0 到 6 之间）

• hr: 序数型；一天中的小时（介于 0 到 23 之间）

• temp: 连续型；摄氏度平均温度（介于-30 到 37 之间）

• rain_1h: 连续型；该小时发生的降雨量（介于 0 到 21 毫米之间）

• snow_1h: 连续型；该小时发生的雪量（当转换为液体形式时）（介于 0 到 2.5 厘米之间）

• cloud_coverage: 连续型；云层覆盖率百分比（介于 0 到 100 之间）

• is_holiday: 二元型；该天是星期一到星期五的国家或州假日吗？（1 表示是，0 表示否）？

• traffic_volume: 连续型；捕获交通量的目标特征

• weather: 分类；该小时天气的简短描述（晴朗 | 云层 | 雾 | 薄雾 | 雨 | 雪 | 未知 | 其他）

理解数据

理解时间序列问题的第一步是理解目标变量。这是因为它决定了你如何处理其他所有事情，从数据准备到建模。目标变量可能与时间有特殊的关系，例如季节性变动或趋势。

理解周

首先，我们可以从每个季节中采样一个 168 小时的时间段，以更好地理解一周中每天之间的方差，然后了解它们如何在季节和假日之间变化：

lb = 168
fig, (ax0,ax1,ax2,ax3) = plt.subplots(4,1, figsize=(15,8))
plt.subplots_adjust(top = 0.99, bottom=0.01, hspace=0.4)
traffic_df[(lb*160):(lb*161)].traffic_volume.plot(ax=ax0)
traffic_df[(lb*173):(lb*174)].traffic_volume.plot(ax=ax1)
traffic_df[(lb*186):(lb*187)].traffic_volume.plot(ax=ax2)
traffic_df[(lb*199):(lb*200)].traffic_volume.plot(ax=ax3) 

前面的代码生成了图 9.2中显示的图表。如果你从左到右阅读它们，你会看到它们都是从星期三开始，以下一周的星期二结束。每周的每一天都是从低点开始和结束，中间有一个高点。工作日通常有两个高峰，对应早晨和下午高峰时段，而周末只有一个下午的峰值：

图 9.2：代表每个季节的几个交通量样本周周期

存在一些主要异常值，例如 10 月 31 日星期六，这基本上是万圣节，并不是官方的节假日。还有 2 月 2 日（星期二）是严重的暴风雪的开始，而夏末的时期比其他样本周要混乱得多。结果发现，那一年州博览会发生了。像万圣节一样，它既不是联邦节日也不是地区节日，但重要的是要注意博览会场地位于明尼阿波利斯和圣保罗之间的一半。你还会注意到，在 7 月 29 日星期五午夜时，交通量有所上升，这可以归因于这是明尼阿波利斯音乐会的大日子。

在比较时间序列中的各个时期时，试图解释这些不一致性是一个很好的练习，因为它有助于你确定需要添加到模型中的变量，或者至少知道缺少了什么。在我们的案例中，我们知道我们的is_holiday变量不包括万圣节或整个州博览会周，也没有针对大型音乐或体育赛事的变量。为了构建一个更稳健的模型，寻找可靠的外部数据源并添加更多覆盖所有这些可能性的特征是明智的，更不用说验证现有变量了。目前，我们将使用我们拥有的数据。

理解日子

对于高速公路扩建项目来说，了解平均工作日的交通状况至关重要。施工队伍只在工作日（周一至周五）工作，除非遇到延误，在这种情况下，他们也会在周末工作。我们还必须区分节假日和其他工作日，因为这些可能有所不同。

为了这个目的，我们将创建一个 DataFrame（weekend_df）并创建一个新列（type_of_day），将小时编码为“假日”、“工作日”或“周末”。然后，我们可以按此列和hr列进行分组，并使用mean和标准差（std）进行聚合。然后我们可以进行pivot，以便我们有一个列，其中包含每个type_of_day类别的平均交通量和标准差，其中行代表一天中的小时数（hr）。然后，我们可以绘制结果 DataFrame。我们可以创建包含标准差的区间：

weekend_df = traffic_df[
    ['hr', 'dow', 'is_holiday', 'traffic_volume']].copy()
weekend_df['type_of_day'] = np.where(
    weekend_df.is_holiday == 1,
    'Holiday',
    np.where(weekend_df.dow >= 5, 'Weekend', 'Weekday')
)
weekend_df = weekend_df.groupby(
['type_of_day','hr']) ['traffic_volume']
    .agg(['mean','std'])
    .reset_index()
    .pivot(index='hr', columns='type_of_day', values=['mean', 'std']
)
weekend_df.columns = [
    ''.join(col).strip().replace('mean','')\
    for col in weekend_df.columns.values
]
fig, ax = plt.subplots(figsize=(15,8))
weekend_df[['Holiday','Weekday','Weekend']].plot(ax=ax)
plt.fill_between(
    weekend_df.index,
    np.maximum(weekend_df.Weekday - 2 * weekend_df.std_Weekday, 0),
    weekend_df.Weekday + 2 * weekend_df.std_Weekday,
    color='darkorange',
    alpha=0.2
)
plt.fill_between(
    weekend_df.index,\
    np.maximum(weekend_df.Weekend - 2 * weekend_df.std_Weekend, 0),
    weekend_df.Weekend + 2 * weekend_df.std_Weekend,
    color='green',
    alpha=0.1
)
plt.fill_between(
    weekend_df.index,\
    np.maximum(weekend_df.Holiday - 2 * weekend_df.std_Holiday, 0),
    weekend_df.Holiday + 2 * weekend_df.std_Holiday,
    color='cornflowerblue',
    alpha=0.1
) 

前面的代码片段产生了以下图表。它表示每小时平均交通量，但变化很大，这就是为什么建筑公司正在谨慎行事。图中绘制了代表每个阈值的水平线：

• 容量满载时为 5,300。

• 半容量时为 2,650，之后建筑公司将因每日指定金额被罚款。

• 无施工阈值是 1,500，之后建筑公司将因每小时指定金额被罚款。

他们只想在通常低于 1500 阈值的小时内工作，周一到周五。这五个小时将是晚上 11 点（前一天）到早上 5 点。如果他们必须周末工作，这个时间表通常会推迟到凌晨 1 点，并在早上 6 点结束。在工作日，变化相对较小，所以建筑公司坚持只在工作日工作是可以理解的。在这些小时里，节假日看起来与周末相似，但节假日的变化甚至比周末更大，这可能是更成问题的情况：

图 9.3：节假日、工作日和周末的平均每小时交通量，以及间隔

通常，对于这样的项目，你会探索预测变量，就像我们对目标所做的那样。这本书是关于模型解释的，所以我们将通过解释模型来了解预测变量。但在我们到达模型之前，我们必须为它们准备数据。

数据准备

第一步数据准备是将数据分割成训练集、验证集和测试集。请注意，测试数据集包括最后 52 周（2184小时），而验证数据集包括之前的 52 周，因此它从4368小时开始，到 DataFrame 最后一行之前的2184小时结束：

train = traffic_df[:-4368]
valid = traffic_df[-4368:-2184]
test = traffic_df[-2184:] 

现在 DataFrame 已经被分割，我们可以绘制它以确保其部分是按照预期分割的。我们可以使用以下代码来完成：

plt.plot(train.index.values, train.traffic_volume.values,
          label='train')
plt.plot(valid.index.values, valid.traffic_volume.values,
           label='validation')
plt.plot(test.index.values, test.traffic_volume.values,
          label='test')
plt.ylabel('Traffic Volume')
plt.legend() 

上述代码生成了图 9.4。它显示，训练数据集分配了近 4 年的数据，而验证和测试各分配了一年。在这个练习中，我们将不再引用验证数据集，因为它只是在训练期间作为工具来评估模型在每个 epoch 后的预测性能。

图 9.4：时间序列分割为训练集、验证集和测试集

下一步是对数据进行 min-max 归一化。我们这样做是因为较大的值会导致所有神经网络的学习速度变慢，而 LSTM 非常容易发生梯度爆炸和消失。相对均匀且较小的数字可以帮助解决这些问题。我们将在本章后面讨论这个问题，但基本上，网络要么在数值上不稳定，要么在达到全局最小值方面无效。

我们可以使用scikit包中的MinMaxScaler进行 min-max 归一化。目前，我们只会对归一化器进行fit操作，以便我们可以在需要时使用它们。我们将为我们的目标（traffic_volume）创建一个名为y_scaler的归一化器，并为其余变量（X_scaler）创建另一个归一化器，使用整个数据集，以确保无论使用哪个部分（train、valid或test），转换都是一致的。所有的fit过程只是保存公式，使每个变量适合在零和一之间：

y_scaler = MinMaxScaler()
y_scaler.fit(traffic_df[['traffic_volume']])
X_scaler = MinMaxScaler()
X_scaler.fit(traffic_df.drop(['traffic_volume'], axis=1)) 

现在，我们将使用我们的缩放器 transform 我们的训练和测试数据集，为每个创建 y 和 X 对：

y_train = y_scaler.transform(train[['traffic_volume']])
X_train = X_scaler.transform(train.drop(['traffic_volume'], axis=1))
y_test = y_scaler.transform(test[['traffic_volume']])
X_test = X_scaler.transform(test.drop(['traffic_volume'], axis=1)) 

然而，对于时间序列模型，我们创建的 y 和 X 对并不有用，因为每个观测值都是一个时间步长。每个时间步长不仅仅是该时间步长发生的特征，而且在一定程度上是它之前发生的事情，称为滞后。例如，如果我们根据 168 个滞后观测值预测交通，对于每个标签，我们将需要每个特征的之前 168 小时的数据。因此，你必须为每个时间步长以及其滞后生成一个数组。幸运的是，keras 有一个名为 TimeseriesGenerator 的函数，它接受你的 X 和 y 并生成一个生成器，该生成器将数据馈送到你的模型。你必须指定一个特定的 length，这是滞后观测值的数量（也称为 lookback window）。默认的 batch_size 是一个，但我们使用 24，因为客户更喜欢一次获取 24 小时的预测，而且使用更大的批次大小进行训练和推理要快得多。

自然地，当你需要预测明天时，你需要明天的天气，但你可以用天气预报来补充时间步长：

gen_train = TimeseriesGenerator(
    X_train,
    y_train,
    length=lb,
    batch_size=24
)
gen_test = TimeseriesGenerator(
    X_test,
    y_test,
    length=lb,
    batch_size=24
)
print(
    "gen_train:%s×%s→%s" % (len(gen_train),
    gen_train[0][0].shape, gen_train[0][1].shape)
)
print(
    "gen_test:%s×%s→%s" % (len(gen_test),
    gen_test[0][0].shape, gen_test[0][1].shape)
) 
gen_train) and the testing generator (gen_test), which use a length of 168 hours and a batch size of 24:

gen_train:  1454 ×   (24, 168, 15)   →   (24, 1)
gen_test:   357  ×   (24, 168, 15)   →   (24, 1) 

任何使用 1 周滞后窗口和 24 小时批次大小训练的模型都需要这个生成器。每个生成器是与每个批次对应的元组的列表。这个元组的索引 0 是 X 特征数组，而索引 1 是 y 标签数组。因此，输出的第一个数字是列表的长度，即批次的数量。X 和 y 数组的维度随后。

例如，gen_train 有 1,454 个批次，每个批次有 24 个时间步长，长度为 168，有 15 个特征。从这些 24 个时间步长中预期的预测标签的形状是 (24,1)。

最后，在继续处理模型和随机解释方法之前，让我们尝试通过初始化我们的随机种子来使事情更具可重复性：

rand = 9
os.environ['PYTHONHASHSEED']=str(rand)
tf.random.set_seed(rand)
np.random.seed(rand) 

加载 LSTM 模型

我们可以快速加载模型并像这样输出其摘要：

model_name = 'LSTM_traffic_168_compact1.hdf5'
model_path = get_file(
    model_name,
    'https://github.com/PacktPublishing/Interpretable-\ 
    Machine-Learning-with-Python-2E/blob/main/models/{}?raw=true'
    .format(model_name)
)
lstm_traffic_mdl = keras.models.load_model(model_path)
lstm_traffic_mdl.summary() 
bidirectional LSTM layer with an output of (24, 168). 24 corresponds to the batch size, while 168 means that there’s not one but two 84-unit LSTMs going in opposite directions and meeting in the middle. It has a dropout of 10%, and then a dense layer with a single ReLu-activated unit. The ReLu ensures that all the predictions are over zero since negative traffic volume makes no sense:

Model: "LSTM_traffic_168_compact1"
_________________________________________________________________
Layer (type)                 Output Shape              Param #   
=================================================================
Bidir_LSTM (Bidirectional)   (24, 168)                 67200     
_________________________________________________________________
Dropout (Dropout)            (24, 168)                 0         
_________________________________________________________________
Dense (Dense)                (24, 1)                   169       
=================================================================
Total params: 67,369
Trainable params: 67,369
Non-trainable params: 0
_________________________________________________________________ 

现在，让我们使用传统的解释方法评估 LSTM_traffic_168_compact1 模型。

使用传统的解释方法评估时间序列模型

时间序列回归模型可以像评估任何回归模型一样进行评估；也就是说，使用来自 均方误差 或 R-squared 分数的指标。当然，在某些情况下，你可能需要使用具有中位数、对数、偏差或绝对值的指标。这些模型不需要任何这些。

使用标准的回归指标

evaluate_reg_mdl 函数可以评估模型，输出一些标准的回归指标，并绘制它们。此模型的参数是拟合的模型 (lstm_traffic_mdl)，X_train (gen_train)，X_test (gen_test)，y_train 和 y_test。

可选地，我们可以指定一个y_scaler，以便模型使用标签的逆变换进行评估，这使得绘图和均方根误差（RMSE）更容易理解。在这种情况下，另一个非常必要的可选参数是y_truncate=True，因为我们的y_train和y_test的维度比预测标签大。这种差异发生是因为由于回望窗口，第一次预测发生在数据集的第一个时间步之后。因此，我们需要从y_train中减去这些时间步，以便与gen_train的长度匹配：

现在我们将使用以下代码评估这两个模型。为了观察预测的进度，我们将使用predopts={"verbose":1}。

y_train_pred, y_test_pred, y_train, y_test =\
    mldatasets.evaluate_reg_mdl(lstm_traffic_mdl,
        gen_train,
        gen_test,
        y_train,
        y_test,
        scaler=y_scaler,
        y_truncate=True,
        predopts={"verbose":1}
) 
Figure 9.5. The *regression plot* is, essentially, a scatter plot of the observed versus predicted traffic volumes, fitted to a linear regression model to show how well they match. These plots show that the model tends to predict zero traffic when it’s substantially higher. Besides that, there are a number of extreme outliers, but it fits relatively well with a test RMSE of 430 and only a slightly better train RMSE:

图 9.5：“LSTM_traffic_168_compact1”模型的预测性能评估

我们还可以通过比较观察到的与预测的交通来评估模型。按小时和类型分解错误可能会有所帮助。为此，我们可以创建包含这些值的 DataFrames - 每个模型一个。但首先，我们必须截断 DataFrame（-y_test_pred.shape[0]），以便它与预测数组的长度匹配，我们不需要所有列，所以我们只提供我们感兴趣的索引：traffic_volume是第 7 列，但我们还希望有dow（第 0 列）、hr（第 1 列）和is_holiday（第 6 列）。我们将traffic_volume重命名为actual_traffic，并创建一个名为predicted_traffic的新列，其中包含我们的预测。然后，我们将创建一个type_of_day列，就像我们之前做的那样，它告诉我们是否是节假日、工作日还是周末。最后，我们可以删除dow和is_holiday列，因为我们不再需要它们：

evaluate_df = test.iloc[-y_test_pred.shape[0]:,[0,1,6,7]]
    .rename(columns={'traffic_volume':'actual_traffic'}
)
evaluate_df['predicted_traffic'] = y_test_pred
evaluate_df['type_of_day'] = np.where(
    evaluate_df.is_holiday == 1,
    'Holiday',
    np.where(evaluate_df.dow >= 5,
    'Weekend', 'Weekday')
)
evaluate_df.drop(['dow','is_holiday'], axis=1, inplace=True) 

您可以通过简单地运行一个带有evaluate_df的单元格来快速查看 DataFrames 的内容。它应该有 4 列。

预测误差聚合

可能是某些日期和时间更容易出现预测误差。为了更好地了解这些误差如何在时间上分布，我们可以按小时分段绘制type_of_day的 RMSE。为此，我们必须首先定义一个rmse函数，然后按type_of_day和hr对每个模型的评估 DataFrame 进行分组，并使用apply函数通过rmse函数进行聚合。然后我们可以通过转置来确保每个type_of_day都有一个按小时显示 RMSE 的列。然后我们可以平均这些列并将它们存储在一个系列中：

def rmse(g):
    rmse = np.sqrt(
    metrics.mean_squared_error(g['actual_traffic'],
                               g['predicted_traffic'])
    )
    return pd.Series({'rmse': rmse})
evaluate_by_hr_df = evaluate_df.groupby(['type_of_day', 'hr'])
    .apply(rmse).reset_index()
    .pivot(index='hr', columns='type_of_day', values='rmse')
mean_by_daytype_s = evaluate_by_hr_df.mean(axis=0) 

现在我们有了包含节假日、工作日和周末每小时 RMSE 的 DataFrames，以及这些“类型”的日平均数，我们可以使用evaluate_by_hr DataFrame 来绘制它们。我们还将创建带有每个type_of_day平均值的虚线水平线，这些平均值来自mean_by_daytype pandas系列：

evaluate_by_hr_df.plot()
ax = plt.gca()
ax.set_title('Hourly RMSE distribution', fontsize=16)
ax.set_ylim([0,2500])
ax.axhline(
    y=mean_by_daytype_s.Holiday,
    linewidth=2,
    color='cornflowerblue',
    dashes=(2,2)
)
ax.axhline(
    y=mean_by_daytype_s.Weekday,
    linewidth=2,
    color='darkorange',
    dashes=(2,2)
)
ax.axhline(
    y=mean_by_daytype_s.Weekend,
    linewidth=2,
    color='green',
    dashes=(2,2)
) 

上述代码生成了图 9.6中显示的图表。正如我们所见，该模型在假日有很高的 RMSE。然而，该模型可能高估了交通量，而在这种特定情况下，高估不如低估糟糕，因为低估可能导致交通延误和额外的罚款成本：

图 9.6：“LSTM_traffic_168_compact1”模型按 type_of_day 类型划分的小时 RMSE

将模型评估视为分类问题

的确，就像分类问题可以有假阳性和假阴性，其中一个是比另一个更昂贵的，你可以用诸如低估和过度估计等概念来构建任何回归问题。这种构建特别有用，当其中一个比另一个更昂贵时。如果你有明确定义的阈值，就像我们在这个项目中做的那样，你可以像评估分类问题一样评估任何回归问题。我们将使用半容量和“无施工”阈值混淆矩阵来评估它。为了完成这项任务，我们可以使用np.where来获取实际值和预测值超过每个阈值的二进制数组。然后我们可以使用compare_confusion_matrices函数来比较模型的混淆矩阵：

actual_over_half_cap = np.where(evaluate_df['actual_traffic'] >\
                                2650, 1, 0)
pred_over_half_cap = np.where(evaluate_df['predicted_traffic'] >\
                              2650, 1, 0)
actual_over_nc_thresh = np.where(evaluate_df['actual_traffic'] >\
                                 1500, 1, 0)
pred_over_nc_thresh = np.where(evaluate_df['predicted_traffic'] >\
                               1500, 1, 0)
mldatasets.compare_confusion_matrices(
    actual_over_half_cap,
    pred_over_half_cap,
    actual_over_nc_thresh,
    pred_over_nc_thresh,
    'Over Half-Capacity',
    'Over No-Construction Threshold'
) 
Figure 9.7.

图 9.7：“LSTM_traffic_168_compact1”模型的超过半容量和“无施工”阈值的混淆矩阵

我们最感兴趣的是假阴性（左下象限）的百分比，因为当实际上交通量超过了阈值时，预测没有超过阈值将导致高额罚款。另一方面，假阳性的成本在于在交通量实际上没有超过阈值的情况下提前离开施工现场。尽管如此，安全总是比后悔好！如果你比较“无施工”阈值的假阴性（0.85%），它不到半容量阈值（3.08%）的三分之一。最终，最重要的是无施工阈值，因为目的是在接近半容量之前停止施工。

现在我们已经利用传统方法来理解模型的决策，让我们继续探讨一些更高级的模型无关方法。

使用集成梯度生成 LSTM 归因

我们第一次在第七章，可视化卷积神经网络中了解到集成梯度（IG）。与该章节中研究的其他基于梯度的归因方法不同，路径集成梯度不依赖于卷积层，也不限于分类问题。

事实上，因为它计算了输出相对于输入沿路径平均的梯度，所以输入和输出可以是任何东西！通常与卷积神经网络（CNNs）和循环神经网络（RNNs）一起使用整合梯度，就像我们在本章中解释的那样。坦白说，当你在网上看到 IG LSTM 的例子时，它有一个嵌入层，是一个 NLP 分类器，但 IG 对于甚至处理声音或遗传数据的 LSTMs 也非常有效！

整合梯度解释器和我们将继续使用的解释器可以访问交通数据集的任何部分。首先，让我们为所有这些创建一个生成器：

y_all = y_scaler.transform(traffic_df[['traffic_volume']])
X_all = X_scaler.transform(
    traffic_df.drop(['traffic_volume'], axis=1)
)
gen_all = TimeseriesGenerator(
    X_all, y_all, length=lb, batch_size=24
) 

整合梯度是一种局部解释方法。所以，让我们获取一些我们可以解释的“感兴趣的样本实例”。我们知道假日可能需要专门的逻辑，所以让我们看看我们的模型是否注意到了is_holiday在某个例子（holiday_afternoon_s）中的重要性。早晨也是一个问题，尤其是由于天气条件，早晨的拥堵时间比平均水平更长，所以我们有一个例子（peak_morning_s）。最后，一个炎热的日子可能会有更多的交通，尤其是在周末（hot_Saturday_s）：

X_df = traffic_df.drop(['traffic_volume'], axis=1 \
    ).reset_index(drop=True)
holiday_afternoon_s = X_df[
    (X_df.index >= 43800) & (X_df.dow==0) &\
    (X_df.hr==16) &(X_df.is_holiday==1)
].tail(1)
peak_morning_s = X_df[
    (X_df.index >= 43800) & (X_df.dow==2) &\
    (X_df.hr==8) & (X_df.weather_Clouds==1) & (X_df.temp<20)
].tail(1)
hot_Saturday_s = X_df[
    (X_df.index >= 43800) & (X_df.dow==5) &\
    (X_df.hr==12) & (X_df.temp>29) & (X_df.weather_Clear==1)
].tail(1) 

现在我们已经创建了一些实例，让我们实例化我们的解释器。来自alibi包的IntegratedGradients只需要一个深度学习模型，但建议为积分近似设置步骤数（n_steps）和内部批次大小。我们将为我们的模型实例化一个解释器：

ig = IntegratedGradients(
    lstm_traffic_mdl, n_steps=25, internal_batch_size=24
) 

在我们对样本和解释器进行迭代之前，重要的是要意识到我们需要如何将样本输入到解释器中，因为它需要一个包含 24 个样本的批次。为此，一旦我们减去回望窗口（nidx），我们就必须获取样本的索引。然后，你可以从生成器（gen_all）中获取该样本的批次。每个批次包含 24 个时间步长，所以你需要将nidx向下取整到 24（nidx//24），以获取该样本的批次位置。一旦你获取了该样本的批次（batch_X）并打印了形状（24, 168, 15），第一个数字是 24 这一点不应该让你感到惊讶。当然，我们还需要获取批次内样本的索引（nidx%24），以获取该样本的数据：

nidx = holiday_afternoon_s.index.tolist()[0] – lb
batch_X = gen_all[nidx//24][0]
print(batch_X.shape) 

for循环将使用之前解释的方法来定位样本的批次（batch_X）。这个batch_X被输入到explain函数中。这是因为这是一个回归问题，没有目标类别；也就是说，target=None。一旦产生了解释，attributions属性将包含整个批次的属性。我们只能获取样本的属性，并将其transpose以产生一个形状为(15, lb)的图像。for循环中的其余代码只是获取用于刻度的标签，然后绘制一个图像，该图像拉伸以适应我们的figure维度，以及其标签：

samples = [holiday_afternoon_s, peak_morning_s, hot_saturday_s]
sample_names = ['Holiday Afternoon', 'Peak Morning' , 'Hot Saturday']
for s in range(len(samples)):
    nidx = samples[s].index.tolist()[0] - lb
    batch_X = gen_all[nidx//24][0]
    explanation = ig.explain(batch_X, target=None)
    attributions = explanation.attributions[0]
    attribution_img = np.transpose(attributions[nidx%24,:,:])
    end_date = traffic_df.iloc[samples[s].index
        ].index.to_pydatetime()[0]
    date_range = pd.date_range(
        end=end_date, periods=8, freq='1D').to_pydatetime().tolist()
    columns = samples[s].columns.tolist()  
    plt.title(
        'Integrated Gradient Attribution Map for "{}"'.\
        format(sample_names[s], lb), fontsize=16
    )
    divnorm = TwoSlopeNorm(
        vmin=attribution_img.min(),
        vcenter=0,
        vmax=attribution_img.max()
    )
    plt.imshow(
        attribution_img,
        interpolation='nearest' ,
        aspect='auto',
        cmap='coolwarm_r',
        norm=divnorm
)
    plt.xticks(np.linspace(0,lb,8).astype(int), labels=date_range)
    plt.yticks([*range(15)], labels=columns)
    plt.colorbar(pad=0.01,fraction=0.02,anchor=(1.0,0.0))
    plt.show() 

上述代码将生成图 9.8中显示的图表。在y轴上，您可以看到变量名称，而在x轴上，您可以看到对应于所讨论样本回望窗口的日期。x轴的最右侧是样本的日期，随着您向左移动，您会向时间后退。例如，假日下午的样本是 9 月 3 日下午 4 点，有一周的回望，所以每个向后的刻度代表该日期前一天。

图 9.8：对于“LSTM_traffic_168_compact1”模型所有样本的注释集成梯度归因图

您可以通过查看图 9.8中的归因图强度来判断哪些小时/变量对预测很重要。每个归因图右侧的颜色条可以用作参考。红色中的负数表示负相关，而蓝色中的正数表示正相关。然而，一个相当明显的是，随着每个图向时间后退，强度往往会减弱。由于它是双向的，所以这种情况发生在两端。令人惊讶的是，这个过程发生得有多快。

让我们从底部开始。对于“热周六”，随着您接近预测时间（周六中午），星期几、小时、温度和晴朗的天气在这个预测中扮演的角色越来越重要。天气开始较凉爽，这解释了为什么在温度特征中红色区域出现在蓝色区域之前。

对于“高峰早晨”，归因是有意义的，因为它在之前有雨和多云之后变得晴朗，这导致高峰时段迅速达到顶峰而不是缓慢增加。在一定程度上，LSTM 已经学会了只有最近的天气才重要——不超过两三天。然而，集成梯度减弱的原因不仅仅是这一点。它们也因为梯度消失问题而减弱。这个问题发生在反向传播过程中，因为梯度值在每一步都要乘以权重矩阵，所以梯度可以指数级减少到零。

LSTM 被组织在一个非常长的序列中，这使得网络在长期捕捉依赖关系方面越来越无效。幸运的是，这些 LSTM 是有状态的，这意味着它们通过利用前一个批次的状态将批次按顺序连接起来。状态性确保了从长序列中学习，尽管存在梯度消失问题。这就是为什么当我们观察“假日下午”的归因图时，对于is_holiday有负归因，这是预料到没有高峰时段的合理原因。结果证明，9 月 3 日（劳动节）距离前一个假日（独立日）近两个月，而独立日是一个更盛大的节日。模型是否可能捕捉到这些模式呢？

我们可以尝试根据交通模式对假日进行子分类，看看这是否能帮助模型识别它们。我们还可以对以前的天气条件进行滚动汇总，以便模型更容易地捕捉到最近的天气模式。天气模式跨越数小时，因此汇总是直观的，而且更容易解释。解释方法可以为我们指明如何改进模型的方向，当然，改进的空间很大。

接下来，我们将尝试一种基于排列的方法！

使用 SHAP 的 KernelExplainer 计算全局和局部归因

排列方法通过对输入进行修改来评估它们将对模型输出产生多大的影响。我们首次在第四章，全局模型无关解释方法中讨论了这一点，但如果你还记得，有一个联盟框架可以执行这些排列，从而为不同特征的联盟产生每个特征的边际贡献的平均值。这个过程的结果是Shapley 值，它具有如加法和对称性等基本数学性质。不幸的是，对于不是特别小的数据集，Shapley 值的计算成本很高，所以 SHAP 库有近似方法。其中一种方法就是KernelExplainer，我们在第四章中也解释了它，并在第五章，局部模型无关解释方法中使用它。它使用加权局部线性回归来近似 Shapley 值，就像 LIME 所做的那样。

为什么使用 KernelExplainer？

我们有一个深度学习模型，那么为什么我们不使用 SHAP 的DeepExplainer，就像我们在第七章，可视化卷积神经网络中使用的 CNN 一样呢？DeepExplainer 将 DeepLIFT 算法改编来近似 Shapley 值。它与任何用于表格数据、CNN 和具有嵌入层的 RNN（例如用于 NLP 分类器或用于检测基因组序列的 RNN）都配合得非常好。对于多元时间序列，它变得更加复杂，因为 DeepExplainer 不知道如何处理输入的三维数组。即使它知道，它还包括了之前时间步的数据，因此你无法在不考虑之前时间步的情况下对单个时间步进行排列。例如，如果排列规定温度降低五度，这不应该影响之前数小时内的所有温度吗？如果温度降低 20 度呢？这不意味着它可能处于不同的季节，并且天气完全不同——也许还有更多的云和雪？

SHAP 的KernelExplainer可以接收任何任意的黑盒predict函数。它还对输入维度做出了一些假设。幸运的是，我们可以在它排列之前更改输入数据，使其对KernelExplainer来说，就像它正在处理一个表格数据集一样。任意的predict函数不必简单地调用模型的predict函数——它可以在输入和输出过程中更改数据！

定义一个策略使其与多元时间序列模型一起工作

为了模仿基于排列输入数据的可能过去天气模式，我们可以创建一个生成模型或类似的东西。这种策略将帮助我们生成适合排列时间步的多种过去时间步，以及为特定类别生成图像。尽管这可能会导致更准确的预测，但我们不会使用这种策略，因为它非常耗时。

相反，我们将使用gen_all生成器中的现有示例来找到最适合排列输入的时间序列数据。我们可以使用距离度量来找到最接近排列输入的那个。然而，我们必须设置一些限制，因为如果排列是在周六早上 5 点，温度为 27 摄氏度，云量为 90%，那么最接近的观察可能是在周五早上 7 点，但无论天气交通如何，它都会完全不同。因此，我们可以实现一个过滤器函数，确保它只找到相同dow、is_holiday和hr的最近观察。过滤器函数还可以清理排列样本，删除或修改模型中任何无意义的部分，例如分类特征的连续值：

图 9.9：排列近似策略

图 9.9展示了使用距离函数找到修改后的排列样本最近观察的过程。此函数返回最近的观察索引，但模型不能对单个观察（或时间步）进行预测，因此它需要其过去直到lookback窗口的小时历史。因此，它从生成器中检索正确的批次并对其进行预测，但预测将处于不同的尺度上，因此它们需要使用y_scaler进行逆变换。一旦predict函数迭代了所有样本并对它们进行了预测和缩放，它将它们发送回KernelExplainer，该工具输出它们的 SHAP 值。

为排列近似策略打下基础

您可以定义一个自定义的过滤器函数（filt_fn）。它接受一个包含整个数据集（X_df）的pandas DataFrame，您希望从中过滤，以及用于过滤的排列样本（x）和lookback窗口的长度。

该函数还可以修改排列后的样本。在这种情况下，我们必须这样做，因为模型中有许多特征是离散的，但排列过程使它们变得连续。正如我们之前提到的，所有过滤操作所做的只是通过限制选项来保护距离函数，防止它找到排列样本的非合理最近样本：

def filt_fn(X_df, x, lookback):
    x_ = x.copy()
    x_[0] = round(x_[0]) #round dow
    x_[1] = round(x_[1]) #round hr
    x_[6] = round(x_[6]) #round is_holiday
    if x_[1] < 0:#if hr < 0
        x_[1] = 24 + x_[1]
        x_[0] = x_[0] – 1  #make it previous day
    if x_[0] < 0:#if dow < 0
        x_[0] = 7 + x_[0] #make it previous week
        X_filt_df = X_df[
            (X_df.index >= lookback) & (X_df.dow==x_[0]) &\
            (X_df.hr==x_[1]) & (X_df.is_holiday==x_[6]) &\
            (X_df.temp-5<=x_[2]) & (X_df.temp+5>=x_[2])
        ]
    return X_filt_df, x_ 

如果你参考 图 9.9，在过滤器函数之后，我们接下来应该定义的是距离函数。我们可以使用 scipy.spatial.distance.cdist 接受的任何标准距离函数，例如“欧几里得”、“余弦”或“汉明”。这些标准距离函数的问题在于，它们要么与连续变量很好地工作，要么与离散变量很好地工作，但不能两者都很好地工作。我们在这个数据集中两者都有！

幸运的是，存在一些可以处理这两种情况的替代方案，例如异构欧几里得-重叠度量（HEOM）和异构值差异度量（HVDM）。这两种方法根据变量的性质应用不同的距离度量。HEOM 使用归一化的欧几里得距离 对连续变量，对离散变量使用“重叠”距离；即如果相同则为零距离，否则为 1。

HVDM 更复杂，因为对于连续变量，它是两个值之间的绝对距离，除以所涉及特征的均方根的四倍 )，这是一个处理异常值的好距离度量。对于离散变量，它使用归一化的值差异度量，这是基于两个值的条件概率之间的差异。

尽管 HVDM 对于具有许多连续值的集合数据集比 HEOM 更好，但在这种情况下却是过度设计。一旦数据集通过星期几 (dow) 和小时 (hr) 过滤，剩余的离散特征都是二进制的，因此“重叠”距离是理想的，而对于剩下的三个连续特征（temp、rain_1h、snow_1h 和 cloud_coverage），欧几里得距离应该足够。distython 有一个 HEOM 距离方法，它只需要一个背景数据集 (X_df.values) 和分类特征的索引 (cat_idxs)。我们可以使用 np.where 命令编程识别这些特征。

如果你想验证这些是否是正确的，请在单元格中运行 print(cat_idxs)。只有索引 2、3、4 和 5 应该被省略：

cat_idxs = np.where(traffic_df.drop(['traffic_volume'],\
                                    axis=1).dtypes != np.float64)[0]
heom_dist = HEOM(X_df.values, cat_idxs)
print(cat_idxs) 

现在，我们可以创建一个 lambda 函数，将 图 9.9 中描述的所有内容放在一起。它利用一个名为 approx_predict_ts 的函数来处理整个流程。它接受我们的过滤器函数 (filt_fn)、距离函数 (heom_dist.heom)、生成器 (gen_all) 和拟合的模型 (lstm_traffic_mdl)，并将它们链接在一起，如 图 9.9 所示。它还使用我们的缩放器 (X_scaler 和 y_scaler) 对数据进行缩放。距离是在转换后的特征上计算的，以提高准确性，并且预测在输出过程中进行反向转换：

predict_fn = lambda X: mldatasets.approx_predict_ts(
    X, X_df,
    gen_all,
    lstm_traffic_mdl,
    dist_metric=heom_dist.heom,
    lookback=lookback,
    filt_fn=filt_fn,
    X_scaler=X_scaler,
    y_scaler=y_scaler
) 

我们现在可以使用KernelExplainer的预测函数，但应该在最能代表施工队预期工作条件的样本上进行；也就是说，他们计划在 3 月到 11 月工作，最好是工作日和交通量低的时间。为此，让我们创建一个只包括这些月份的 DataFrame（working_season_df），并使用predict_fn和 DataFrame 的 k-means 作为背景数据初始化一个KernelExplainer：

working_season_df =\
    traffic_df[lookback:].drop(['traffic_volume'], axis=1).copy()
working_season_df =\
    working_season_df[(working_season_df.index.month >= 3) &\
                      (working_season_df.index.month <= 11)]
explainer = shap.KernelExplainer(
    predict_fn, shap.kmeans(working_season_df.values, 24)
) 

我们现在可以为working_season_df DataFrame 的随机观测值集生成 SHAP 值。

计算 SHAP 值

我们将从其中采样 48 个观测值。KernelExplainer相当慢，尤其是在使用我们的近似方法时。为了获得最佳的全球解释，最好使用大量的观测值，同时也要使用高nsamples，这是在解释每个预测时需要重新评估模型次数的数量。不幸的是，如果每种都有 50 个，那么解释器运行起来将需要数小时，这取决于你的可用计算资源，所以我们将会使用nsamples=10。你可以查看 SHAP 的进度条并相应地调整。一旦完成，它将生成包含 SHAP 值的特征重要性summary_plot：

X_samp_df = working_season_df.sample(80, random_state=rand)
shap_values = explainer.shap_values(X_samp_df, nsamples=10)
shap.summary_plot(shap_values, X_samp_df) 

上述代码绘制了以下图形中显示的摘要。不出所料，hr和dow是最重要的特征，其次是某些天气特征。奇怪的是，温度和降雨似乎并没有在预测中起到作用，但晚春到秋季可能不是一个显著因素。或者，也许更多的观测值和更高的nsample将产生更好的全球解释：

图 9.10：基于 48 个采样观测值产生的 SHAP 摘要图

我们可以用上一节中选择的感兴趣实例进行相同的操作，以进行局部解释。让我们遍历所有这些数据点。然后，我们可以生成一个单一的shap_values，但这次使用nsamples=80，然后为每个生成一个force_plot：

for s in range(len(samples)):
    print('Local Force Plot for "{}"'.format(sample_names[s]))
    shap_values_single = explainer.shap_values(
        datapoints[i], nsamples=80)
    shap.force_plot(
    explainer.expected_value,
    shap_values_single[0],
    samples[s],
    matplotlib=True
)
    plt.show() 

上述代码生成了图 9.11中显示的图形。“假日午后”的小时数(hr=16)推动预测值升高，而它是星期一(dow=0)和假日(is_holiday=1)的事实则推动预测值向相反方向移动。另一方面，“高峰早晨”主要由于小时数(hr=8.0)而处于高峰状态，但它有高cloud_coverage，肯定的weather_Clouds，而且没有降雨(rain_1h=0.0)。最后，“炎热的周六”由于星期数(dow=5)推动值降低，但异常高的值主要由于它是中午没有降雨和云层。奇怪的是，高于正常温度不是影响因素之一：

图 9.11：使用 SHAP 值和 nsamples=80 生成的力图，用于假日午后、高峰早晨和炎热周六

使用 SHAP 基于博弈论的方法，我们可以衡量现有观察值的排列如何使预测结果在许多可能特征联盟中边际变化。然而，这种方法可能非常有限，因为我们的背景数据中现有的方差塑造了我们对于结果方差的理解。

在现实世界中，变异性通常由数据中未表示的内容决定——但可能性极小。例如，在明尼阿波利斯夏季凌晨 5 点之前达到 25°C（77°F）并不常见，但随着全球变暖，它可能会变得频繁，因此我们想要模拟它如何影响交通模式。预测模型特别容易受到风险的影响，因此模拟是评估这种不确定性的关键解释组成部分。对不确定性的更好理解可以产生更稳健的模型，并直接指导决策。接下来，我们将讨论我们如何使用敏感性分析方法产生模拟。

使用因素优先级识别有影响力的特征

莫里斯方法是几种全局敏感性分析方法之一，范围从简单的分数因子到复杂的蒙特卡洛过滤。莫里斯位于这个光谱的某个位置，分为两个类别。它使用一次一个采样，这意味着在连续模拟之间只有一个值发生变化。它也是一个基本效应（EE）方法，这意味着它不量化模型中因素的确切效应，而是衡量其重要性和与其他因素的关系。顺便说一句，因素只是另一个在应用统计学中常用的特征或变量的名称。为了与相关理论保持一致，我们将在本节和下一节中使用这个词汇。

莫里斯的另一个特性是，它比我们接下来将要研究的基于方差的计算方法更节省计算资源。它可以提供比回归、导数或基于因子的简单且成本较低的方法更多的见解。它不能精确量化效应，但可以识别那些具有可忽略或交互效应的效应，这使得它成为在因素数量较少时筛选因素的理想方法。筛选也被称为因素优先级，因为它可以根据它们的分类来优先考虑因素。

计算莫里斯敏感性指数

莫里斯方法推导出与单个因素相关联的基本效应分布。每个基本效应分布都有一个平均值（µ）和标准差（σ）。这两个统计数据有助于将因素映射到不同的分类中。当模型非单调时，平均值可能是负数，因此莫里斯方法的一个变体通过绝对值（µ^*）进行调整，以便更容易解释。我们在这里将使用这种变体。

现在，让我们将此问题的范围限制在更易于管理的范围内。施工队将面临的道路交通不确定性将持续从 5 月到 10 月，周一至周五，晚上 11 点到凌晨 5 点。因此，我们可以从working_season_df DataFrame 中进一步提取子集，以生成一个工作小时 DataFrame（working_hrs_df），我们可以对其进行describe。我们将包括 1%、50%和 99%的百分位数，以了解中位数和异常值所在的位置：

working_hrs_df = working_season_df[
    (working_season_df.dow < 5)
    & ((working_season_df.hr < 5) | (working_season_df.hr > 22))
]
working_hrs_df.describe(percentiles=[.01,.5,.99]).transpose() 

上述代码生成了图 9.12中的表格。我们可以使用这张表格来提取我们在模拟中使用的特征范围。通常，我们会使用超过现有最大值或最小值的合理值。对于大多数模型，任何特征值都可以在其已知限制之外增加或减少，并且由于模型学习到了单调关系，它可以推断出合理的结局。例如，它可能学习到超过某个点的降雨量将逐渐减少交通。那么，假设你想模拟每小时 30 毫米的严重洪水；它可以准确预测无交通：

图 9.12：施工队计划工作期间的汇总统计

然而，因为我们使用的是从历史值中采样的预测近似方法，所以我们受到如何将边界推到已知范围之外的限制。因此，我们将使用 1%和 99%的百分位数值作为我们的限制。我们应该注意，这对于任何发现来说都是一个重要的注意事项，特别是对于可能超出这些限制的特征，例如temp、rain_1h和snow_1h。

从图 9.12的总结中，我们还需要注意的一点是，许多与天气相关的二元特征非常稀疏。你可以通过它们的极低平均值来判断。每个添加到敏感性分析模拟中的因素都会减慢其速度，因此我们只会选择前三个；即weather_Clear、weather_Clouds和weather_Rain。这些因素与其他六个因素一起在“问题”字典（morris_problem）中指定，其中包含它们的对应names、bounds和groups。现在，bounds是关键，因为它表示每个因素将模拟哪些值范围。我们将使用[0,4]（周一至周五）作为dow的值，以及[-1,4]（晚上 11 点到凌晨 4 点）作为hr的值。过滤器函数自动将负小时转换为前一天的小时，因此周二的一 1 相当于周一的 23。其余的界限是由百分位数确定的。请注意，groups中的所有因素都属于同一组，除了三个天气因素：

morris_problem = {
    # There are nine variables
    'num_vars': 10,
    # These are their names
    'names': ['dow', 'hr', 'temp', 'rain_1h', 'snow_1h',\
              'cloud_coverage', 'is_holiday', 'weather_Clear',\
              'weather_Clouds', 'weather_Rain'],
    # Plausible ranges over which we'll move the variables
    'bounds': [
        [0, 4], # dow Monday - Firday
        [-1, 4], # hr
        [-12, 25.], # temp (C)
        [0., 3.1], # rain_1h
        [0., .3], # snow_1h
        [0., 100.], # cloud_coverage
        [0, 1], # is_holiday
        [0, 1], # weather_Clear
        [0, 1], # weather_Clouds
        [0, 1] # weather_Rain
    ],
    # Only weather is grouped together
    'groups': ['dow', 'hr', 'temp', 'rain_1h', 'snow_1h',\
                'cloud_coverage', 'is_holiday', 'weather', 'weather',\
                'weather']
} 

一旦定义了字典，我们就可以使用 SALib 的 sample 方法生成 Morris 方法样本。除了字典外，它还需要轨迹数量（256）和级别（num_levels=4）。该方法使用因素和级别的网格来构建输入随机逐个移动的轨迹（OAT）。这里需要注意的重要一点是，更多的级别会增加这个网格的分辨率，可能使分析更好。然而，这可能会非常耗时。最好从轨迹数量和级别之间的比例 25:1 或更高开始。

然后，你可以逐步降低这个比例。换句话说，如果你有足够的计算能力，你可以让 num_levels 与轨迹数量相匹配，但如果你有这么多可用的计算能力，你可以尝试 optimal_trajectories=True。然而，鉴于我们有组，local_optimization 必须设置为 False。sample 的输出是一个数组，每个因素一列，(G + 1) × T 行（其中 G 是组数，T 是轨迹数）。我们有八个组，256 个轨迹，所以 print 应该输出一个 2,304 行 10 列的形状：

morris_sample = ms.sample(morris_problem, 256,\
                          num_levels=4, seed=rand)
print(morris_sample.shape) 

由于 predict 函数只与 15 个因素一起工作，我们应该修改样本，用零填充剩余的五个因素。我们使用零，因为这这些特征的中位数。中位数最不可能增加交通量，但你应该根据具体情况调整默认值。如果你还记得我们 第二章 的 心血管疾病 （CVD）示例，可解释性关键概念，增加 CVD 风险的特性值有时是最小值或最大值。

np.hstack 函数可以将数组水平拼接，使得前八个因素之后跟着三个零因素。然后，有一个孤独的第九个样本因素对应于 weather_Rain，接着是两个零因素。结果数组应该与之前一样行数，但列数为 15：

morris_sample_mod = np.hstack(
    (
        morris_sample[:,0:9],
        np.zeros((morris_sample.shape[0],3)),
        morris_sample[:,9:10],
        np.zeros((morris_sample.shape[0],2))
    )
)
print(morris_sample_mod.shape) 

被称为 morris_sample_mod 的 numpy 数组现在以我们的 predict 函数可以理解的形式包含了 Morris 样本。如果这是一个在表格数据集上训练过的模型，我们就可以直接利用模型的 predict 函数。然而，就像我们使用 SHAP 一样，我们必须使用近似方法。这次，我们不会使用 predict_fn，因为我们想在 approx_predict_ts 中设置一个额外的选项，progress_bar=True。其他一切都将保持不变。进度条将很有用，因为这可能需要一段时间。运行单元格，休息一下喝杯咖啡：

morris_preds = mldatasets.approx_predict_ts(
    morris_sample_mod,
    X_df,
    gen_all,
    lstm_traffic_mdl,
    filt_fn=filt_fn,
    dist_metric=heom_dist.heom,
    lookback=lookback,
    X_scaler=X_scaler,
    y_scaler=y_scaler,
    progress_bar=True
) 

要使用SALib的analyze函数进行敏感性分析，你需要你的问题字典（morris_problem），原始的 Morris 样本（morris_sample），以及我们用这些样本生成的预测（morris_preds）。还有一个可选的置信区间水平参数（conf_level），但默认的 0.95 是好的。它使用重采样来计算这个置信水平，默认为 1,000。这个设置也可以通过可选的num_resamples参数来改变：

morris_sensitivities = ma.analyze(
    morris_problem, morris_sample, morris_preds,\
    print_to_console=False
) 

分析基本影响

analyze将返回一个包含 Morris 敏感性指数的字典，包括平均值（µ）和标准差（σ）的基本影响，以及平均值（µ^）的绝对值。在表格格式中更容易欣赏这些值，这样我们就可以将它们放入 DataFrame 中，并根据µ*^*排序和着色，µ^*可以解释为因素的整体重要性。另一方面，σ表示因素与其它因素的交互程度：

morris_df = pd.DataFrame(
    {
        'features':morris_sensitivities['names'],
        'μ':morris_sensitivities['mu'],
        'μ*':morris_sensitivities['mu_star'],
        'σ':morris_sensitivities['sigma']
    }
)
morris_df.sort_values('μ*', ascending=False).style\
    .background_gradient(cmap='plasma', subset=['μ*']) 

前面的代码输出了图 9.13中展示的 DataFrame。你可以看出is_holiday是其中最重要的因素之一，至少在问题定义中指定的范围（morris_problem）内是这样。还有一点需要注意，天气确实有绝对的基本影响，但交互效应并不确定。组别很难评估，尤其是当它们是稀疏的二进制因素时：

图 9.13：因素的基本影响分解

前面的图中的 DataFrame 不是可视化基本影响的最佳方式。当因素不多时，更容易绘制它们。SALib提供了两种绘图方法。水平条形图（horizontal_bar_plot）和协方差图（covariance_plot）可以并排放置。协方差图非常好，但它没有注释它所界定的区域。我们将在下一节中了解这些。因此，仅出于教学目的，我们将使用text来放置注释：

fig, (ax0, ax1) = plt.subplots(1,2, figsize=(12,8))
mp.horizontal_bar_plot(ax0, morris_sensitivities, {})
mp.covariance_plot(ax1, morris_sensitivities, {})
ax1.text(
    ax1.get_xlim()[1] * 0.45, ax1.get_ylim()[1] * 0.75,\
    'Non-linear and/or-monotonic', color='gray',\
    horizontalalignment='center'
)
ax1.text(ax1.get_xlim()[1] * 0.75, ax1.get_ylim()[1] * 0.5,\
    'Almost Monotonic', color='gray', horizontalalignment='center')
ax1.text(ax1.get_xlim()[1] * 0.83, ax1.get_ylim()[1] * 0.2,\
    'Monotonic', color='gray', horizontalalignment='center')
ax1.text(ax1.get_xlim()[1] * 0.9, ax1.get_ylim()[1] * 0.025,
    'Linear', color='gray', horizontalalignment='center') 

前面的代码生成了图 9.14中显示的图表。左边的条形图按µ^对因素进行排序，而每根从条形中伸出的线表示相应的置信区间。右边的协方差图是一个散点图，µ^位于x轴上，σ位于y轴上。因此，点越往右，它就越重要，而它在图中越往上，它与其它因素的交互就越多，就越不单调。自然地，这意味着那些交互不多且主要单调的因素符合线性回归的假设，如线性性和多重共线性。然而，线性与非线性或非单调之间的范围由σ和µ^*的比率对角确定：

图 9.14：表示基本效应的条形图和协方差图

你可以通过前面的协方差图看出，所有因素都是非线性或非单调的。hr无疑是其中最重要的，其次是接下来的两个因素（dow和temp）相对靠近，然后是weather和is_holiday。weather组没有在图中显示，因为交互性结果不确定，但cloud_coverage、rain_1h和snow_1h的交互性比它们单独重要得多。

基本效应帮助我们理解如何根据它们对模型结果的影响来分类我们的因素。然而，这并不是一个稳健的方法来正确量化它们的影响或由因素相互作用产生的影响。为此，我们必须转向使用概率框架分解输出方差并将其追溯到输入的基于方差的全局方法。这些方法包括傅里叶振幅敏感性测试（FAST）和Sobol。我们将在下一节研究后一种方法。

使用因素固定量化不确定性和成本敏感性

使用 Morris 指数，很明显，所有因素都是非线性或非单调的。它们之间有很高的交互性——正如预期的那样！气候因素（temp、rain_1h、snow_1h和cloud_coverage）很可能与hr存在多重共线性。在hr、is_holiday、dow和目标之间也存在一些模式。许多这些因素肯定与目标没有单调关系。我们已经知道了这一点。例如，交通在一天中的小时数增加时并不总是增加。情况对一周中的某一天也是如此！

然而，我们不知道is_holiday和temp对模型的影响程度，尤其是在机组人员的工作时间内，这是一个重要的见解。话虽如此，使用 Morris 指数进行因素优先级排序通常被视为起点或“第一设置”，因为一旦确定存在交互效应，最好是解开它们。为此，有一个“第二设置”，称为因素固定。我们可以量化方差，通过这样做，量化所有因素带来的不确定性。

只有基于方差的方法才能以统计严谨的方式量化这些效应。Sobol 敏感性分析是这些方法之一，这意味着它将模型的输出方差分解成百分比，并将其归因于模型的输入和交互。像 Morris 一样，它有一个采样步骤，以及一个敏感性指数估计步骤。

与 Morris 不同，采样不遵循一系列级别，而是遵循输入数据的分布。它使用准蒙特卡洛方法，在超空间中采样点，这些点遵循输入的概率分布。蒙特卡洛方法是一系列执行随机采样的算法，通常用于优化或模拟。它们寻求在用蛮力或完全确定性的方法无法解决的问题上的捷径。蒙特卡洛方法在敏感性分析中很常见，正是出于这个原因。准蒙特卡洛方法有相同的目标。然而，它们收敛得更快，因为它们使用确定性低偏差序列而不是使用伪随机序列。Sobol 方法使用Sobol 序列，由同一位数学家设计。我们将使用另一种从 Sobol 派生出的采样方案，称为 Saltelli 的。

一旦生成样本，蒙特卡洛估计器就会计算基于方差的敏感性指数。这些指数能够量化非线性非加性效应和第二阶指数，这些指数与两个因素之间的相互作用相关。Morris 可以揭示模型中的交互性，但不能精确地说明它是如何表现的。Sobol 可以告诉你哪些因素在相互作用以及相互作用的程度。

生成和预测 Saltelli 样本

要使用SALib开始 Sobol 敏感性分析，我们必须首先定义一个问题。我们将与 Morris 做同样的事情。这次，我们将减少因素，因为我们意识到weather分组导致了不确定的结果。我们应该包括所有天气因素中最稀疏的；即weather_Clear。由于 Sobol 使用概率框架，将temp、rain_1h和cloud_coverage的范围扩展到它们的最大和最小值是没有害处的，如图9.12所示：

sobol_problem = {
    'num_vars': 8,
    'names': ['dow', 'hr', 'temp', 'rain_1h', 'snow_1h',
              'cloud_coverage', 'is_holiday', 'weather_Clear'],
    'bounds': [
        [0, 4], # dow Monday through Friday
        [-1, 4], # hr
        [-3., 31.], # temp (C)
        [0., 21.], # rain_1h
        [0., 1.6], # snow_1h
        [0., 100.], # cloud_coverage
        [0, 1], # is_holiday
        [0, 1] # weather_Clear
      ],
    'groups': None
} 

生成样本看起来也应该很熟悉。Saltelli 的sample函数需要以下内容：

• 问题陈述（sobol_problem）

• 每个因素要生成的样本数量（300）

• 第二阶索引以进行计算（calc_second_order=True）

由于我们想要交互作用，sample的输出是一个数组，其中每一列代表一个因素，有行（其中N是样本数量，F是因素数量）。我们有八个因素，每个因素有 256 个样本，所以print应该输出 4,608 行和 8 列的形状。首先，我们将像之前一样使用hstack修改它，添加 7 个空因素以进行预测，从而得到 15 列：

saltelli_sample = ss.sample(
    sobol_problem, 256, calc_second_order=True, seed=rand
)
saltelli_sample_mod = np.hstack(
    (saltelli_sample, np.zeros((saltelli_sample.shape[0],7)))
)
print(saltelli_sample_mod.shape) 

现在，让我们对这些样本进行预测。这可能需要一些时间，所以又是咖啡时间：

saltelli_preds = mldatasets.pprox._predict_ts(
    saltelli_sample_mod,
    X_df,
    gen_all,
    lstm_traffic_mdl,
    filt_fn=filt_fn,
    dist_metric=heom_dist.heom,
    lookback=lookback,
    X_scaler=X_scaler,
    y_scaler=y_scaler,
    progress_bar=True
) 

执行 Sobol 敏感性分析

对于 Sobol 敏感性分析（analyze），你所需要的只是一个问题陈述（sobol_problem）和模型输出（saltelli_preds）。但是预测并不能讲述不确定性的故事。当然，预测的交通流量有方差，但只有当交通量超过 1,500 时，这个问题才会出现。不确定性是你想要与风险或回报、成本或收入、损失或利润相关联的东西——一些你可以与你问题相关联的实质性东西。

首先，我们必须评估是否存在任何风险。为了了解样本中的预测交通量是否在工作时间内超过了无建设阈值，我们可以使用print(max(saltelli_preds[:,0]))。最大交通水平应该在 1,800-1,900 左右，这意味着至少存在一些风险，即建筑公司将会支付罚款。我们不必使用预测（saltelli_preds）作为模型的输出，我们可以创建一个简单的二进制数组，当它超过 1,500 时为 1，否则为 0。我们将称之为costs，然后使用它运行analyze函数。注意，这里也设置了calc_second_order=True。如果sample和analyze没有一致的设置，它将抛出一个错误。与 Morris 一样，有一个可选的置信区间水平参数（conf_level），但默认的 0.95 是好的：

costs = np.where(saltelli_preds > 1500, 1,0)[:,0]
factor_fixing_sa = sa.analyze(
    sobol_problem,
    costs,
    calc_second_order=True,
    print_to_console=False
) 

analyze将返回一个包含 Sobol 敏感性指数的字典，包括一阶（S1）、二阶（S2）和总阶（ST）指数，以及总置信区间（ST_conf）。这些指数对应于百分比，但除非模型是加性的，否则总数不一定相加。在表格格式中更容易欣赏这些值，这样我们可以将它们放入 DataFrame 中，并根据总数进行排序和着色，总数可以解释为因素的整体重要性。然而，我们将省略二阶指数，因为它们是二维的，类似于相关图：

sobol_df = pd.DataFrame(
    {
        'features':sobol_problem['names'],
        '1st':factor_fixing_sa['S1'],
        'Total':factor_fixing_sa['ST'],
        'Total Conf':factor_fixing_sa['ST_conf'],
        'Mean of Input':saltelli_sample.mean(axis=0)[:8]
    }
)
sobol_df.sort_values('Total', ascending=False).style
    .background_gradient(cmap='plasma', subset=['Total']) 

上一段代码输出了图 9.15中展示的 DataFrame。你可以看出temp和is_holiday至少在问题定义中指定的边界内排在前面四位。另一个需要注意的事情是weather_Clear确实对其自身有更大的影响，但rain_1h和cloud_coverage似乎对潜在成本没有影响，因为它们的总一阶指数为零：

图 9.15：八个因素的 Sobol 全局敏感性指数

关于一阶值的一些有趣之处在于它们有多低，这表明交互作用占模型输出方差的大部分。我们可以很容易地使用二阶索引来证实这一点。这些索引和一阶索引的组合加起来就是总数：

S2 = factor_fixing_sa['S2']
divnorm = TwoSlopeNorm(vmin=S2.min(), vcenter=0, vmax=S2.max())
sns.heatmap(S2, center=0.00, norm=divnorm, cmap='coolwarm_r',\
            annot=True, fmt ='.2f',\
            xticklabels=sobol_problem['names'],\
            yticklabels=sobol_problem['names']) 

上一段代码输出了图 9.16中的热图：

图 9.16：八个因素的 Sobol 二阶指数

在这里，您可以知道is_holiday和weather_Clear是两个对输出方差贡献最大的因素，其绝对值最高为 0.26。dow和hr与所有因素都有相当大的相互作用。

引入一个现实成本函数

现在，我们可以创建一个成本函数，它接受我们的输入（saltelli_sample）和输出（saltelli_preds），并计算双城将对建筑公司罚款多少，以及额外的交通可能产生的任何额外成本。

如果输入和输出都在同一个数组中，这样做会更好，因为我们需要从两者中获取详细信息来计算成本。我们可以使用hstack将样本及其对应的预测结果连接起来，生成一个包含八个列的数组（saltelli_sample_preds）。然后我们可以定义一个成本函数，它可以计算包含这些九个列的数组的成本（cost_fn）：

#Join input and outputs into a sample+prediction array
saltelli_sample_preds = np.hstack((saltelli_sample, saltelli_preds)) 

我们知道，对于任何样本预测，半容量阈值都没有超过，所以我们甚至不需要在函数中包含每日罚款。除此之外，罚款是每辆超过每小时无施工阈值的车辆 15 美元。除了这些罚款之外，为了能够按时离开，建筑公司估计额外的成本：如果凌晨 4 点超过阈值，额外工资为 1,500 美元，周五额外 4,500 美元以加快设备移动速度，因为周末不能停在高速公路的路肩上。一旦我们有了成本函数，我们就可以遍历组合数组（saltelli_sample_preds），为每个样本计算成本。列表推导可以有效地完成这项工作：

#Define cost function
def cost_fn(x):
    cost = 0
    if x[8] > 1500:
        cost = (x[8] - 1500) * 15
    if round(x[1]) == 4:
        cost = cost + 1500
        if round(x[0]) == 4:
            cost = cost + 4500
    return cost
#Use list comprehension to compute costs for sample+prediction array
costs2 = np.array([cost_fn(xi) for xi in saltelli_sample_preds])
#Print total fines for entire sample predictions
print('Total Fines: $%s' % '{:,.2f}'.format(sum(costs2))) 

print语句应该输出一个介于 17 万美元和 20 万美元之间的成本。但不必担心！建筑队每年只计划在现场工作大约 195 天，每天 5 小时，总共 975 小时。然而，有 4,608 个样本，这意味着由于交通过多，几乎有 5 年的预测成本。无论如何，计算这些成本的目的在于了解它们与模型输入的关系。更多的样本年意味着更紧密的置信区间：

factor_fixing2_sa = sa.analyze(
    sobol_problem, costs2, calc_second_order=True,
    print_to_console=False
) 

现在，我们可以再次进行分析，但使用costs2，并将分析保存到factor_fixing2_sa字典中。最后，我们可以使用这个字典的值生成一个新的排序和彩色编码的 DataFrame，就像我们之前为图 9.15所做的那样，这将生成图 9.17中的输出。

如您从图 9.17中可以看出，一旦实际成本被考虑在内，dow、hr和is_holiday成为更具风险的因素，而与图 9.15相比，snow_1h和temp变得不那么相关：

图 9.17：使用现实成本函数计算八个因素的 Sobol 全局敏感性指数

用表格难以欣赏的是敏感性指数的置信区间。为此，我们可以使用条形图，但首先，我们必须将整个字典转换成一个 DataFrame，以便SALib的绘图函数可以绘制它：

factor_fixing2_df = factor_fixing2_sa.to_df()
fig, (ax) = plt.subplots(1,1, figsize=(15, 7))
sp.plot(factor_fixing2_df[0], ax=ax) 

前面的代码生成了图 9.18中的条形图。dow的 95%置信区间比其他重要因素大得多，考虑到一周中各天之间的差异很大，这并不令人惊讶。另一个有趣的见解是weather_Clear具有负一阶效应，因此正的总阶指数完全归因于二阶指数，这扩大了置信区间：

图 9.18：使用现实成本函数绘制的条形图，包含 Sobol 敏感性总阶指数及其置信区间

要了解如何，让我们再次绘制图 9.16所示的散点图，但这次使用factor_fixing2_sa而不是factor_fixing_sa。图 9.19中的散点图应该描绘出模型中成本的现实反映：

图 9.19：在考虑更现实的成本函数时，七个因素的 Sobol 二阶指数

前面的散点图显示了与图 9.16中相似的显著交互，但由于有更多的阴影，它们更加细腻。很明显，weather_Clear与is_holiday结合时具有放大作用，而对dow和hr则有调和作用。

任务完成

任务是训练一个交通预测模型，并了解哪些因素会创造不确定性，并可能增加建筑公司的成本。我们可以得出结论，潜在的 35,000 美元/年的罚款中有很大一部分可以归因于is_holiday因素。因此，建筑公司应该重新考虑工作假日。三月至十一月之间只有七个或八个假日，由于罚款，它们可能比在几个星期日工作成本更高。考虑到这个警告，任务已经成功，但仍有很多改进的空间。

当然，这些结论是针对LSTM_traffic_168_compact1模型——我们可以将其与其他模型进行比较。尝试将笔记本开头的model_name替换为LSTM_traffic_168_compact2，这是一个同样小巧但显著更稳健的模型，或者LSTM_traffic_168_optimal，这是一个更大但表现略好的模型，并重新运行笔记本。或者浏览名为Traffic_compact2和Traffic_optimal的笔记本，这些笔记本已经使用相应的模型重新运行。你会发现，可以训练和选择能够更好地管理不确定输入的模型。话虽如此，改进并不总是通过简单地选择更好的模型就能实现。

例如，可以进一步深入探讨的是temp、rain_1h和snow_1h的真正影响。我们的预测近似方法排除了 Sobol 测试极端天气事件的影响。如果我们修改模型以在单个时间步长上训练聚合的天气特征，并内置一些安全措施，我们就可以使用 Sobol 模拟天气极端情况。而且，敏感性分析的“第三设置”，即因素映射，可以帮助精确指出某些因素值如何影响预测结果，从而进行更稳健的成本效益分析，但这一点我们不会在本章中涉及。

在本书的第二部分，我们探讨了多种解释方法的生态系统：全局和局部；针对特定模型和非特定模型；基于排列和基于敏感度的。对于任何机器学习用例，可供选择的方法并不缺乏。然而，必须强调的是，没有任何方法是完美的。尽管如此，它们可以相互补充，以更接近地理解您的机器学习解决方案及其旨在解决的问题。

本章关注预测中的确定性，旨在揭示机器学习社区中的一个特定问题：过度自信。在“可解释性的商业案例”部分的第一章，《解释、可解释性、可解释性；以及这一切为什么都重要？》，描述了充斥在人类决策中的许多偏见。这些偏见通常是由对领域知识或我们模型令人印象深刻的成果的过度自信所驱动的。而这些令人印象深刻的成果使我们无法理解我们模型的局限性，随着公众对 AI 的不信任增加，这一点变得更加明显。

正如我们在第一章中讨论的，解释、可解释性、可解释性；以及为什么这一切都很重要？，机器学习仅用于解决不完整问题。否则，我们不如使用在闭环系统中发现的确定性程序编程。解决不完整问题的最佳方法是一个不完整的解决方案，它应该被优化以尽可能多地解决它。无论是通过梯度下降、最小二乘估计还是分割和修剪决策树，机器学习不会产生一个完美泛化的模型。机器学习中的这种不完整性正是我们需要解释方法的原因。简而言之：模型从我们的数据中学习，我们可以从我们的模型中学到很多，但只有当我们解释它们时才能做到！

然而，可解释性并不止于此。模型解释可以驱动决策并帮助我们理解模型的优势和劣势。然而，数据或模型本身的问题有时会使它们变得难以解释。在本书的第三部分中，我们将学习如何通过降低复杂性、减轻偏差、设置护栏和增强可靠性来调整模型和训练数据以提高可解释性。

统计学家 George E.P. Box 曾著名地开玩笑说，“所有模型都是错误的，但有些是有用的。”也许它们并不总是错误的，但机器学习从业者需要谦逊地接受，即使是高性能模型也应受到审查，并且我们对它们的假设也应受到质疑。机器学习模型的不确定性是可以预期的，不应该是羞耻或尴尬的来源。这使我们得出本章的另一个结论：不确定性伴随着后果，无论是成本还是利润提升，我们可以通过敏感性分析来衡量这些后果。

摘要

阅读本章后，你应该了解如何评估时间序列模型的预测性能，知道如何使用集成梯度对他们进行局部解释，以及如何使用 SHAP 产生局部和全局归因。你还应该知道如何利用敏感性分析因子优先级和因子固定来优化任何模型。

在下一章中，我们将学习如何通过特征选择和工程来降低模型的复杂性，使其更具可解释性。

数据集和图像来源

• TomTom，2019 年，交通指数：nonews.co/wp-content/uploads/2020/02/TomTom2019.pdf

• UCI 机器学习仓库，2019 年，都市州际交通流量数据集：archive.ics.uci.edu/ml/datasets/Metro+Interstate+Traffic+Volume

进一步阅读

• Wilson, D.R. 和 Martinez, T.，1997 年，改进的异构距离函数。J. Artif. Int. Res. 6-1. 第 1-34 页：arxiv.org/abs/cs/9701101

• Morris, M., 1991, 《初步计算实验的因子抽样计划》. Quality Engineering, 37, 307-310: doi.org/10.2307%2F1269043

• Saltelli, A., Tarantola, S., Campolongo, F., and Ratto, M., 2007, 《实践中的敏感性分析：评估科学模型指南》. Chichester: John Wiley & Sons.

• Sobol, I.M., 2001, 《非线性数学模型的全球敏感性指数及其蒙特卡洛估计》. MATH COMPUT SIMULAT, 55(1–3), 271-280: doi.org/10.1016/S0378-4754(00)00270-6

• Saltelli, A., P. Annoni, I. Azzini, F. Campolongo, M. Ratto, and S. Tarantola, 2010, 《模型输出的方差敏感性分析：总敏感性指数的设计和估计器》. Computer Physics Communications, 181(2):259-270: doi.org/10.1016/j.cpc.2009.09.018

在 Discord 上了解更多

要加入这本书的 Discord 社区——在那里您可以分享反馈、向作者提问，并了解新书发布——请扫描下面的二维码：

packt.link/inml