Kaggle-笔记本开发指南-二-

Kaggle 笔记本开发指南（二）

原文：annas-archive.org/md5/24dfa97e36ead23596f0ef3b74ce3f05

译者：飞龙

协议：CC BY-NC-SA 4.0

第六章：你能预测蜜蜂的亚种吗？

在本章中，我们将学习如何处理图像数据并开始构建用于图像分类的模型。多年来，计算机视觉在数据科学和数据分析中的应用呈指数增长。在 Kaggle 上一些最引人注目的（拥有大量点赞和复制的，例如复制和编辑）笔记本并不是 探索性数据分析（EDA）笔记本，而是构建模型的笔记本。

在本章中，我们将演示如何使用您深入的数据分析来准备构建模型，并且我们还将向您介绍模型迭代优化过程的一些见解。这不仅仅是为了比赛，而是为了一个图像数据集。数据集是 BeeImage Dataset: Annotated Honey Bee Images（参见 参考文献 1）。在前一章中，我们也开始使用 Plotly 作为可视化库。在本章中，我们将继续使用 Plotly 来可视化数据集特征。我们将一些有用的可视化函数与 Plotly 一起放在一个实用脚本中，名为 plotly-utils（参见 参考文献 2）。与本章相关的笔记本是 Honeybee Subspecies Classification（参见 参考文献 3）。

本章将涵盖以下主题：

• 对 BeeImage Dataset: Annotated Honey Bee Images 的全面数据探索。

• 在准备模型基线之后，逐步优化模型，分析对训练和验证指标变化的影响，并采取新措施进一步改进模型。

数据探索

BeeImage Dataset: Annotated Honey Bee Images 包含一个 逗号分隔格式 (.csv) 文件，bee_data.csv，包含 5172 行和 9 列，以及一个包含 5172 张图片的文件夹：

图 6.1：bee_data.csv 数据文件的样本

正如您所看到的，前面的数据框包含以下列：

• file: 图片文件名

• date: 拍摄图片的日期

• time: 拍摄图片的时间

• location: 美国位置，包括城市、州和国家名称

• zip code: 与位置相关的邮政编码

• subspecies: 当前图像中蜜蜂所属的亚种

• health: 当前图像中蜜蜂的健康状态

• pollen_carrying: 表示图片中蜜蜂是否带有花粉附着在其腿上

• caste: 蜜蜂的社会阶层

我们将开始数据探索之旅，进行一些质量检查，重点关注 bee_data.csv 文件，然后是图像。对于数据质量检查，我们将使用在 第四章 中之前介绍的一个实用脚本，data_quality_stats。

数据质量检查

如下所示，数据集没有任何缺失值。所有特征都是 string 类型。

图 6.2：bee_data.csv 文件中的缺失值。结果使用 data_quality_stats 函数获得

在图 6.3中，我们展示了数据集特征的唯一值。数据是在以下情况下收集的：

• 在 6 个不同的日期和 35 个不同的时间

• 在 8 个地点，7 个不同的邮编

在数据中，有七个亚种，用六个不同的健康问题表示。

图 6.3：bee_data.csv 文件中的唯一值。结果使用 data_quality_stats 函数获得

从图 6.4所示的数据中可以看出，21%的图像来自单一日期（16 个不同日期中的一个）。曾经有 11%的图像是在某个时间收集的。有一个单一的位置（加利福尼亚州的萨拉托加，邮编 95070），在那里收集了 2000 张（或 39%）的图像。意大利蜜蜂是最常见的物种。65%的图像代表健康的蜜蜂。几乎所有图像都显示了不带花粉的蜜蜂，而且所有图像都来自工蜂群体。

图 6.4：bee_data.csv 文件中最频繁的值。结果使用 data_quality_stats 函数获得

接下来，我们将与bee_data.csv中的特征并行回顾图像数据。我们还将介绍读取和可视化图像的函数。

探索图像数据

首先，我们检查数据集中存在的所有图像名称是否也存在于图像文件夹中：

file_names = list(honey_bee_df['file'])
print("Matching image names: {}".format(len(set(file_names).intersection(image_files)))) 

结果是，所有在.csv文件中索引的图像都存在于images文件夹中。接下来，我们检查图像大小。为此，我们可以使用以下代码读取图像：

def read_image_sizes(file_name):
    """
    Read images size using skimage.io
    Args:
        file_name: the name of the image file
    Returns:
        A list with images shape
    """
    image = skimage.io.imread(config['image_path'] + file_name)
    return list(image.shape) 

或者，我们可以使用以下代码根据 OpenCSV（cv2）库读取图像：

def read_image_sizes_cv(file_name):
    """
    Read images size using OpenCV
    Args:
        file_name: the name of the image file
    Returns:
        A list with images shape
    """
    image = cv2.imread(config['image_path'] + file_name)
    return list(image.shape) 
skimage.io:

%timeit m = np.stack(subset.apply(read_image_sizes)) 

下面的代码用于使用基于 opencv 的方法测量读取图像的执行时间：

%timeit m = np.stack(subset.apply(read_image_sizes_cv)) 

比较显示，使用基于opencv的方法执行更快：

• 使用skimage.io：

129 ms ± 4.12 ms 每循环（7 次运行的平均值±标准差，每次循环 1 次）

• 使用opencv：

127 ms ± 6.79 ms 每循环（7 次运行的平均值±标准差，每次循环 10 次）

然后，我们应用最快的方法提取每个图像的形状（宽度、高度和深度，或颜色维度的数量）并将其添加到每个图像的数据集中：

t_start = time.time()
m = np.stack(honey_bee_df['file'].apply(read_image_sizes_cv))
df = pd.DataFrame(m,columns=['w','h','c'])
honey_bee_df = pd.concat([honey_bee_df,df],axis=1, sort=False)
t_end = time.time()
print(f"Total processing time (using OpenCV): {round(t_end-t_start, 2)} sec.") 

执行前面代码的输出是：

Total processing time (using OpenCV): 34.38 sec. 
boxplot. In the first, we show the image width distribution, and in the second, the image height distribution. The boxplot shows the minimum, first quartile, median, third quartile, and maximum values in the distribution of the value we plot. We also show the outliers’ values as points on each of the traces:

traceW = go.Box(
    x = honey_bee_df['w'],
    name="Width",
     marker=dict(
                color='rgba(238,23,11,0.5)',
                line=dict(
                    color='red',
                    width=1.2),
            ),
    orientation='h')
traceH = go.Box(
    x = honey_bee_df['h'],
    name="Height",
    marker=dict(
                color='rgba(11,23,245,0.5)',
                line=dict(
                    color='blue',
                    width=1.2),
            ),
    orientation='h')
data = [traceW, traceH]
layout = dict(title = 'Width & Heights of images',
          xaxis = dict(title = 'Size', showticklabels=True), 
          yaxis = dict(title = 'Image dimmension'),
          hovermode = 'closest',
         )
fig = dict(data=data, layout=layout)
iplot(fig, filename='width-height') 

结果绘制在图 6.5中。宽度和高度的中间值分别为 61 和 62。宽度和高度都有许多异常值（宽度最大值为 520，高度最大值为 392）。

图 6.5：图像的宽度和高度分布

在我们的分析中，我们包括了数据集中的所有特征，而不仅仅是与图像相关的特征。在我们开始构建预测模型的基线之前，我们希望了解与蜜蜂图像数据集：标注的蜜蜂图像相关的所有方面。

位置

通过按拍摄图片的位置和 ZIP 代码对数据集中的数据进行分组，我们可以观察到有一个位置具有相同的 ZIP 代码和类似的名字：

图 6.6：拍摄带有蜜蜂的图片的位置和 ZIP 代码

我们可以观察到，美国乔治亚州的雅典出现了两个略有不同的名称。我们只是使用以下代码将它们合并：

honey_bee_df = honey_bee_df.replace({'location':'Athens, Georgia, USA'}, 'Athens, GA, USA') 

现在，让我们使用 Plotly 实用脚本模块中的一个函数来可视化结果位置数据的分布：

tmp = honey_bee_df.groupby(['zip code'])['location'].value_counts()
df = pd.DataFrame(data={'Images': tmp.values}, index=tmp.index).reset_index()
df['code'] = df['location'].map(lambda x: x.split(',', 2)[1])
plotly_barplot(df, 'location', 'Images', 'Tomato', 'Locations', 'Number of images', 'Number of bees images per location') 

函数plotly_barplot的代码如下：

def plotly_barplot(df, x_feature, y_feature, col, x_label, y_label, title):
    """
    Plot a barplot with number of y for category x
    Args:
        df: dataframe
        x_feature: x feature
        y_feature: y feature
        col: color for markers
        x_label: x label
        y_label: y label
        title: title

    Returns:
        None
    """
    trace = go.Bar(
            x = df[x_feature],
            y = df[y_feature],
            marker=dict(color=col),
            #text=df['location']
        )
    data = [trace]
    layout = dict(title = title,
              xaxis = dict(title = x_label, showticklabels=True, tickangle=15), 
              yaxis = dict(title = y_label),
              hovermode = 'closest'
             )
    fig = dict(data = data, layout = layout)
    iplot(fig, filename=f'images-{x_feature}-{y_feature}') 

在图 6.7中，我们展示了拍摄蜜蜂图像的位置分布。大多数图像来自加利福尼亚州的萨拉托加（2000 张图像），其次是乔治亚州的雅典和爱荷华州的迪莫因。

图 6.7：位置分布

我们还基于一个选定的标准构建了一个用于可视化图像子集的函数。以下代码是根据位置选择图像并显示它们的子集（一行五张，来自同一位置）：

#list of locations
locations = (honey_bee_df.groupby(['location'])['location'].nunique()).index
def draw_category_images(var,cols=5):
    categories = (honey_bee_df.groupby([var])[var].nunique()).index
    f, ax = plt.subplots(nrows=len(categories),ncols=cols, figsize=(2*cols,2*len(categories)))
    # draw a number of images for each location
for i, cat in enumerate(categories):
        sample = honey_bee_df[honey_bee_df[var]==cat].sample(cols)
        for j in range(0,cols):
            file=config['image_path'] + sample.iloc[j]['file']
            im=imageio.imread(file)
            ax[i, j].imshow(im, resample=True)
            ax[i, j].set_title(cat, fontsize=9)  
    plt.tight_layout()
    plt.show() 

在图 6.8中，展示了这个选择的一部分（仅限于前两个位置）。完整的图像可以在相关的笔记本中查看：

图 6.8：来自两个地点的蜜蜂图像（从完整图片中选择，使用前面的代码获取）

日期和时间

让我们继续详细分析我们数据集中的特征。我们现在开始分析date和time数据。我们将date转换为datetime并提取年、月和日。我们还转换time并提取小时和分钟：

honey_bee_df['date_time'] = pd.to_datetime(honey_bee_df['date'] + ' ' + honey_bee_df['time'])
honey_bee_df["year"] = honey_bee_df['date_time'].dt.year
honey_bee_df["month"] = honey_bee_df['date_time'].dt.month
honey_bee_df["day"] = honey_bee_df['date_time'].dt.day
honey_bee_df["hour"] = honey_bee_df['date_time'].dt.hour
honey_bee_df["minute"] = honey_bee_df['date_time'].dt.minute 

在图 6.9中展示了每天和大约的小时及位置的蜜蜂图像数量的可视化。这个可视化的代码首先通过date_time和hour对数据进行分组，并计算每个日期和一天中的时间收集到的图像数量：

tmp = honey_bee_df.groupby(['date_time', 'hour'])['location'].value_counts()
df = pd.DataFrame(data={'Images': tmp.values}, index=tmp.index).reset_index() 

然后，我们构建当我们在图中的某个点上悬停时显示的文本。这个文本将包括小时、位置和图像数量。然后我们将悬停文本作为新数据集中的一个新列添加：

hover_text = []
for index, row in df.iterrows():
    hover_text.append(('Date/time: {}<br>'+
                      'Hour: {}<br>'+
                      'Location: {}<br>'+
                      'Images: {}').format(row['date_time'],
                                            row['hour'],
                                            row['location'],
                                            row['Images']))
df['hover_text'] = hover_text 

然后，我们为每个位置绘制一个散点图，表示图片收集的时间和小时。每个点的尺寸与在该位置、一天中的某个时间点和某个日期拍摄的图像数量成比例：

locations = (honey_bee_df.groupby(['location'])['location'].nunique()).index
data = []
for location in locations:
    dfL = df[df['location']==location]
    trace = go.Scatter(
        x = dfL['date_time'],y = dfL['hour'],
        name=location,
        marker=dict(
            symbol='circle',
            sizemode='area',
            sizeref=0.2,
            size=dfL['Images'],
            line=dict(
                width=2
            ),),
        mode = "markers",
        text=dfL['hover_text'],
    )
    data.append(trace)

layout = dict(title = 'Number of bees images per date, approx. hour and location',
          xaxis = dict(title = 'Date', showticklabels=True), 
          yaxis = dict(title = 'Hour'),
          hovermode = 'closest'
         )
fig = dict(data = data, layout = layout)
iplot(fig, filename='images-date_time') 

在下一张图像，图 6.9中，我们看到运行上述代码的结果。大多数图片是在八月份拍摄的。大多数图片也是在下午时段拍摄的。

图 6.9：每天和大约的小时及位置的蜜蜂图像数量

亚种

我们使用相同的 plotly_barplot 函数来可视化亚种的分布。大多数蜜蜂是意大利蜜蜂，其次是俄罗斯蜜蜂和卡尼奥兰蜜蜂（见 图 6.10）。其中 428 张图像未被分类（标签值为**-1**）。我们将未分类的图像归为一个亚种类别。

图 6.10：按日期和大约的小时及位置划分的蜜蜂图像数量

在 图 6.11 中，我们展示了一些图像的选择，其中只包含少数亚种样本：

图 6.11：几个亚种蜜蜂图像的样本

现在，让我们表示每个亚种和位置的图像数量，以及每个亚种和小时的图像数量（见 图 6.12）。收集到的图像数量最多的是来自加利福尼亚州萨拉托加（1972 张图像），所有图像都是意大利蜜蜂。在 13 点收集到的图像数量最多，所有图像也都是意大利蜜蜂（909 张图像）。

图 6.12：按亚种和位置划分的图像数量（上方）以及按亚种和小时划分的图像数量（下方）

Subspecies 图像在重量和高度上具有很大的多样性。图 6.13 使用箱线图展示了重量和高度的分布。

图 6.13：每个亚种的图像大小分布 – 宽度（上方）和高度（下方）

VSH 意大利蜜蜂在宽度和高度上都有最大的平均值和最大的方差。西方蜜蜂、卡尼奥兰蜜蜂和混合本地种群 2在重量和高度上的分布最为紧凑（方差较低）。数量最多的亚种意大利蜜蜂显示出较小的中位数和较大的方差，有很多异常值。在下图中，我们将在同一个散点图上展示重量和高度分布：

图 6.14：每个亚种的图像大小分布 – 散点图

前面的图示使用散点图展示了重量和高度分布，下面的代码展示了这一可视化的实现。首先，我们定义一个函数来绘制散点图，其中图像宽度在 x 轴上，图像高度在 y 轴上：

def draw_trace_scatter(dataset, subspecies):
    dfS = dataset[dataset['subspecies']==subspecies];
    trace = go.Scatter(
        x = dfS['w'],y = dfS['h'],
        name=subspecies,
        mode = "markers",
        marker = dict(opacity=0.8),
        text=dfS['subspecies'], 
    )
    return trace 

我们现在使用上面定义的函数为每个亚种绘制散点图。每次函数调用都会创建一个轨迹，我们将这些轨迹添加到 Plotly 图中：

subspecies = (honey_bee_df.groupby(['subspecies'])['subspecies'].nunique()).index
def draw_group(dataset, title,height=600):
    data = list()
    for subs in subspecies:
        data.append(draw_trace_scatter(dataset, subs))

    layout = dict(title = title,
              xaxis = dict(title = 'Width',showticklabels=True),
              yaxis = dict(title = 'Height', showticklabels=True, tickfont=dict(
                family='Old Standard TT, serif',
                size=8,
                color='black'),), 
              hovermode = 'closest',
              showlegend=True,
                  width=800,
                  height=height,
             )
    fig = dict(data=data, layout=layout)
    iplot(fig, filename='subspecies-image')
draw_group(honey_bee_df,  "Width and height of images per subspecies") 

健康

图 6.15 展示了具有各种健康问题的图像分布。大多数图像是健康蜜蜂（3384），其次是少量瓦拉罗，蜂箱甲虫（579），瓦拉罗，小蜂箱甲虫（472），以及蚂蚁问题（457）：

图 6.15：不同健康问题的图像数量

如果我们分析每个亚种和健康问题的图像数量（见图 6.16），我们可以观察到只有少量健康和亚种值组合存在。大多数图像是健康的意大利蜜蜂（1972），其次是少量瓦螨、蜂箱甲虫，然后是意大利蜜蜂（579），最后是健康的俄罗斯蜜蜂（527）。未知亚种要么是健康的（177）要么是蜂群被盗（251）。

图 6.16：每个亚种和不同健康问题的蜜蜂图像数量

在图 6.17中，我们绘制了每个地点、亚种和健康问题的图像数量：

图 6.17：每个地点、亚种和健康问题的图像数量

其他

携带花粉的蜜蜂图像数量很少。图 6.18显示了一些携带花粉和不携带花粉的蜜蜂图像。所有蜜蜂都来自同一个等级：工蜂等级。

图 6.18：携带和不携带花粉的蜜蜂图像的选择

结论

我们使用plotly_sankey脚本从plotly_utils实用脚本模块中绘制的桑基图，来绘制图 6.19中的摘要图。桑基图主要用于可视化流程或流动，例如，在经济学中能源的生产及其来源和消费者。我这里用它来达到另一个目的，即总结具有多个特征的数据分布。它显示了同一图表中按日期、时间、地点、ZIP 代码、亚种和健康分布的图像。由于空间限制，这里没有给出桑基图的适配代码（请参考参考文献 2获取与本章相关的代码示例）；我们只包含了适配蜜蜂数据以使用此功能的代码：

tmp = honey_bee_df.groupby(['location', 'zip code', 'date', 'time', 'health'])['subspecies'].value_counts()
df = pd.DataFrame(data={'Images': tmp.values}, index=tmp.index).reset_index()
fig = plotly_sankey(df,cat_cols=['date', 'time', 'location', 'zip code', 'subspecies', 'health'],value_cols='Images',
                    title='Honeybee Images: date | time | location | zip code | subspecies | health', 
                 color_palette=[ "darkgreen", "lightgreen", "green", "gold", "black", "yellow"],
           height=800)
iplot(fig, filename='Honeybee Images') 

图 6.19中的可视化，一个漏斗形图，使我们能够在一个单一的图表中捕捉多个特征之间的关系：

图 6.19：图像摘要

到目前为止，我们分析了数据集中特征分布。现在，我们对数据集中的数据有了更好的理解。在本章接下来的部分，我们将开始准备构建一个机器学习模型，以对亚种进行图像分类，这是本章的第二大且更为重要的目标。

亚种分类

本节的目标将是使用迄今为止调查的图像构建一个机器学习模型，该模型可以正确预测亚种。由于我们只有一个数据集，我们将首先将数据分割成三个子集：用于训练、验证和测试数据。我们将在训练过程中使用训练和验证数据：训练数据用于向模型提供数据，验证数据用于验证模型如何使用新数据预测类别（即亚种）。然后，训练和验证后的模型将用于预测测试集中的类别，该类别既未用于训练也未用于验证。

数据分割

首先，我们将数据分割成train和test，使用 80%–20%的分割。然后，再次将train数据分割成训练和验证，使用相同的 80%–20%分割。分割使用stratify和subspecies作为参数执行，确保平衡的子集，同时尊重训练、验证和测试子样本集中类的整体分布。这里选择的训练/验证/测试分割的百分比是任意选择的，并不是研究或优化的结果。在您的实验中，您可以处理不同的训练/验证/测试子集值，也可以选择不使用stratify：

train_df, test_df = train_test_split(honey_bee_df, test_size=config['test_size'], random_state=config['random_state'], 
             stratify=honey_bee_df['subspecies'])
train_df, val_df = train_test_split(train_df, test_size=config['val_size'], random_state=config['random_state'], 
             stratify=train_df['subspecies']) 

最终，我们将有三个子集，如下所示：

• 训练集行数：3309

• 验证集行数：828

• 测试集行数：1035

我们将图像分割成子集，对应于图像名称的子集。我们创建了读取图像并将它们全部调整到配置中定义的相同维度的函数，使用skimage.io和opencv。我们决定将所有图像调整到 100 x 100 像素。我们的决定是基于对图像尺寸分布的分析。您可以选择修改笔记本中提供的代码（参考 3）并尝试不同的图像尺寸。

以下代码使用skimage.io读取图像并根据配置中设置的尺寸调整大小。您可以更改配置并调整图像大小，使用不同的图像高度和宽度值：

def read_image(file_name):
    """
    Read and resize the image to image_width x image_height
    Args:
        file_name: file name for current image
    Returns:
        resized image
    """
    image = skimage.io.imread(config['image_path'] + file_name)
    image = skimage.transform.resize(image, (config['image_width'], config['image_height']), mode='reflect')
    return image[:,:,:config['image_channels']] 

下面的代码使用 OpenCV 读取并调整图像大小。该函数与上面展示的之前的函数不同之处仅在于读取图像文件的方法：

def read_image_cv(file_name):
    """
    Read and resize the image to image_width x image_height
    Args:
        file_name: file name for current image
    Returns:
        resized image
    """
    image = cv2.imread(config['image_path'] + file_name)
    image = cv2.resize(image, (config['image_width'], config['image_height']))
    return image[:,:,:config['image_channels']] 

然后，我们将这些函数应用于所有数据集文件，以读取和调整数据集中的图像大小。

我们还创建了与分类目标变量对应的虚拟变量。我们更喜欢使用这种方法，因为我们将为多类分类准备一个模型，该模型为每个类别输出概率：

def categories_encoder(dataset, var='subspecies'):
    X = np.stack(dataset['file'].apply(read_image))
    y = pd.get_dummies(dataset[var], drop_first=False)
    return X, y
s_time = time.time()
X_train, y_train = categories_encoder(train_df)
X_val, y_val = categories_encoder(val_df)
X_test, y_test = categories_encoder(test_df)
e_time = time.time()
print(f"Total time: {round(e_time-s_time, 2)} sec.") 

通过这一点，我们已经展示了如何读取和调整我们的图像大小。接下来，我们将看到如何通过乘以我们的图像来增加训练集中的数据量，以便向模型展示更多种类的数据。

数据增强

我们将使用深度学习模型来对图像中的亚种进行分类。通常，深度学习模型在训练数据量较大时表现更好。使用数据增强，我们还创建了更多样化的数据，这对模型质量也有益。如果我们在训练过程中让模型接触到更多样化的数据，我们的模型将提高其泛化能力。

我们基于keras.preprocessing.image中的ImageDataGenerator定义了一个数据增强组件。在本笔记本中，我们将使用 Keras 构建模型的各个组件。ImageDataGenerator组件通过应用以下参数初始化，以创建训练数据集的随机变化：

• 对原始图像进行旋转（0 到 180 度范围内）

• 缩放（10%）

• 水平和垂直方向上的平移（10%）

• 水平和垂直方向的平移（10%）

这些变化可以分别控制。并非所有用例都允许或从应用上述所有转换中受益（例如，考虑具有建筑或其他地标图像的情况，对于这些图像，旋转并不太有意义）。以下代码可以用于我们的情况来初始化和拟合图像生成器：

image_generator = ImageDataGenerator(
        featurewise_center=False,
        samplewise_center=False,
        featurewise_std_normalization=False,
        samplewise_std_normalization=False,
        zca_whitening=False,
        rotation_range=180,
        zoom_range = 0.1, 
        width_shift_range=0.1,
        height_shift_range=0.1, 
        horizontal_flip=True,
        vertical_flip=True)
image_generator.fit(X_train) 

然后，我们继续构建和训练基线模型。

构建基线模型

几乎总是建议您从一个简单的模型开始，然后进行错误分析。根据错误分析，您将需要进一步细化您的模型。例如，如果您观察到您的基线模型在训练中获得了很大的误差，您需要首先改进训练。您可以通过添加更多数据、改进您的数据标注或创建更好的特征来实现这一点。如果您的训练误差很小，但您有较高的验证误差，这意味着您的模型可能过度拟合了训练数据。在这种情况下，您需要尝试提高模型泛化能力。您可以尝试各种技术来提高模型泛化能力。关于此类分析，请参阅本章末尾的参考 4。

我们将使用Keras库来定义我们的模型。Keras（见参考 5）是 TensorFlow 机器学习平台（见参考 6）的包装器。它允许您通过定义具有专用层的顺序结构来创建强大的深度学习模型。我们将向我们的模型添加以下层：

• 一个具有 3 维度的 16 个滤波器的Conv2D层

• 一个具有 2 倍缩减因子的MaxPooling2D层

• 一个具有 3 维度的 16 个滤波器的卷积层

• 一个Flatten层

• 一个具有亚种目标特征类别数量的Dense层

上述架构是一个非常简单的卷积神经网络的例子。convolutional2d层的作用是对 2D 输入应用滑动卷积滤波器。maxpool2d层将通过在输入窗口（参见参考文献 5获取更多详细信息）上取最大值来沿其空间维度（宽度和高度）对输入进行下采样。

构建所述架构的代码如下：

model1=Sequential()
model1.add(Conv2D(config['conv_2d_dim_1'], 
                  kernel_size=config['kernel_size'], 
                  input_shape=(config['image_width'], config['image_height'],config['image_channels']), 
                  activation='relu', padding='same'))
model1.add(MaxPool2D(config['max_pool_dim']))
model1.add(Conv2D(config['conv_2d_dim_2'],  kernel_size=config['kernel_size'],
                  activation='relu', padding='same'))
model1.add(Flatten())
model1.add(Dense(y_train.columns.size, activation='softmax'))
model1.compile(optimizer='adam', loss='categorical_crossentropy', metrics=['accuracy']) 

在图 6.20中，我们展示了模型的摘要信息。如您所见，可训练参数的总数是 282,775：

图 6.20：基准模型的摘要

对于基准模型，我们从一个小模型开始，并减少 epoch 的数量进行训练。输入图像的大小是 100 x 100 x 3（如我们之前解释的）。我们将对这个模型进行五次 epoch 的训练。批大小设置为 32。运行训练的代码如下：

train_model1  = model1.fit_generator(image_generator.flow(X_train, y_train, batch_size=config['batch_size']),
                        epochs=config['no_epochs_1'],
                        validation_data=[X_val, y_val],
                        steps_per_epoch=len(X_train)/config['batch_size']) 

图 6.21显示了基准模型的训练日志。在训练过程中，我们没有保存最佳模型版本；最后一步的模型权重将用于测试。

图 6.21：基准模型的训练日志。显示了每个步骤的训练损失和准确率以及验证损失和准确率。

在每个批次之后更新训练损失和准确率，在每个 epoch 结束时计算验证损失和准确率。接下来，在模型训练和验证之后，我们将评估测试集的损失和准确率：

score = model1.evaluate(X_test, y_test, verbose=0)
print('Test loss:', score[0])
print('Test accuracy:', score[1]) 

在图 6.22中，我们展示了训练和验证损失以及训练和验证准确率：

图 6.22：基准模型 – 训练和验证准确率（左）和训练和验证损失（右）

得到的结果如下：

• 测试损失：0.42

• 测试准确率：0.82

测试损失指的是损失函数，这是一个数学函数，用于衡量预测值与真实值之间的差异。在训练过程中，通过测量这个值，对于训练集和验证集，我们可以监控模型的学习和预测的改进情况。

使用sklearn中的metrics.classification_report指标分类报告，我们计算了训练数据中每个类的精确度、召回率、f1 分数和准确率。相应的代码如下：

def test_accuracy_report(model):
    predicted = model.predict(X_test)
    test_predicted = np.argmax(predicted, axis=1)
    test_truth = np.argmax(y_test.values, axis=1)
    print(metrics.classification_report(test_truth, test_predicted, target_names=y_test.columns)) 
    test_res = model.evaluate(X_test, y_test.values, verbose=0)
    print('Loss function: %s, accuracy:' % test_res[0], test_res[1]) 

在图 6.23中，我们展示了测试集的分类报告，其中我们使用了基于训练数据的基线模型拟合。精确度、召回率和 f1 分数的宏平均分别为0.78、0.72和0.74（支持数据为1035）。加权平均精确度、召回率和 f1 分数分别为0.82、0.83和0.82。

图 6.23：使用基线模型对测试数据的分类报告

这些加权平均分数较高，因为与所有类别分数的简单平均值不同，这些是加权平均值，所以与更好表示的类别相关的较高分数将对整体平均有更高的贡献。1 Mixed local stock 2（0.52）和VSH 意大利蜜蜂（0.68）的精确度/类别分数是最差的。最差的整体分数是VSH 意大利蜜蜂，其中召回率为 0.33。

逐步优化模型

如果我们现在回到训练和验证错误，我们可以看到验证和训练准确率大约是 0.81 和 0.82。

我们将继续训练模型，为了避免过拟合，我们还将引入两个Dropout层，每个层的系数为 0.4。Dropout层在神经网络中用作正则化方法。其目的是防止过拟合并提高模型的泛化能力。作为参数给出的系数是在每个训练 epoch 中随机选择的输入百分比，将其设置为 0。模型的结构在图 6.24中描述。可训练参数的数量将保持不变。

图 6.24：优化模型摘要。添加了两个 Dropout 层

我们还将 epoch 的数量扩展到 10。让我们看看图 6.25中的结果，其中我们展示了训练和验证准确率以及训练和验证损失。

图 6.25：优化后的模型（版本 2）- 训练和验证准确率（左）和训练和验证损失（右）

最终训练损失为 0.32，最终训练准确率为 0.87。最终验证损失为 0.28，最终验证准确率为 0.88。这些都是改进的结果。当然，训练准确率主要是由于我们训练了更多的 epoch。

图 6.26：使用第二次优化的模型（训练 epoch 增加到 10 并添加 Dropout 层）的测试数据分类报告

根据验证准确度，结果是更多的训练周期以及添加Dropout层的结果，这些层保持了过拟合的控制。现在让我们检查测试损失和准确度，以及查看测试数据的整个分类报告。

宏平均和加权平均的度量标准在精确度、召回率和 f1 分数上都有所提高。我们还可以看到，对于使用基线获得的小分数类别，精确度、召回率和 f1 分数有了显著提高。1 混合本地股票 2的精确度为 0.52，现在精确度为 0.63。至于VSH 意大利蜜蜂，精确度为 0.68，现在为 0.97。我们注意到西方蜜蜂的精确度有所下降，但这个少数类的支持只有 7，所以这个结果是预期的。

我们继续优化我们的模型以提高验证和测试度量标准——换句话说，提高模型性能。在下一个迭代中，我们将训练周期数增加到 50。此外，我们将添加三个回调函数，如下所示：

• 一个学习率调度器，用于实现学习率变化的非线性函数。通过在每个周期改变学习率，我们可以改善训练过程。我们引入的用于控制学习率的函数将逐渐降低学习函数的值。

• 一个早期停止器，基于损失函数的演变（如果损失在一定数量的周期内没有改善）和一个耐心因子（在监控函数没有看到任何改善之后的周期数，我们停止训练）来停止训练周期。

• 每次获得最佳准确度时，都会有一个检查指针来保存表现最佳的模型。这将使我们能够使用的不是最后一次训练周期的模型参数，而是所有周期中表现最佳的模型。

三个回调函数的代码如下：

annealer3 = LearningRateScheduler(lambda x: 1e-3 * 0.995 ** (x+config['no_epochs_3']))
earlystopper3 = EarlyStopping(monitor='loss', patience=config['patience'], verbose=config['verbose'])
checkpointer3 = ModelCheckpoint('best_model_3.h5',
                                monitor='val_accuracy',
                                verbose=config['verbose'],
                                save_best_only=True,
                                save_weights_only=True) 

适配模型的代码也如下所示：

train_model3  = model3.fit_generator(image_generator.flow(X_train, y_train, batch_size=config['batch_size']),
                        epochs=config['no_epochs_3'],
                        validation_data=[X_val, y_val],
                        steps_per_epoch=len(X_train)/config['batch_size'],
                        callbacks=[earlystopper3, checkpointer3, annealer3]) 

训练可能需要分配的最大周期数，或者如果满足早期停止标准（即，在等于耐心因子的周期数之后没有损失函数的改善），它可能会提前结束。无论如何，实现了最佳验证准确度的模型将被保存并用于测试。

在图 6.27中，我们展示了此进一步优化的模型的训练和验证准确度以及训练和验证损失的变化。最终获得的训练损失为 0.18，最终训练准确度为 0.93。对于验证，最后的验证损失为 0.21，最后的验证准确度为 0.91。在最后一个周期，学习率为 6.08e-4。最佳验证准确度是在第 46 个周期获得的，为 0.92。

图 6.27：经过优化的模型（版本 3，包含学习率调度器、提前停止和检查点）——训练和验证准确率（左）和训练和验证损失（右）

我们使用保存的模型检查点（用于第 46 个 epoch）来预测测试数据。在图 6.28中，我们展示了第三个模型的分类报告。

宏平均指标进一步改进，分别达到精确度、召回率和 f1 分数的 0.88、0.90 和 0.89。精确度、召回率和 f1 分数的加权平均值也分别提高到 0.91、0.90 和 0.90。

图 6.28：使用第三个优化模型的测试数据分类报告（训练 epochs 增加到 50，并添加了学习率调度器、提前停止器和模型检查点）

在这里，我们将停止迭代改进模型的过程。您可以继续对其进行优化。您可以尝试添加更多的卷积和 maxpool 层，使用不同的内核数量和步长值来处理不同的超参数，包括不同的批量大小或学习率调度器。您还可以更改优化方案。改变模型的另一种方法是，通过数据增强控制类别图像的平衡（目前，蜜蜂图像在亚种类别上不平衡）。

您还可以尝试各种数据增强参数，并尝试使用不同的数据增强解决方案。参见参考文献 5，了解一个目前非常流行的图像数据增强库的示例，Albumentations，由一群数据科学家、研究人员和计算机视觉工程师创建，其中包括著名的 Kaggle 大师 Vladimir Iglovikov。

摘要

在本章中，我们首先介绍了一个新的数据集，其中包含了在不同日期和不同地点收集的图像元数据，包括各种患有不同疾病的蜜蜂亚种。我们还介绍了一些函数，用于基于skimage.io和 opencv (cv2)读取、缩放和从图像中提取特征。

我们使用了一个新创建的实用脚本，基于 Plotly 可视化表格数据，并利用 Plotly 的灵活性创建定制的图形，从而创建了有洞察力的可视化。我们还创建了图像的可视化函数。

在详细的数据探索分析（EDA）之后，我们转向构建蜜蜂亚种的预测模型。在这里，我们介绍了一种数据增强方法，通过从原始图像集中创建变化（旋转、缩放、平移和镜像）来增加初始可用训练数据。我们使用 stratify 将数据分为训练、验证和测试子集，以在随机采样三个子集时考虑类别不平衡。我们首先训练和验证了一个基线模型，然后，在执行错误分析后，我们逐步改进了初始模型，通过添加更多步骤，引入 Dropout 层，然后使用几个回调：学习率调度器、早期停止器和模型检查点。我们分析了训练、验证和测试错误的迭代改进，不仅关注训练和验证损失和准确率，还关注测试数据的分类报告。

在下一章中，我们将介绍文本数据分析的技术和工具，展示您如何准备数据以创建使用文本数据的基线模型。

参考文献

1. BeeImage 数据集：标注的蜜蜂图像：www.kaggle.com/datasets/jenny18/honey-bee-annotated-images

2. plotly-script 和 Kaggle 工具脚本：github.com/PacktPublishing/Developing-Kaggle-Notebooks/blob/develop/Chapter-06/plotly-utils.ipynb

3. 蜜蜂亚种分类，Kaggle 笔记本：github.com/PacktPublishing/Developing-Kaggle-Notebooks/blob/develop/Chapter-06/honeybee-subspecies-classification.ipynb

4. 安德鲁·吴，机器学习渴望：info.deeplearning.ai/machine-learning-yearning-book

5. Keras：keras.io/

6. TensorFlow：www.tensorflow.org/

7. 使用 Albumentations 与 Tensorflow：github.com/albumentations-team/albumentations_examples/blob/master/notebooks/tensorflow-example.ipynb

加入我们书籍的 Discord 空间

加入我们的 Discord 社区，与志同道合的人相聚，并在以下地点与超过 5000 名成员一起学习：

packt.link/kaggle

第七章：文本分析一切所需

在本章中，我们将学习如何分析文本数据并创建机器学习模型来帮助我们。我们将使用Jigsaw 无意中在毒性分类中的偏差数据集（参见参考文献 1）。这个竞赛的目标是构建检测毒性和减少可能错误地与毒性评论相关联的少数族裔的不当偏见的模型。通过这个竞赛，我们引入了自然语言处理（NLP）领域。

竞赛中使用的数据来源于由 Aja Bogdanoff 和 Christa Mrgan 于 2015 年创立的 Civil Comments 平台（参见参考文献 2），旨在解决在线讨论中的礼貌问题。当平台于 2017 年关闭时，他们选择保留约 200 万条评论，供研究人员理解并改善在线对话中的礼貌。Jigsaw 组织赞助了这项工作，并随后开始了一场语言毒性分类竞赛。在本章中，我们将把纯文本转换为有意义的、模型准备好的数字，以便根据评论的毒性将它们分类到不同的组别。

简而言之，本章将涵盖以下主题：

• 对Jigsaw 无意中在毒性分类中的偏差竞赛数据集的数据进行探索

• 介绍 NLP 特定的处理和分析技术，包括词频、分词、词性标注、命名实体识别和词嵌入

• 对文本数据的预处理进行迭代优化以准备模型基线

• 为这次文本分类竞赛创建一个模型基线

数据中有什么？

来自Jigsaw 无意中在毒性分类中的偏差竞赛数据集的数据包含训练集中的 180 万行和测试集中的 97,300 行。测试数据仅包含一个评论列，不包含目标（预测值）列。训练数据除了评论列外，还包括另外 43 列，包括目标特征。目标是介于 0 和 1 之间的数字，代表预测此竞赛的目标注释。此目标值表示评论的毒性程度（0表示零/无毒性，1表示最大毒性），其他 42 列是与评论中存在某些敏感主题相关的标志。主题与五个类别相关：种族和民族、性别、性取向、宗教和残疾。更详细地说，以下是每个类别的标志：

• 种族和民族：亚洲人、黑人、犹太人、拉丁美洲人、其他种族或民族和白人

• 性别：女性、男性、跨性别和其他性别

• 性取向：双性恋、异性恋、同性恋/女同性恋和其他性取向

• 宗教：无神论者、佛教徒、基督徒、印度教徒、穆斯林和其他宗教

• 残疾：智力或学习障碍、其他残疾、身体残疾和精神或心理疾病

还有一些特征（即数据集中的列）用于识别评论：created_data、publication_id、parent_id和article_id。还提供了与评论相关的一些用户反馈信息特征：rating、funny、wow、sad、likes、disagree和sexual_explicit。最后，还有两个与注释相关的字段：identity_annotator_count和toxicity_annotator_count。

让我们从对目标特征和敏感特征的快速分析开始。

目标特征

我们首先想看看目标特征的分布。让我们看看图 7.1中这些值的分布直方图：

图 7.1：目标值分布（训练数据，190 万列）

对于这个直方图，我们在Y轴上使用了对数刻度；这样做的原因是我们想看到值向 0 的偏斜分布。当我们这样做时，我们观察到我们有一个双峰分布：峰值在约 0.1 的间隔处，振幅减小，并且出现频率较低，趋势缓慢上升，叠加。大部分的目标值（超过 100 万）都是0。

敏感特征

我们将查看之前列出的敏感特征的分布（种族和民族、性别、性取向、宗教和残疾）。由于分布的偏斜（与目标类似，我们在0处有集中），我们将在Y轴上再次使用对数刻度。

图 7.2展示了种族和民族特征值的分布。这些值看起来不连续且非常离散，直方图显示了几个分离的峰值：

图 7.2：种族和民族特征值的分布

我们可以在图 7.3中观察到性别特征值有类似的分布：

图 7.3：性别特征值分布

在图 7.4中，我们展示了额外毒性特征值（severe_toxicity、obscene、identity_attack、insult或threat）的分布。如您所见，分布更加均匀，并且insult有增加的趋势：

图 7.4：额外毒性特征值的分布

让我们也看看目标值与种族或民族、性别、性取向、宗教和残疾特征值之间的相关性。我们在此不展示所有特征的相关矩阵，但您可以在与本章节相关的笔记本中检查它。在这里，我们只展示了与目标相关的第一个 16 个特征，按相关性因素排序：

train_corr['target'].sort_values(ascending=False)[1:16] 

让我们看看图 7.5中按相关性因素与目标排序的前 15 个特征：

图 7.5：其他特征与目标特征的前 15 个相关因素

接下来，在图 7.6中，我们展示了这些选定特征与目标特征之间的相关矩阵：

图 7.6：目标与相关性最高的 15 个特征之间的相关矩阵

我们可以观察到target与insult（0.93）、obscene（0.49）和identity_attack（0.45）高度相关。此外，severe_toxicity与insult和obscene呈正相关。identity_attack与white、black、muslim和homosexual_gay_or_lesbian有轻微的相关性。

我们研究了target特征（预测特征）和敏感特征的分布。现在我们将转向本章分析的主要主题：评论文本。我们将应用几种 NLP 特定的分析技术。

分析评论文本

NLP 是人工智能的一个领域，它涉及使用计算技术使计算机能够理解、解释、转换甚至生成人类语言。NLP 使用多种技术、算法和模型来处理和分析大量文本数据集。在这些技术中，我们可以提到：

• 分词：将文本分解成更小的单元，如单词、词的一部分或字符

• 词形还原或词干提取：将单词还原为词典形式或去除最后几个字符以得到一个共同形式（词干）

• 词性标注（POS）标记：将语法类别（例如，名词、动词、专有名词和形容词）分配给序列中的每个单词

• 命名实体识别（NER）：识别和分类实体（例如，人名、组织名和地名）

• 词嵌入：使用高维空间来表示单词，在这个空间中，每个单词的位置由其与其他单词的关系决定

• 机器学习模型：在标注数据集上训练模型以学习语言数据中的模式和关系

NLP 应用可以包括情感分析、机器翻译、问答、文本摘要和文本分类。

在对 NLP 进行快速介绍之后，让我们检查实际的评论文本。我们将构建几个词云图（使用整个数据集的 20,000 条评论子集）。我们首先查看整体单词分布（参见本章相关的笔记本），然后查看目标值高于 0.75 和低于 0.25 的单词分布：

图 7.7：低目标分数<0.25（左）和高目标分数>0.75（右）的评论中的常见单词（1-gram）

目标与侮辱高度相关，我们预计在两个特征的词云中会看到相当接近的单词分布。这一假设得到了证实，图 7.8很好地说明了这一点（无论是低分还是高分）：

图 7.8：低侮辱分数<0.25（左）和高侮辱分数>0.75（右）的评论中的常见单词（1-gram）

如您所见，分布显示了目标和侮辱的低分和高分中高频相似单词。

在相关的笔记本中还有更多关于威胁、淫秽和其他特征的词云。这些词云为我们提供了对最频繁单词的良好初始直觉。我们将在构建词汇表和检查词汇覆盖范围部分对整个语料库词汇中的单词频率进行更详细的分析。目前，我们可以观察到，我们进行的分析仅限于单个单词的频率，而没有捕捉到这些单词在整个语料库中的分组情况——换句话说，各种单词是如何一起使用的，以及基于此，如何识别语料库中的主要主题。这种旨在揭示整个语料库潜在语义结构的处理称为主题建模。在此方法中，对单词共现模式的分析使我们能够揭示文本中存在的潜在主题。

相关笔记本中主题建模实现灵感的来源是一系列关于使用潜在狄利克雷分配进行主题建模的文章和教程（见参考文献5-10）。

主题建模

我们首先通过预处理评论文本，使用gensim库来消除特殊字符、常用词、连接词（或停用词）以及长度小于 2 的单词：

def preprocess(text):
    result = []
    for token in gensim.utils.simple_preprocess(text):
        if token not in gensim.parsing.preprocessing.STOPWORDS and len(token) > 2:
            result.append(token)
    return result 

以下代码将定义的preprocess函数应用于所有评论：

%%time
preprocessed_comments = train_subsample['comment_text'].map(preprocess) 

然后，我们使用gensim/corpora中的dictionary创建一个单词字典。我们还过滤极端值，以消除不常见的单词并限制词汇表的大小：

%%time
dictionary = gensim.corpora.Dictionary(preprocessed_comments)
dictionary.filter_extremes(no_below=10, no_above=0.5, keep_n=75000) 

在这些限制下，我们接下来生成一个从字典中生成的词袋（bow）语料库。然后，我们对这个语料库应用TF-IDF（词频-逆文档频率），它为文档在文档集合或语料库中的重要性提供了一个数值表示。

tf组件衡量一个单词在文档中出现的频率。idf组件显示一个单词在整个文档语料库中的重要性（在我们的案例中，是整个评论集合）。这个因素随着一个术语在文档中出现的频率增加而降低。因此，在tfidf转换之后，一个单词和一个文档的系数对于在语料库级别不常见但在当前文档中频率较高的单词来说更大：

%%time
bow_corpus = [dictionary.doc2bow(doc) for doc in preprocessed_comments]
tfidf = models.TfidfModel(bow_corpus)
corpus_tfidf = tfidf[bow_corpus] 

然后，我们应用潜在狄利克雷分配（lda），这是一种基于该语料库中单词频率生成主题的主题模型，使用gensim的并行处理实现（LdaMulticore）：

%%time
lda_model = gensim.models.LdaMulticore(corpus_tfidf, num_topics=20,
                                    id2word=dictionary, passes=2, workers=2) 

让我们用 5 个单词来表示前 10 个主题：

topics = lda_model.print_topics(num_words=5)
for i, topic in enumerate(topics[:10]):
    print("Train topic {}: {}".format(i, topic)) 

主题单词以与主题相关的相对权重显示，如下所示：

图 7.9：前 10 个主题，每个主题选择 5 个（最相关）单词

一旦我们提取了主题，我们就可以遍历文档并识别当前文档（在我们的案例中，是评论）中包含哪些主题。在图 7.10中，我们展示了单个文档的占主导地位的主题（以及相对权重）（以下是用以生成选定评论的主题列表的代码）：

for index, score in sorted(lda_model[bd5], key=lambda tup: -1*tup[1]):
    print("\nScore: {}\t \nTopic: {}".format(score, lda_model.print_topic(index, 5))) 

图 7.10：与一个评论相关联（每个都有相对重要性）的主题

我们更喜欢使用pyLDAvis可视化工具来表示主题。在图 7.11中，我们展示了该工具的屏幕截图（在笔记本中，我们为训练数据生成了 20 个主题，并为测试数据单独生成了主题）。图 7.11中的仪表板显示了主题间距离图。在这里，主题的相对维度（或影响力）由圆盘的大小表示，主题的相对距离由它们的相互距离表示。在右侧，对于当前选定的主题（在左侧面板中），我们可以看到每个主题前 30 个最相关的术语。

色彩较浅的圆盘（笔记本中的蓝色）代表整体单词频率。颜色较深的圆盘（笔记本中的红色）代表所选主题内的估计单词频率。我们可以使用滑块来调整相关性指标（在图片中，这设置为1）。我们可以通过改进预处理步骤（例如，我们可以添加更多特定于该语料库的停用词）、调整字典形成的参数以及控制tfidf和lda的参数来进一步细化此分析。由于 LDA 过程的复杂性，我们还通过从训练数据中子采样来减小语料库的大小。

如以下屏幕截图所示，在 pyLDAvis 工具生成的主题建模仪表板的左侧面板中，我们看到主题间距离图——在语料库中主题影响力的相对维度和相对主题的距离：

图 7.11：使用 pyLDAvis 工具生成的主题建模仪表板（左侧面板）

在使用 pyLDAvis 工具生成的主题建模仪表板的右侧面板中，对于选定的主题，我们可以看到每个主题中最相关的 30 个术语，蓝色代表语料库中的整体词频，红色代表所选主题内的估计词频：

图 7.12：使用 pyLDAvis 工具生成的主题建模仪表板（右侧面板）

我们可以对整个语料库进行重复分析，但这将需要比 Kaggle 上可用的计算资源更多的资源。

通过以上内容，我们已经使用lda方法探索了评论文本语料库中的主题。通过这个程序，我们揭示了文本语料库中的一个隐藏（或潜在）结构。现在我们可以更好地理解不仅单词的频率，而且单词在评论中的关联方式，从而形成评论者讨论的主题。让我们继续从不同的角度探索语料库。我们将对每个评论进行分析，使用 NER 来确定文本中存在哪些类型的概念。然后，我们将开始关注句法元素，并使用 POS 标记来提取名词、动词、形容词和其他词性。

我们回顾这些 NLP 技术的原因有两个。首先，我们希望让你一窥 NLP 中可用的工具和技术之丰富。其次，对于更复杂的机器学习模型，你可以包括使用这些方法导出的特征。例如，除了从文本中提取的其他特征外，你还可以添加通过 NER 或 POS 标记获得的特征。

命名实体识别

让我们在一组评论上执行 NER。NER 是一种信息提取任务，旨在识别和提取非结构化数据（文本）中的命名实体。命名实体包括人名、组织、地理位置、日期和时间、数量和货币。有几种可用的方法来识别和提取命名实体，其中最常用的是spacy和transformers。在我们的案例中，我们将使用spacy来执行 NER。我们更喜欢spacy，因为它与 transformers 相比需要的资源更少，但仍然能给出良好的结果。在此需要注意的是，spacy也支持包括英语、葡萄牙语、西班牙语、俄语和中文在内的 23 种语言。

首先，我们使用spacy.load函数初始化一个nlp对象：

import spacy
nlp = spacy.load('en_core_web_sm') 

这将加载'en_core_web_sm'（sm代表小型）spacy管道，该管道包括tok2vec、tagger、parser、senter、ner、attribute_ruler和lemmatizer组件。我们不会使用此管道提供的所有功能；我们感兴趣的是nlp组件。

然后，我们创建了一个评论选择。我们筛选出包含名称 奥巴马 或 特朗普 且字符数少于 100 的文档。出于演示目的，我们不希望操作大型句子；如果我们使用较小的句子进行操作，将更容易跟随演示。下一个代码片段将执行选择：

selected_text = train.loc[train['comment_text'].str.contains("Trump") | train['comment_text'].str.contains("Obama")] 
selected_text["len"] = selected_text['comment_text'].apply(lambda x: len(x))
selected_text = selected_text.loc[selected_text.len < 100]
selected_text.shape 

我们可以通过两种方式可视化应用 NER 的结果。一种方式是打印出文本中每个识别实体的起始和结束字符以及实体标签。另一种方式是使用来自 spacy 的 displacy 渲染，它将为每个实体添加选定的颜色，并在实体文本旁边添加实体名称（见 图 7.13）。

以下代码用于使用 nlp 提取实体以及使用 displacy 准备可视化。在显示使用 displacy 注释的文本之前，我们打印出每个实体文本，然后是其位置（起始和结束字符位置）和实体标签：

for sentence in selected_text["comment_text"].head(5):
    print("\n")
    doc = nlp(sentence)
    for ent in doc.ents:
        print(ent.text, ent.start_char, ent.end_char, ent.label_)
    displacy.render(doc, style="ent",jupyter=True) 

spacy nlp 中预定义了多个标签。我们可以用一段简单的代码提取每个标签的含义：

import spacy
nlp = spacy.load("en_core_web_sm")
labels = nlp.get_pipe("ner").labels
for label in labels:
    print(f"{label} - {spacy.explain(label)}") 

这里是标签列表及其含义：

• CARDINAL: 不属于其他类型的数字

• DATE: 绝对或相对日期或时间段

• EVENT: 命名的飓风、战役、战争、体育赛事等

• FAC: 建筑物、机场、高速公路、桥梁等

• GPE: 国家、城市或州

• LANGUAGE: 任何命名语言

• LAW: 被制成法律的命名文件

• LOC: 非 GPE 地点、山脉或水体

• MONEY: 货币价值，包括单位

• NORP: 国籍或宗教或政治团体

• ORDINAL: 第一、第二等

• ORG: 公司、机构、机构等

• PERCENT: 百分比，包括 %

• PERSON: 人，包括虚构人物

• PRODUCT: 物体、车辆、食品等（不包括服务）

• QUANTITY: 重量或距离等度量

• TIME: 低于一天的时间

• WORK_OF_ART: 书籍、歌曲等的标题

图 7.13：使用 spacy 和 displacy 进行 NER 结果的可视化

在前面的截图顶部示例中，Bernie（Sanders）被正确识别为人物（PERSON），而（Donald）Trump 被识别为组织（ORG）。这可能是因为前总统 Trump 在作为商人时，经常将他的名字作为他创立的几个组织的名称的一部分。

在底部示例中，奥巴马（也是一位前总统，在争议性政治辩论中经常成为话题）被正确识别为人物。在两种情况下，我们还展示了提取的实体列表，包括每个识别实体的起始和结束位置。

词性标注

通过 NER 分析，我们识别了各种实体（如人、组织、地点等）的特定名称。这些帮助我们将各种术语与特定的语义组关联起来。我们可以进一步探索评论文本，以便了解每个单词的 POS（如名词或动词），并理解每个短语的句法。

让我们首先使用nltk（一个替代的nlp库）从我们在 NER 实验中使用的相同短语小集中提取词性。我们在这里选择nltk是因为，除了比 spacy 更节省资源外，它还提供了高质量的结果。我们还想能够比较两者的结果（spacy和nltk）：

for sentence in selected_text["comment_text"].head(5):
    print("\n")
    tokens = twt().tokenize(sentence)
    tags = nltk.pos_tag(tokens, tagset = "universal")
    for tag in tags:
        print(tag, end=" ") 

结果将如下所示：

图 7.14：使用 nltk 进行 POS 标记

我们也可以使用 spacy 执行相同的分析：

for sentence in selected_text["comment_text"].head(5):
    print("\n")
    doc = nlp(sentence)
    for token in doc:
        print(token.text, token.pos_, token.ent_type_, end=" | ") 

结果将如下所示：

图 7.15：使用 spacy 进行 POS 标记

让我们比较一下图 7.14和图 7.15中的两个输出。这两个库生成的 POS 结果略有不同。其中一些差异是由于实际词性在类别上的映射不同。对于nltk，单词“is”代表一个AUX（辅助），而spacy中的相同“is”则是一个动词。spacy区分专有名词（人名、地名等）和普通名词（NOUN），而nltk则不进行区分。

对于一些具有非标准结构的短语，两个输出都将动词“Go”错误地识别为名词（nltk）和专有名词（spacy）。在spacy的情况下，这是可以预料的，因为“Go”在逗号后面被大写。spacy 区分了并列连词（CCONJ）和从属连词（SCONJ），而nltk只会识别出存在连词（CONJ）。

使用与我们在前一小节中用于突出 NER 的相同 spacy 库扩展，我们也可以表示短语和段落的句法结构。在图 7.16中，我们展示了一个这样的表示示例。在笔记本中，我们使用带有“dep”（依存）标志的displacy展示了所有注释（短语集）的可视化。

图 7.16：使用 spacy 进行 POS 标记，并使用依存关系显示词性之间的短语结构

我们看到了如何使用依存关系来显示实体及其各自的类别，以及如何使用相同的函数来显示词性和短语结构。从参考文献 11中汲取灵感，我们扩展了那里给出的代码示例（并从使用nltk转换为使用spacy进行 POS 提取，因为 nltk 与 spacy 不完全兼容）以便我们可以以与表示命名实体相同的方式显示词性。

这里给出了从参考文献 11修改后的代码（代码解释将在代码块之后进行）：

import re
def visualize_pos(sentence):
    colors = {"PRON": "blueviolet",
              "VERB": "lightpink",
              "NOUN": "turquoise",
              "PROPN": "lightgreen",
              "ADJ" : "lime",
              "ADP" : "khaki",
              "ADV" : "orange",
              "AUX" : "gold",
              "CONJ" : "cornflowerblue",
              "CCONJ" : "magenta",
              "SCONJ" : "lightmagenta",
              "DET" : "forestgreen",
              "NUM" : "salmon",
              "PRT" : "yellow",
              "PUNCT": "lightgrey"}

    pos_tags = ["PRON", "VERB", "NOUN", "PROPN", "ADJ", "ADP",
                "ADV", "AUX", "CONJ", "CCONJ", "SCONJ",  "DET", "NUM", "PRT", "PUNCT"]

    # Fix for issues in the original code
    sentence = sentence.replace(".", " .")
    sentence = sentence.replace("'", "")
    # Replace nltk tokenizer with spacy tokenizer and POS tagging
    doc = nlp(sentence)
    tags = []
    for token in doc:
        tags.append((token.text, token.pos_))

    # Get start and end index (span) for each token
    span_generator = twt().span_tokenize(sentence)
    spans = [span for span in span_generator]
    # Create dictionary with start index, end index, 
# pos_tag for each token
    ents = []
    for tag, span in zip(tags, spans):
        if tag[1] in pos_tags:
            ents.append({"start" : span[0], 
                         "end" : span[1], 
                         "label" : tag[1] })
    doc = {"text" : sentence, "ents" : ents}
    options = {"ents" : pos_tags, "colors" : colors}

    displacy.render(doc, 
                    style = "ent", 
                    options = options, 
                    manual = True,
                   ) 

让我们更好地理解前面的代码。在visualise_pos函数中，我们首先定义了词性与颜色之间的映射（词性如何被突出显示）。然后，我们定义了我们将考虑的词性。接着，我们使用一些特殊字符的替换来纠正原始代码（来自参考文献 11）中存在的一个错误。我们还使用 spacy 分词器，并在tags列表中添加了使用nlp从spacy提取的每个pos的文本和词性。然后，我们计算每个pos的位置，并创建一个字典，包含pos标记和它们在文本中的位置，以便能够用不同的颜色突出显示它们。最后，我们使用displacy渲染所有突出显示的pos的文档。

在图 7.17中，我们展示了将此程序应用于我们的样本评论的结果。现在我们可以更容易地看到 spacy 的一些错误。在第二条评论中，它错误地将第二个“Go”解释为专有名词（PROPN）。这在某种程度上是可以解释的，因为通常在逗号之后，只有专有名词才会用大写字母书写。

图 7.17：使用 spacy 进行词性标注和依赖关系，并使用从参考文献 8 修改的程序来在文本中突出显示词性

我们还可以观察到其他错误。在第一条评论中，“Trump”被标记为名词——即一个简单的名词。术语“republicans”被分类为PROPN，这在美国政治背景下可能是准确的，因为在美国政治中，“Republicans”被视为专有名词。然而，在我们的上下文中，这是不准确的，因为它代表了一个复数形式的简单名词，指代一群倡导共和政府的人。

我们回顾了几种自然语言处理（NLP）技术，这些技术帮助我们更好地理解文本中存在的词语分布、主题、词性（POS）和概念。可选地，我们也可以使用这些技术来生成特征，以便包含在机器学习模型中。

在下一节中，我们将开始分析，目的是为评论分类准备一个监督式 NLP 模型。

准备模型

模型准备，根据我们将要实施的方法，可能更复杂或更简单。在我们的情况下，我们选择从简单的深度学习架构开始构建第一个基线模型（这是比赛时的标准方法），包括一个词嵌入层（使用预训练的词嵌入）和一个或多个双向长短期记忆（LSTM）层。这种架构在这次比赛发生时是一个常见的选择，对于文本分类问题的基线来说仍然是一个好的选择。LSTM代表长短期记忆。它是一种循环神经网络架构，旨在捕捉和记住序列数据中的长期依赖关系。它特别适用于文本分类问题，因为它能够处理和模拟文本序列中的复杂关系和依赖关系。

为了做到这一点，我们需要对评论数据进行一些预处理（在准备构建主题建模模型时我们也进行了预处理）。这次，我们将逐步执行预处理步骤，并监控这些步骤如何影响结果，不是模型的结果，而是表现良好的语言模型的一个先决条件，即词嵌入的词汇覆盖范围。

然后，我们将使用第一个基线模型中的词嵌入来扩展我们模型的泛化能力，以便不在训练集中但在测试集中的单词能够从存在于词嵌入中的单词的邻近性中受益。最后，为了确保我们的方法有效，我们需要预训练的词嵌入具有尽可能大的词汇覆盖范围。因此，我们还将测量词汇覆盖范围并建议改进它的方法。

现在，我们首先构建初始词汇表。

构建词汇表

我们之前使用单词频率、与目标值相关的单词分布以及其他特征、主题建模、命名实体识别（NER）和词性标注（POS tagging）对整个评论语料库的一个子集进行了实验。在接下来的实验中，我们将开始使用整个数据集。我们将使用具有 300 个嵌入大小的词嵌入。

词嵌入是单词的数值表示。它们将单词映射到向量。嵌入大小指的是这些向量的组成部分（或维度）的数量。这个过程使计算机能够理解和比较单词之间的关系。因为所有单词首先都使用词嵌入进行转换（在词嵌入空间中，单词之间的关系由向量之间的关系表示），所以具有相似意义的单词将在词嵌入空间中由对齐的向量表示。

在测试时，新词，不在训练数据中出现的词，也将被表示在词嵌入空间中，并且算法将利用它们与训练数据中存在的其他词的关系。这将增强我们用于文本分类的算法。

此外，我们将字符数（或注释长度）设置为固定值；我们选择了 220 这个维度。对于较短的注释，我们将填充注释序列（即添加空格），而对于较长的注释序列，我们将截断它们（到 220 个字符）。这个程序将确保我们会有相同维度的机器学习模型的输入。让我们首先定义一个用于构建词汇表的函数。在构建这些函数时，我们使用了来自参考文献 12和13的来源。

以下代码用于构建词汇表（即注释中存在的单词的语料库）。我们对每个注释进行分割，并将所有数据收集到一个句子列表中。然后我们解析所有这些句子以创建一个包含词汇的字典。每次解析到的单词如果在字典中作为键存在，我们就增加与该键关联的值。我们得到的是一个包含词汇表中每个单词计数（或总体频率）的词汇表字典：

def build_vocabulary(texts):
    """
    Build the vocabulary from the corpus
    Credits to: [9] [10]
    Args:
        texts: list of list of words
    Returns:
        dictionary of words and their count
    """
    sentences = texts.apply(lambda x: x.split()).values
    vocab = {}
    for sentence in tqdm(sentences):
        for word in sentence:
            try:
                vocab[word] += 1
except KeyError:
                vocab[word] = 1
return vocab 

我们通过连接train和test创建整体词汇表：

# populate the vocabulary
df = pd.concat([train ,test], sort=False)
vocabulary = build_vocabulary(df['comment_text']) 

我们可以检查词汇表中的前 10 个元素，以了解这个词汇表的样子：

# display the first 10 elements and their count
print({k: vocabulary[k] for k in list(vocabulary)[:10]}) 

以下图像显示了运行前面代码的结果。它显示了文本中最常用的单词。正如预期的那样，最常用的单词是英语中最常用的一些单词。

图 7.18：未经任何预处理的词汇表 – 大写和小写单词，以及可能拼写错误的表达

我们将在每次执行额外的（有时重复的）文本转换时重复使用之前引入的build_vocabulary函数。我们执行连续的文本转换，以确保在使用预训练的词嵌入时，我们对注释中词汇表的预训练词嵌入中的单词有良好的覆盖。更大的覆盖范围将确保我们构建的模型有更高的准确性。让我们继续加载一些预训练的词嵌入。

嵌入索引和嵌入矩阵

我们现在将构建一个字典，以词嵌入中的单词作为键，以它们的嵌入表示的数组作为值。我们称这个字典为嵌入索引。然后我们也将构建嵌入矩阵，它是嵌入的矩阵表示。我们将使用 GloVe 的预训练嵌入（300 维）进行我们的实验。GloVe代表全局词表示向量，是一种无监督算法，它产生词嵌入。它通过分析非常大的文本语料库中的全局文本统计来创建向量表示，并捕捉单词之间的语义关系。

以下代码加载预训练的词嵌入：

def load_embeddings(file):
    """
    Load the embeddings
    Credits to: [9] [10]
    Args:
        file: embeddings file
    Returns:
        embedding index
    """
def get_coefs(word,*arr): 
        return word, np.asarray(arr, dtype='float32')
    embeddings_index = dict(get_coefs(*o.split(" ")) for o in open(file, encoding='latin'))
    return embeddings_index
%%time
GLOVE_PATH = '../input/glove840b300dtxt/'
print("Extracting GloVe embedding started")
embed_glove = load_embeddings(os.path.join(GLOVE_PATH,'glove.840B.300d.txt'))
print("Embedding completed") 

获得的嵌入结构大小为 2.19 百万项。接下来，我们将使用我们刚刚创建的词索引和嵌入索引来创建嵌入矩阵：

def embedding_matrix(word_index, embeddings_index):
    '''
    Create the embedding matrix
    credits to: [9] [10] 
    Args:
        word_index: word index (from vocabulary)
        embedding_index: embedding index (from embeddings file)
    Returns:
        embedding matrix
    '''
    all_embs = np.stack(embeddings_index.values())
    emb_mean, emb_std = all_embs.mean(), all_embs.std()
    EMBED_SIZE = all_embs.shape[1]
    nb_words = min(MAX_FEATURES, len(word_index))
    embedding_matrix = np.random.normal(emb_mean, emb_std, (nb_words, EMBED_SIZE))
    for word, i in tqdm(word_index.items()):
        if i >= MAX_FEATURES:
            continue
        embedding_vector = embeddings_index.get(word)
        if embedding_vector is not None:
            embedding_matrix[i] = embedding_vector
    return embedding_matrix 

我们使用MAX_FEATURES参数来限制嵌入矩阵的维度。

检查词汇覆盖率

我们介绍了读取词嵌入和计算嵌入矩阵的功能。现在我们将继续介绍用于评估词嵌入中词汇覆盖率的函数。词汇覆盖率越大，我们构建的模型精度越好。

为了检查嵌入对词汇的覆盖率，我们将使用以下函数：

def check_coverage(vocab, embeddings_index):
    '''
    Check the vocabulary coverage by the embedding terms
    credits to: [9] [10]
    Args:
        vocab: vocabulary
        embedding_index: embedding index (from embeddings file)
    Returns:
        list of unknown words; also prints the vocabulary coverage of embeddings and 
        the % of comments text covered by the embeddings
    '''
    known_words = {}
    unknown_words = {}
    nb_known_words = 0
    nb_unknown_words = 0
for word in tqdm(vocab.keys()):
        try:
            known_words[word] = embeddings_index[word]
            nb_known_words += vocab[word]
        except:
            unknown_words[word] = vocab[word]
            nb_unknown_words += vocab[word]
            pass
print('Found embeddings for {:.3%} of vocabulary'.format(len(known_words)/len(vocab)))
    print('Found embeddings for {:.3%} of all text'.format(nb_known_words/(nb_known_words + nb_unknown_words)))
    unknown_words = sorted(unknown_words.items(), key=operator.itemgetter(1))[::-1]
    return unknown_words 

上述代码遍历所有词汇项（即评论文本中存在的单词）并计算未知单词（即文本中但不在嵌入单词列表中的单词）。然后，它计算词汇中存在于词嵌入索引中的单词百分比。这个百分比以两种方式计算：对词汇中的每个单词进行无权重计算，以及根据文本中的单词频率进行加权计算。

我们将反复应用此函数，在每个预处理步骤之后检查词汇覆盖率。让我们从检查初始词汇的词汇覆盖率开始，此时我们还没有对评论文本进行任何预处理。

逐步提高词汇覆盖率

我们应用check_coverage函数来检查词汇覆盖率，传递两个参数：词汇和嵌入矩阵。在以下符号中，oov代表词汇外：

print("Verify the initial vocabulary coverage")
oov_glove = check_coverage(vocabulary, embed_glove) 

第一次迭代的成果并不理想。尽管我们覆盖了几乎 90%的文本，但词汇表中只有 15.5%的单词被词嵌入覆盖：

图 7.19：词汇覆盖率 – 第一次迭代

我们还可以查看未覆盖术语的列表。因为在oov_glove中，我们按在语料库中出现的次数降序存储了未覆盖的术语，因此，通过选择此列表中的前几个术语，我们可以看到未包含在词嵌入中的最重要的单词。在图 7.20中，我们显示了此列表中的前 10 个术语——前 10 个未覆盖的单词。在这里，“未覆盖”指的是出现在词汇表中的单词（存在于评论文本中），但不在词嵌入索引中（不在预训练的词嵌入中）：

图 7.20：第一轮迭代中，词汇表中未覆盖的词汇中最常见的 10 个单词

通过快速检查图 7.20中的列表，我们看到一些常用词要么是缩写词，要么是非标准口语英语的口语化形式。在网上评论中看到这样的形式是正常的。我们将执行几个预处理步骤，试图通过纠正我们发现的问题来提高词汇覆盖率。在每一步之后，我们还将再次测量词汇覆盖率。

转换为小写

我们将首先将所有文本转换为小写，并将其添加到词汇表中。在词嵌入中，单词将全部为小写：

def add_lower(embedding, vocab):
    '''
    Add lower case words
    credits to: [9] [10]
    Args:
        embedding: embedding matrix
        vocab: vocabulary
    Returns:
        None
        modify the embeddings to include the lower case from vocabulary
    '''
    count = 0
for word in tqdm(vocab):
        if word in embedding and word.lower() not in embedding:  
            embedding[word.lower()] = embedding[word]
            count += 1
print(f"Added {count} words to embedding") 

我们将此小写转换应用于训练集和测试集，然后我们重新构建词汇表并计算词汇覆盖率：

train['comment_text'] = train['comment_text'].apply(lambda x: x.lower())
test['comment_text'] = test['comment_text'].apply(lambda x: x.lower())
print("Check coverage for vocabulary with lower case")
oov_glove = check_coverage(vocabulary, embed_glove)
add_lower(embed_glove, vocabulary) # operates on the same vocabulary
oov_glove = check_coverage(vocabulary, embed_glove) 

图 7.21显示了应用小写转换后新的词汇覆盖率：

图 7.21：词汇覆盖率——第二次迭代，所有单词的小写

我们可以在单词百分比和文本百分比覆盖率上观察到一些小的改进。让我们继续通过从注释文本中移除缩写词来继续操作。

移除缩写词

接下来，我们将移除缩写词。这些是单词和表达式的修改形式。我们将使用一个预定义的通常遇到的缩写词字典。这些将映射到存在于嵌入中的单词上。由于空间有限，我们在这里只包括缩写词字典中的一些项目示例，但整个资源可以在与本章相关的笔记本中找到：

contraction_mapping = {"ain't": "is not", "aren't": "are not","can't": "cannot", "'cause": "because", "could've": "could have", "couldn't": "could not", "didn't": "did not",  "doesn't": "does not", "don't": "do not", "hadn't": "had not", "hasn't": "has not",
...
} 

使用以下函数，我们可以获取 GloVe 嵌入中已知缩写词的列表：

def known_contractions(embed):
    '''
    Add know contractions
    credits to: [9] [10]
    Args:
        embed: embedding matrix
    Returns:
        known contractions (from embeddings)
    '''
    known = []
    for contract in tqdm(contraction_mapping):
        if contract in embed:
            known.append(contract)
    return known 

我们可以使用下一个函数从词汇表中清除已知的缩写词——即，使用缩写词字典来替换它们：

def clean_contractions(text, mapping):
    '''
    Clean the contractions

    credits to: [9] [10]
    Args:
        text: current text
        mapping: contraction mappings
    Returns: modify the comments to use the base form from contraction mapping
    '''
    specials = ["’", "‘", "´", "`"]
    for s in specials:
        text = text.replace(s, "'")
    text = ' '.join([mapping[t] if t in mapping else t for t in text.split(" ")])
    return text 

在我们对训练集和测试集应用clean_contractions并再次应用构建词汇表和测量词汇覆盖率的函数后，我们得到了关于词汇覆盖率的新统计数据：

图 7.22：词汇覆盖率——第三轮迭代，使用缩写词字典替换缩写词后的结果

通过检查没有覆盖的表达式并增强它，使其等于语料库中未覆盖的表达式与单词或单词组相等，可以进一步细化缩写字典。每个单词都表示在嵌入向量中。

移除标点和特殊字符

接下来，我们将移除标点和特殊字符。以下列表和函数对此步骤很有用。首先，我们列出未知标点：

punct_mapping = "/-'?!.,#$%\'()*+-/:;<=>@[\\]^_`{|}~" + '""“”’' + '∞θ÷α•à−βø³π'₹´°£€\×™√²—–&'
punct_mapping += '©^®` <→°€™' ♥←ק″′Â█½à…"✶"–●â►−¢²¬░¶↑±¿▾═¦║―¥▓—'─▒: ¼⊕▼▪†■’▀¨▄♫⭐é¯♦¤▲踾Ñ'∞∙）↓、│（»，♪╩╚³・╦╣╔╗▬❤ïع≤‡√'
def unknown_punct(embed, punct):
    '''
    Find the unknown punctuation
    credits to: [9] [10] 
    Args:
        embed: embedding matrix
        punct: punctuation
    Returns:
        unknown punctuation
    '''
    unknown = ''
for p in punct:
        if p not in embed:
            unknown += p
            unknown += ' '
return unknown 

然后我们清理特殊字符和标点：

puncts = {"‘": "'", "´": "'", "°": "", "€": "euro", "—": "-", "–": "-", "’": "'", "_": "-", "`": "'", '“': '"', '”': '"', '“': '"', "£": "pound",
          '∞': 'infinity', 'θ': 'theta', '÷': '/', 'α': 'alpha', '•': '.', 'à': 'a', '−': '-', 'β': 'beta', '∅': '', '³': '3', 'π': 'pi', '…': ' '}
def clean_special_chars(text, punct, mapping):
    '''
    Clean special characters
    credits to: [9] [10]
    Args:
        text: current text
        punct: punctuation
        mapping: punctuation mapping
    Returns:
        cleaned text
    '''
for p in mapping:
        text = text.replace(p, mapping[p])
    for p in punct:
        text = text.replace(p, f' {p} ') 
    return text 

让我们在图 7.23中再次检查词汇覆盖率。这次，通过词嵌入，我们将词汇覆盖率从大约 15%提高到 54%。此外，文本覆盖率从 90%提高到 99.7%。

图 7.23：词汇覆盖率 – 第四次迭代，在清理标点和特殊字符之后

查看未覆盖的前 20 个单词，我们发现有一些带重音的小词、特殊字符和习语。我们将标点字典扩展到包括最频繁的特殊字符，并在我们再次运行build_vocabulary和check_coverage之后，我们得到了词汇覆盖率的新状态：

more_puncts = {'▀': '.', '▄': '.', 'é': 'e', 'è': 'e', 'ï': 'i','⭐': 'star', 'ᴀ': 'A',  'ᴀɴᴅ': 'and', '»': ' '}
train['comment_text'] = train['comment_text'].apply(lambda x: clean_special_chars(x, punct_mapping, more_puncts))
test['comment_text'] = test['comment_text'].apply(lambda x: clean_special_chars(x, punct_mapping, more_puncts))
%%time
df = pd.concat([train ,test], sort=False)
vocab = build_vocabulary(df['comment_text'])
print("Check coverage after additional punctuation replacement")
oov_glove = check_coverage(vocab, embed_glove) 

这次有一个微小的改进，但我们可以继续处理频繁的表达式或频繁的特殊字符替换，直到我们得到显著的改进。

通过添加额外的嵌入源到当前的预训练嵌入中，可以进一步改善注释语料库的词汇覆盖率。让我们试试这个。我们使用了来自GloVe的预训练嵌入。我们也可以使用 Facebook 的FastText。FastText是一个非常实用的行业标准库，在许多公司的日常搜索和推荐引擎中广泛使用。让我们加载嵌入并使用组合嵌入向量重新创建嵌入索引。

在我们将两个词嵌入字典合并后，一个包含 2.19 百万个条目，另一个包含 2.0 百万个条目（两者都具有 300 维的向量维度），我们得到了一个包含 2.8 百万个条目的字典（由于两个字典中有很多共同词汇）。然后我们重新计算了词汇覆盖率。在图 7.24中，我们展示了这一操作的结果。

图 7.24：词汇覆盖率 – 第五次迭代，在将 FastText 预训练词嵌入添加到初始 GloVe 嵌入字典之后

总结一下我们的过程，我们的意图是构建一个基于预训练词嵌入的基线解决方案。我们介绍了两种预训练词嵌入算法，GloVe 和 FastText。预训练意味着我们使用了已经训练好的算法；我们没有从我们的数据集的评论语料库中计算词嵌入。为了有效，我们需要确保我们用这些词嵌入覆盖了评论文本词汇表的良好范围。最初，覆盖率相当低（词汇表的 15% 和整个文本的 86%）。通过转换为小写、删除缩写、删除标点符号和替换特殊字符，我们逐渐提高了这些统计数据。在最后一步，我们通过添加来自替代源的预训练嵌入来扩展嵌入字典。最终，我们能够确保词汇表的 56% 覆盖率和整个文本的 99.75% 覆盖率。

下一步是创建一个基线模型，并在单独的笔记本中完成。我们只会重用当前笔记本中实验的一部分函数。

构建基线模型

这些天，每个人至少都会通过微调 Transformer 架构来构建一个基线模型。自从 2017 年的论文《Attention Is All You Need》（参考文献 14）以来，这些解决方案的性能一直在持续提升，对于像 Jigsaw Unintended Bias in Toxicity Classification 这样的比赛，一个基于 Transformer 的最新解决方案可能会轻易地将你带入金牌区域。

在这个练习中，我们将从一个更经典的基线开始。这个解决方案的核心是基于 Christof Henkel（Kaggle 昵称：Dieter）、Ane Berasategi（Kaggle 昵称：Ane）、Andrew Lukyanenko（Kaggle 昵称：Artgor）、Thousandvoices（Kaggle 昵称）和 Tanrei（Kaggle 昵称）的贡献；参见参考文献 12、13、15、16、17 和 18。

该解决方案包括四个步骤。在第一步，我们将训练数据和测试数据加载为 pandas 数据集，然后对这两个数据集进行预处理。预处理主要基于我们之前执行的预处理步骤，因此我们在这里不会重复这些步骤。

在第二步，我们执行分词并准备数据以呈现给模型。分词过程如下所示（我们这里不展示整个过程）：

 logger.info('Fitting tokenizer')
    tokenizer = Tokenizer() 
    tokenizer.fit_on_texts(list(train[COMMENT_TEXT_COL]) + list(test[COMMENT_TEXT_COL]))
    word_index = tokenizer.word_index
    X_train = tokenizer.texts_to_sequences(list(train[COMMENT_TEXT_COL]))
    X_test = tokenizer.texts_to_sequences(list(test[COMMENT_TEXT_COL]))
    X_train = pad_sequences(X_train, maxlen=MAX_LEN)
    X_test = pad_sequences(X_test, maxlen=MAX_LEN) 

我们在这里使用了一个基本的分词器来自 keras.preprocessing.text。分词后，每个输入序列都会用预定义的 MAX_LEN 进行填充，这个长度是根据整个评论语料库的平均/中位长度以及可用的内存和运行时限制选定的。

在第三步，我们构建嵌入矩阵和模型结构。构建嵌入矩阵的代码主要基于我们在前几节中已经介绍过的过程。在这里，我们只是将其系统化：

def build_embedding_matrix(word_index, path):
    '''
     Build embeddings
    '''
    logger.info('Build embedding matrix')
    embedding_index = load_embeddings(path)
    embedding_matrix = np.zeros((len(word_index) + 1, EMB_MAX_FEAT))
    for word, i in word_index.items():
        try:
            embedding_matrix[i] = embedding_index[word]
        except KeyError:
            pass
except:
            embedding_matrix[i] = embeddings_index["unknown"]

    del embedding_index
    gc.collect()
    return embedding_matrix
def build_embeddings(word_index):
    '''
     Build embeddings
    '''
    logger.info('Load and build embeddings')
    embedding_matrix = np.concatenate(
        [build_embedding_matrix(word_index, f) for f in EMB_PATHS], axis=-1) 
    return embedding_matrix 

该模型是一个具有词嵌入层、SpatialDropout1D层、两个双向 LSTM 层、GlobalMaxPooling1D与GlobalAveragePooling1D的拼接、两个具有'relu'激活的密集层，以及一个具有'sigmoid'激活的目标输出的密集层的深度学习架构。

在词嵌入层中，输入被转换，使得每个词都由其对应的向量表示。经过这种转换后，模型捕捉到了输入中词语之间的语义距离信息。SpatialDropout1D层通过在训练过程中随机停用神经元来帮助防止过拟合（系数给出了每个 epoch 停用神经元的百分比）。双向 LSTM 层的作用是处理输入序列的前向和反向，增强上下文理解以获得更好的预测。GlobalAveragePooling1D层的作用是计算整个序列中每个特征的均值，在 1D（顺序）数据中降低维度同时保留关键信息。这相当于揭示了序列的潜在表示。密集层的输出是模型的预测。有关实现的更多细节，请参阅参考文献 17和参考文献 18：

def build_model(embedding_matrix, num_aux_targets, loss_weight):
    '''
        Build model
    '''
    logger.info('Build model')
    words = Input(shape=(MAX_LEN,))
    x = Embedding(*embedding_matrix.shape, weights=[embedding_matrix], trainable=False)(words)
    x = SpatialDropout1D(0.3)(x)
    x = Bidirectional(CuDNNLSTM(LSTM_UNITS, return_sequences=True))(x)
    x = Bidirectional(CuDNNLSTM(LSTM_UNITS, return_sequences=True))(x)
    hidden = concatenate([GlobalMaxPooling1D()(x),GlobalAveragePooling1D()(x),])
    hidden = add([hidden, Dense(DENSE_HIDDEN_UNITS, activation='relu')(hidden)])
    hidden = add([hidden, Dense(DENSE_HIDDEN_UNITS, activation='relu')(hidden)])
    result = Dense(1, activation='sigmoid')(hidden)
    aux_result = Dense(num_aux_targets, activation='sigmoid')(hidden)

    model = Model(inputs=words, outputs=[result, aux_result])
    model.compile(loss=[custom_loss,'binary_crossentropy'], loss_weights=[loss_weight, 1.0], optimizer='adam')
    return model 

在第四步，我们进行训练，准备提交，并提交。为了减少运行时使用的内存，我们使用临时存储，并在删除未使用的分配数据后进行垃圾回收。我们运行模型两次，指定数量的NUM_EPOCHS（代表训练数据通过算法的一个完整遍历），然后使用可变权重平均测试预测。然后我们提交预测：

def run_model(X_train, y_train, y_aux_train, embedding_matrix, word_index, loss_weight):
    '''
        Run model
    '''
    logger.info('Run model')

    checkpoint_predictions = []
    weights = []
    for model_idx in range(NUM_MODELS):
        model = build_model(embedding_matrix, y_aux_train.shape[-1], loss_weight)
        for global_epoch in range(NUM_EPOCHS):
            model.fit(
                X_train, [y_train, y_aux_train],
                batch_size=BATCH_SIZE, epochs=1, verbose=1,
                callbacks=[LearningRateScheduler(lambda epoch: 1.1e-3 * (0.55 ** global_epoch))]
            )
            with open('temporary.pickle', mode='rb') as f:
                X_test = pickle.load(f) # use temporary file to reduce memory
            checkpoint_predictions.append(model.predict(X_test, batch_size=1024)[0].flatten())
            del X_test
            gc.collect()
            weights.append(2 ** global_epoch)
        del model
        gc.collect()

    preds = np.average(checkpoint_predictions, weights=weights, axis=0)
    return preds
def submit(sub_preds):
    logger.info('Prepare submission')
    submission = pd.read_csv(os.path.join(JIGSAW_PATH,'sample_submission.csv'), index_col='id')
    submission['prediction'] = sub_preds
    submission.reset_index(drop=False, inplace=True)
    submission.to_csv('submission.csv', index=False) 

使用这种解决方案（完整代码请见参考文献 16），我们可以通过延迟提交获得核心得分 0.9328，从而在私有排行榜的上半部分获得排名。接下来，我们将展示如何通过使用基于 Transformer 的解决方案，我们可以获得更高的分数，进入银牌甚至金牌区域。

基于 Transformer 的解决方案

在竞赛期间，BERT 和一些其他 Transformer 模型已经可用，并提供了几个高分解决方案。在这里，我们不会尝试复制它们，但我们会指出最易于实现的实现。

在参考文献 20中，齐申·哈结合了几种解决方案，包括 BERT-Small V2、BERT-Large V2、XLNet 和 GPT-2（使用竞赛数据作为数据集进行微调的模型）来获得 0.94656 的私有排行榜得分（延迟提交），这将使你进入前 10 名（包括金牌和奖项区域）。

仅使用 BERT-Small 模型（参见参考文献 21）的解决方案将产生 0.94295 的私有排行榜分数。使用 BERT-Large 模型（参见参考文献 22）将导致 0.94388 的私有排行榜分数。这两个解决方案都将位于银牌区域（在私有排行榜中分别位于 130 和 80 左右，作为晚提交）。

摘要

在本章中，我们学习了如何处理文本数据，使用各种方法来探索这类数据。我们首先分析了我们的目标和文本数据，对文本数据进行预处理以便将其包含在机器学习模型中。我们还探讨了各种 NLP 工具和技术，包括主题建模、命名实体识别（NER）和词性标注（POS tagging），然后准备文本以构建基线模型，通过迭代过程逐步提高目标集（在这种情况下，目标是提高竞赛数据集中文本语料库中词汇的词嵌入覆盖范围）的数据质量。

我们介绍了并讨论了一个基线模型（基于几位 Kaggle 贡献者的工作）。这个基线模型架构包括一个词嵌入层和双向 LSTM 层。最后，我们查看了一些基于 Transformer 架构的最先进解决方案，无论是作为单一模型还是组合使用，以获得排行榜上部的分数（银牌和金牌区域）。

在下一章中，我们将开始处理信号数据。我们将介绍各种信号模态（声音、图像、视频、实验或传感器数据）特定的数据格式。我们将分析来自LANL 地震预测 Kaggle 竞赛的数据。

参考文献

1. Jigsaw 毒性分类中的无意偏差，Kaggle 竞赛数据集：www.kaggle.com/c/jigsaw-unintended-bias-in-toxicity-classification/

2. Aja Bogdanoff，告别文明评论，Medium：medium.com/@aja_15265/saying-goodbye-to-civil-comments-41859d3a2b1d

3. Gabriel Preda，Jigsaw 评论文本探索：github.com/PacktPublishing/Developing-Kaggle-Notebooks/blob/develop/Chapter-07/jigsaw-comments-text-exploration.ipynb

4. Gabriel Preda，Jigsaw 简单基线：github.com/PacktPublishing/Developing-Kaggle-Notebooks/blob/develop/Chapter-07/jigsaw-simple-baseline.ipynb

5. Susan Li, Python 中的主题建模和潜在狄利克雷分配（LDA）：towardsdatascience.com/topic-modeling-and-latent-dirichlet-allocation-in-python-9bf156893c24

6. Aneesha Bakharia, 提高主题模型的解释能力：towardsdatascience.com/improving-the-interpretation-of-topic-models-87fd2ee3847d

7. Carson Sievert, Kenneth Shirley, LDAvis：一种可视化和解释主题的方法：www.aclweb.org/anthology/W14-3110

8. Lucia Dosin, 主题建模实验 – PyLDAvis：www.objectorientedsubject.net/2018/08/experiments-on-topic-modeling-pyldavis/

9. Renato Aranha, 在 Elon 推文上的主题建模（LDA）：www.kaggle.com/errearanhas/topic-modelling-lda-on-elon-tweets

10. 潜在狄利克雷分配，维基百科：en.wikipedia.org/wiki/Latent_Dirichlet_allocation

11. Leonie Monigatti, 使用 NLTK 和 SpaCy 可视化词性标签：towardsdatascience.com/visualizing-part-of-speech-tags-with-nltk-and-spacy-42056fcd777e

12. Ane, Quora 预处理 + 模型：www.kaggle.com/anebzt/quora-preprocessing-model

13. Christof Henkel (Dieter), 如何：在使用嵌入时进行预处理：www.kaggle.com/christofhenkel/how-to-preprocessing-when-using-embeddings

14. Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N. Gomez, Lukasz Kaiser, Illia Polosukhin, 注意力即是全部：arxiv.org/abs/1706.03762

15. Christof Henkel (Dieter), keras 基线 lstm + attention 5 折：www.kaggle.com/christofhenkel/keras-baseline-lstm-attention-5-fold

16. Andrew Lukyanenko, 在 keras 上使用 CNN 的折叠：www.kaggle.com/code/artgor/cnn-in-keras-on-folds

17. Thousandvoices, 简单 LSTM：www.kaggle.com/code/thousandvoices/simple-lstm/s

18. Tanrei, 使用身份参数的简单 LSTM 解决方案：www.kaggle.com/code/tanreinama/simple-lstm-using-identity-parameters-solution

19. Gabriel Preda, Jigsaw Simple Baseline: www.kaggle.com/code/gpreda/jigsaw-simple-baseline

20. Qishen Ha, Jigsaw_predict: www.kaggle.com/code/haqishen/jigsaw-predict/

21. Gabriel Preda, Jigsaw_predict_BERT_small: www.kaggle.com/code/gpreda/jigsaw-predict-bert-small

22. Gabriel Preda, Jigsaw_predict_BERT_large: www.kaggle.com/code/gpreda/jigsaw-predict-bert-large

加入我们书籍的 Discord 空间

加入我们的 Discord 社区，与志同道合的人相聚，并在以下地点与超过 5000 名成员一起学习：

packt.link/kaggle

第八章：分析声学信号以预测下一次模拟地震

在前面的章节中，我们探讨了基本的表格格式数据，包括分类、有序和数值数据，以及文本、地理坐标和图像。当前章节将我们的关注点转向不同的数据类别，特别是模拟或实验信号数据。这种数据类型通常以多种格式出现，而不仅仅是标准的 CSV 文件格式。

我们的主要案例研究将是LANL 地震预测 Kaggle 竞赛的数据（参见参考文献 1）。我为此竞赛贡献了一个广为人知且经常被分叉的笔记本，名为LANL 地震 EDA 和预测（参见参考文献 2），它将作为本章主要笔记本的基础资源。然后我们将深入研究特征工程，采用各种对开发竞赛预测模型至关重要的信号分析技术。我们的目标是构建一个初始模型，该模型可以预测竞赛的目标变量：故障时间，即下一次模拟实验室地震前的剩余时间。

在地震预测领域的的研究表明，在地震发生之前，板块的运动会在低频声学频谱中产生信号。通过研究这些信号，研究人员试图理解信号的轮廓与故障（即地震）发生时刻之间的关系。在实验室中，模拟板块的滑动和剪切。这个竞赛使用实验室测量数据，包括声学信号以及故障发生的时间。

总结来说，本章将涵盖以下主题：

• 用于各种信号数据的数据格式

• 探索LANL 地震预测 Kaggle 竞赛数据

• 特征工程

• 训练LANL 地震预测竞赛的模型

介绍 LANL 地震预测竞赛

LANL 地震预测竞赛的核心是利用地震信号来确定实验室诱导地震的确切时间。目前，预测自然地震仍然超出了我们的科学知识和技术能力。科学家理想的情景是预测此类事件的时间、地点和震级。

然而，在高度控制的模拟环境中创造的模拟地震，模仿了真实的地震活动。这些模拟允许使用在自然环境中观察到的相同类型的信号来预测实验室生成的地震。在这场比赛中，参与者使用声学数据输入信号来估计下一次人工地震发生的时间，如参考文献 3中详细说明。挑战在于预测地震的时间，解决地震预测中的三个关键未知因素之一：它将在何时发生，将在哪里发生，以及它将有多强大。

训练数据是一个包含两列的单个文件：声学信号幅度和时间至失效。测试数据由多个文件组成（总共 2,526 个），包含声学信号幅度段，我们将需要预测每个段的时间至失效。一个样本提交文件包含一个列，即段 ID seg_id 和要预测的值：time_to_failure。

竞赛者需要使用训练文件中的声学信号和时间至失效数据来训练他们的模型，并预测测试文件夹中每个文件的每个段的时间至失效。这些竞赛数据格式非常方便，即逗号分隔值（CSV）格式，但这不是必需的。Kaggle 上其他带有信号数据的竞赛或数据集使用不同的、不太常见的格式。因为本章是关于分析信号数据的，所以这里是回顾这种格式的正确位置。让我们首先了解一下这些格式中的一些。

信号数据的格式

在 Kaggle 上举办的几场比赛使用了声音数据作为常规表格特征的补充。2021 年、2022 年和 2023 年，康奈尔鸟类实验室的 BirdCLEF（LifeCLEF 鸟类识别挑战）组织了三场比赛，用于从鸟鸣样本中预测鸟种（参见参考文献 4以了解这些比赛中的一个示例）。这些比赛使用的格式是.ogg。.ogg格式用于以更少的带宽存储音频数据，技术上被认为优于.mp3格式。

我们可以使用librosa库（参见参考文献 5）读取这些类型的文件格式。以下代码可以用来加载一个.ogg文件并显示声音波形：

import matplotlib.pyplot as plt
import librosa
def display_sound_wave(sound_path=None,
               text="Test", 
               color="green"):
    """
    Display a sound wave
    Args
        sound_path: path to the sound file
        text: text to display
        color: color for text to display
    Returns
        None
    """
if not sound_path:
        return    
    y_sound, sr_sound = librosa.load(sound_path)
    audio_sound, _ = librosa.effects.trim(y_sound)
    fig, ax = plt.subplots(1, figsize = (16, 3))
    fig.suptitle(f'Sound Wave: {text}', fontsize=12)
    librosa.display.waveshow(y = audio_sound, sr = sr_sound, color = color) 

当librosa库加载音频声音时，将返回一个时间序列，包含浮点数值（参见参考文献 6）。它不仅支持.ogg格式；它还可以与 soundfile 或 Audioread 支持的任何代码一起工作。默认采样率为 22050，但也可以在加载时设置，使用参数sr。在加载音频波形时还可以使用的其他参数是偏移量和持续时间（两者都以秒为单位 - 一起，它们允许您选择要加载的声音波形的时问间隔）。

在 BirdCLEF 竞赛的早期版本中，康奈尔鸟鸣识别（见参考文献 7），数据集中的音频声音以.mp3格式给出。对于这种格式，我们可以使用 librosa 来加载、转换或可视化声音波。波形音频文件格式（或WAV），另一种常用格式，也可以使用 librosa 加载。

对于.wav格式，我们可以使用scipy.io模块的wavfile来加载数据。以下代码将加载并显示一个.wav格式的文件。在这种情况下，振幅没有缩放到-1:1 的区间（最大值为 32K）：

import matplotlib.pyplot as plt
from scipy.io import wavfile
def display_wavefile(sound_path=None,
               text="Test", 
               color="green"):
    """
    Display a sound wave - load using wavefile
    sr: sample rate
    y_sound: sound samples
    Args
        sound_path: path to the sound file
        text: text to display
        color: color for text to display
    Returns
        None
    """
if not sound_path:
        return    
    sr_sound, y_sound = wavfile.load(sound_path)
    fig, ax = plt.subplots(1, figsize = (16, 3))
    fig.suptitle(f'Sound Wave: {text}', fontsize=12)
    ax.plot(np.linspace(0, sr_sound/len(y_sound), sr_sound), y_sound) 

信号，不仅仅是音频信号，也可以存储在.npy或.npz格式中，这两种格式都是numpy格式，用于存储数组数据。这些格式可以使用numpy函数加载，如下面的代码片段所示。对于.npy格式，这将加载一个多列数组：

import numpy as np
f = np.load('data_path/file.npy', allow_pickle=True)
columns_, data_  = f
data_df = pd.DataFrame(data_, columns = columns_) 

对于.npz格式，以下代码将加载一个类似的结构，之前已压缩（只有一个文件）：

import numpy as np
f = np.load('data_path/file.npz', allow_pickle=True)
columns_, data_ = f['arr_0']
data_df = pd.DataFrame(data_, columns = columns_) 

对于存储在.rds格式中的数据，这是一种用于保存数据的 R 特定格式，我们可以使用以下代码加载数据：

!pip install pyreadr
import pyreadr
f = pyreadr.read_r('data_path/file.rds')
data_df = f[None] 
for CO, focusing on the COCL dimension (Column Burden kg m-2), and includes values for latitude, longitude, and time:

from netCDF4 import Dataset
data = Dataset(file_path, more="r")
lons = data.variables['lon'][:]
lats = data.variables['lat'][:]
time = data.variables['time'][:]
COCL = data.variables['COCL'][:,:,:]; COCL = COCL[0,:,:] 

更多详情，请参阅参考文献 9。现在，让我们回到我们的竞赛数据，它以 CSV 格式存储，尽管它代表音频信号（声波），正如我们之前所阐明的。

探索我们的竞赛数据

LANL 地震预测数据集包含以下数据：

• 一个只有两列的train.csv文件：

  • acoustic_data: 这是声学信号的振幅。

  • time_to_failure: 这是指当前数据段对应的故障时间。

• 一个包含 2,624 个文件的小段声学数据的测试文件夹。

• 一个sample_submission.csv文件；对于每个测试文件，参赛者需要提供一个故障时间的估计。

训练数据（9.56 GB）包含 6.92 亿行。训练数据中样本的实际时间常数是由time_to_failure值的连续变化引起的。声学数据是整数，从-5,515 到 5,444，平均值为 4.52，标准差为 10.7（值在 0 附近波动）。time_to_failure值是实数，范围从 0 到 16，平均值为 5.68，标准差为 3.67。为了减少训练数据的内存占用，我们以减少的维度读取声学数据和time_to_failure：

%%time
train_df = pd.read_csv(os.path.join(PATH,'train.csv'), dtype={'acoustic_data': np.int16, 'time_to_failure': np.float32}) 

让我们检查training数据中的第一个值。我们不会使用所有的time_to_failure数据（只有与时间间隔结束相关的值，我们将对这些值进行聚合以获取声学数据）；因此，将时间到故障的大小从双精度浮点数缩小到浮点数不是很重要：

图 8.1. 训练数据中的第一行数据

让我们在同一张图上可视化声音信号值和故障时间。我们将使用 1/100 的子采样率（每 100 行取一个样本）来表示完整训练数据（参见图 8.2）。我们将使用以下代码来表示这些图表：

def plot_acc_ttf_data(idx, train_ad_sample_df, train_ttf_sample_df, title="Acoustic data and time to failure: 1% sampled data"):
    """
    Plot acoustic data and time to failure
    Args:
        train_ad_sample_df: train acoustic data sample
        train_ttf_sample_df: train time to failure data sample
        title: title of the plot
    Returns:
        None
    """
    fig, ax1 = plt.subplots(figsize=(12, 8))
    plt.title(title)
    plt.plot(idx, train_ad_sample_df, color='r')
    ax1.set_ylabel('acoustic data', color='r')
    plt.legend(['acoustic data'], loc=(0.01, 0.95))
    ax2 = ax1.twinx()
    plt.plot(idx, train_ttf_sample_df, color='b')
    ax2.set_ylabel('time to failure', color='b')
    plt.legend(['time to failure'], loc=(0.01, 0.9))
    plt.grid(True) 

图 8.2. 整个训练集的声音信号数据和故障时间数据，以 1/100 的比例进行子采样

让我们放大时间间隔的第一部分。我们将展示前 1%的数据（不进行子采样）。在图 8.3中，我们在同一张图上展示了前 6.29 百万行数据的声音信号和故障时间。我们可以观察到在故障前（但时间上不是非常接近），有一个大的振荡，既有负峰也有正峰。这个振荡还由几个较小的不规则振荡在非规则的时间间隔前导。

图 8.3：数据前 1%的声音信号数据和故障时间数据

让我们也看看训练数据的下一个 1%（不进行子采样）。在图 8.4中，我们展示了声音信号值和故障时间的时序图。在这个时间间隔内没有发生故障。我们观察到许多不规则的小振荡，既有负峰也有正峰。

图 8.4：训练集中第 2%的数据的声音信号数据和故障时间

让我们也看看训练集中数据的最后几个百分比（最后 5%的时间）在训练集中。在图 8.5中，我们观察到几个较大的振荡叠加在较小的不规则振荡上，并且在故障前有一个主要振荡：

图 8.5：训练集中最后 5%的数据的声音信号数据和故障时间

现在，让我们也看看测试数据样本中声音信号变化的一些例子。测试数据中有 2,624 个数据段文件。我们将从中选择一些进行可视化。由于测试数据中我们只有声音信号，我们将使用修改后的可视化函数：

def plot_acc_data(test_sample_df, segment_name):
    """
    Plot acoustic data for a train segment
    Args:
        test_sample_df: test acoustic data sample
        segment_name: title of the plot
    Returns:
        None
    """
    fig, ax1 = plt.subplots(figsize=(12, 8))
    plt.title(f"Test segment: {segment_name}")
    plt.plot(test_sample_df, color='r')
    ax1.set_ylabel('acoustic data', color='r')
    plt.legend([f"acoustic data: {segment_name}"], loc=(0.01, 0.95))
    plt.grid(True) 

在图 8.6中，我们展示了seg_00030f段的声音信号图：

图 8.6：测试段 seg_00030f 的声音信号数据

在下一张图中，我们展示了seg_0012b5段的声音信号图：

图 8.7：测试段 seg_0012b5 的声音信号数据

在与本章相关的笔记本中，您可以查看更多此类测试声学信号的示例。测试段显示了相当大的信号波形多样性，描述了相同的小振荡序列，其中穿插着不同振幅的峰值，类似于我们在之前下采样的训练数据中可以看到的。

解决方案方法

竞赛中的任务是准确预测测试数据集中每个段落的单一time_to_failure值。test集的每个段由 150,000 个数据行组成。相比之下，training数据集非常庞大，包含 6.92 亿行，其中一列专门用于我们的目标变量：失效时间。我们计划将训练数据分成均匀的段，每个段包含 150,000 行，并使用每个段的最终失效时间值作为该段的目标变量。这种方法旨在使训练数据与测试数据的格式相匹配，从而促进更有效的模型训练。

此外，我们还将通过聚合训练和测试数据集中的值来构建新的特征，从而得到一个包含每个数据段多个特征的单一行。下一节将深入探讨用于特征生成的信号处理技术。

特征工程

我们将使用几个特定于信号处理的库来生成大多数特征。从 SciPy（Python 科学库）中，我们使用了signal模块的一些函数。Hann 函数返回一个 Hann 窗口，它修改信号以平滑采样信号末尾的值到 0（使用余弦“钟形”函数）。Hilbert 函数通过Hilbert变换计算解析信号。Hilbert 变换是信号处理中使用的数学技术，具有将原始信号的相位移动 90 度的特性。

其他使用的库函数来自numpy：快速傅里叶变换（FFT）、mean、min、max、std（标准差）、abs（绝对值）、diff（信号中两个连续值之间的差值）和quantile（将样本分为相等大小、相邻的组）。我们还使用了一些来自pandas的统计函数：mad（中值绝对偏差）、kurtosis、skew和median。我们正在实现计算趋势特征和经典 STA/LTA 的函数。经典 STA/LTA 表示 STA 长度短时间窗口信号振幅与长时间窗口 LTA 的比率。让我们深入探讨吧！

趋势特征和经典 STA/LTA

我们首先定义两个函数，用于计算趋势特征和经典的短期平均/长期平均（STA/LTA）。STA/LTA 是地震信号分析技术，在地震学中应用。它测量短期信号平均值与长期信号平均值的比率。在地震检测中很有用，因为它可以识别地震数据中的独特模式。因此，它也将是我们模型中一个有用的特征。

我们在此处展示了计算趋势特征的代码。这是使用线性回归模型（用于 1D 数据）来检索结果回归线的斜率。我们在进行回归（即计算数据的绝对值的斜率/趋势）之前将所有采样数据转换为正值。趋势数据包含有关整体信号的重要信息：

def add_trend_feature(arr, abs_values=False):
    """
    Calculate trend features
    Uses a linear regression algorithm to extract the trend
    from the list of values in the array (arr)
    Args:
        arr: array of values
        abs_values: flag if to use abs values, default is False

    Returns:
        trend feature
    """
    idx = np.array(range(len(arr)))
    if abs_values:
        arr = np.abs(arr)
    lr = LinearRegression()
    lr.fit(idx.reshape(-1, 1), arr)
    return lr.coef_[0] 

接下来，我们计算经典的 STA/LTA，它表示长度为STA的短时间窗口信号幅度与长时间窗口LTA的比值。该函数接收信号和短时间平均窗口以及长时间平均窗口的长度作为参数：

def classic_sta_lta(x, length_sta, length_lta):
    """
    Calculate classic STA/LTA
    STA/LTA represents the ratio between amplitude of the
    signal on a short time window of length LTA and on a
    long time window LTA
    Args:
        length_sta: length of short time average window
        length_lta: length of long time average window
    Returns:
        STA/LTA
    """
    sta = np.cumsum(x ** 2)
    # Convert to float
    sta = np.require(sta, dtype=np.float)
    # Copy for LTA
    lta = sta.copy()
    # Compute the STA and the LTA
    sta[length_sta:] = sta[length_sta:] - sta[:-length_sta]
    sta /= length_sta
    lta[length_lta:] = lta[length_lta:] - lta[:-length_lta]
    lta /= length_lta
    # Pad zeros
    sta[:length_lta - 1] = 0
# Avoid division by zero by setting zero values to tiny float
    dtiny = np.finfo(0.0).tiny
    idx = lta < dtiny
    lta[idx] = dtiny
    return sta / lta 

接下来，我们实现一个计算特征的功能，该功能接收样本索引、数据子样本和转换后的训练数据句柄作为参数。此功能将使用各种信号处理算法从每个段的时间变化声学信号中构建聚合特征。在训练数据的情况下，我们使用训练集的 150K 行窗口（没有步长）。在测试集的情况下，每个测试文件代表 150K 的段。在以下小节中，我们将回顾将包含在模型中的工程特征。

由 FFT 派生的特征

模型的特征之一是对整个段应用快速傅里叶变换（FFT）；这并不是直接用作特征，而是作为计算多个聚合函数的基础（参见下一小节）。FFT 使用快速实现的离散傅里叶变换来计算。

我们使用numpy实现的一维数组FFT（fft.fft），这非常快，因为numpy基于BLAS（基本线性代数子程序）和Lapack（线性代数包），这两个库提供了执行基本向量和矩阵操作以及求解线性代数方程的程序。这里使用的函数的输出是一个复数值的一维数组。然后，我们从复数值数组中提取实部和虚部的向量，并计算以下特征：

• 提取 FFT 的实部和虚部；这是进一步处理声学信号快速傅里叶变换的第一步。

• 计算 FFT 的实部和虚部的平均值、标准差、最小值和最大值。从之前的变换中，该变换将 FFT 的实部和虚部分开，然后我们计算这些聚合函数。

• 计算 FFT 向量末尾 5K 和 15K 数据点的 FFT 的实部和虚部的平均值、标准差、最小值和最大值。

创建文件段以及 FFT 和由 FFT 派生的特征的代码如下。首先，我们计算声学数据子集的 FFT。然后，我们计算 FFT 的实部和虚部。从实 FFT 分量中，我们使用 pandas 的聚合函数计算平均值、标准差、最大值和最小值。然后，我们从 FFT 信号的虚部计算类似值：

def create_features(seg_id, seg, X):
    """
    Create features
    Args:
        seg_id: the id of current data segment to process
        seg: the current selected segment data
        X: transformed train data
    Returns:
        None
    """
    xc = pd.Series(seg['acoustic_data'].values)
    zc = np.fft.fft(xc)

    #FFT transform values
    realFFT = np.real(zc)
    imagFFT = np.imag(zc)
    X.loc[seg_id, 'Rmean'] = realFFT.mean()
    X.loc[seg_id, 'Rstd']  = realFFT.std()
    X.loc[seg_id, 'Rmax']  = realFFT.max()
    X.loc[seg_id, 'Rmin']  = realFFT.min()
    X.loc[seg_id, 'Imean'] = imagFFT.mean()
    X.loc[seg_id, 'Istd']  = imagFFT.std()
    X.loc[seg_id, 'Imax']  = imagFFT.max()
    X.loc[seg_id, 'Imin']  = imagFFT.min()
    X.loc[seg_id, 'Rmean_last_5000']  = realFFT[-5000:].mean()
    X.loc[seg_id, 'Rstd__last_5000']  = realFFT[-5000:].std()
    X.loc[seg_id, 'Rmax_last_5000']   = realFFT[-5000:].max()
    X.loc[seg_id, 'Rmin_last_5000']   = realFFT[-5000:].min()
    X.loc[seg_id, 'Rmean_last_15000'] = realFFT[-15000:].mean()
    X.loc[seg_id, 'Rstd_last_15000'] = realFFT[-15000:].std()
    X.loc[seg_id, 'Rmax_last_15000'] = realFFT[-15000:].max()
    X.loc[seg_id, 'Rmin_last_15000'] = realFFT[-15000:].min() 

我们接着计算由各种聚合函数派生的特征。

由聚合函数派生的特征

使用 pandas 的聚合函数mean、std、max和min计算整个段落的平均值、标准差、最大值和最小值的代码如下：

 xc = pd.Series(seg['acoustic_data'].values)
    zc = np.fft.fft(xc)

    X.loc[seg_id, 'mean'] = xc.mean()
    X.loc[seg_id, 'std'] = xc.std()
    X.loc[seg_id, 'max'] = xc.max()
    X.loc[seg_id, 'min'] = xc.min() 

我们继续计算额外的聚合特征。对于我们的模型，我们将包括各种信号处理技术，正如您将注意到的，然后，在训练我们的基线模型后，通过测量特征重要性，我们将确定哪些特征对我们的模型预测贡献更大。

接下来，我们计算整个段落的平均变化量；这里的“段落”指的是原始的声学数据子集。“变化”是通过numpy函数diff和参数1计算的。此函数接收一个值数组，并计算数组中每个连续值之间的差异。然后我们计算差异值数组的平均值。我们还计算整个声学数据段的变化率的平均值。这是通过将新变化向量中的非零值的平均值除以数据段中的原始值来计算的。这些特征的代码如下：

 X.loc[seg_id, 'mean_change_abs'] = np.mean(np.diff(xc))
    X.loc[seg_id, 'mean_change_rate'] = np.mean(nonzero(np.diff(xc) / xc[:-1])[0]) 

此外，我们还计算每个整个段落的绝对值（最大值和最小值）。在计算绝对值之后，我们再计算最小值和最大值。

当我们聚合时间信号时，我们试图包括一个多样化的特征范围，以尽可能多地捕捉信号模式。这个代码如下：

 X.loc[seg_id, 'abs_max'] = np.abs(xc).max()
    X.loc[seg_id, 'abs_min'] = np.abs(xc).min() 

可以计算每个声学数据段前 10K、50K 和最后 10K 值的聚合函数集，如下所示：

• 每个声学数据段前 50K 和最后 10K 值的平均值

• 每个声学数据段前 50K 和最后 10K 值的平均值

• 每个声学数据段前 50K 和最后 10K 值的最小值

• 每个声学数据段前 50K 和最后 10K 值的最大值

这些特征正在聚合信号的一部分较小部分，因此它们将只从故障前更小的时间间隔中捕获信号特征。在整个信号长度和信号较小部分上的聚合特征组合将增加更多关于信号的信息。这些特征的代码如下：

 X.loc[seg_id, 'std_first_50000'] = xc[:50000].std()
    X.loc[seg_id, 'std_last_50000'] = xc[-50000:].std()
    X.loc[seg_id, 'std_first_10000'] = xc[:10000].std()
    X.loc[seg_id, 'std_last_10000'] = xc[-10000:].std()

    X.loc[seg_id, 'avg_first_50000'] = xc[:50000].mean()
    X.loc[seg_id, 'avg_last_50000'] = xc[-50000:].mean()
    X.loc[seg_id, 'avg_first_10000'] = xc[:10000].mean()
    X.loc[seg_id, 'avg_last_10000'] = xc[-10000:].mean()

    X.loc[seg_id, 'min_first_50000'] = xc[:50000].min()
    X.loc[seg_id, 'min_last_50000'] = xc[-50000:].min()
    X.loc[seg_id, 'min_first_10000'] = xc[:10000].min()
    X.loc[seg_id, 'min_last_10000'] = xc[-10000:].min()

    X.loc[seg_id, 'max_first_50000'] = xc[:50000].max()
    X.loc[seg_id, 'max_last_50000'] = xc[-50000:].max()
    X.loc[seg_id, 'max_first_10000'] = xc[:10000].max()
    X.loc[seg_id, 'max_last_10000'] = xc[-10000:].max() 

接下来，我们包括整个声学数据段的极大值与极小值之比以及极大值与极小值之间的差值。我们还添加了超过一定振荡幅度（超过 500 个单位）的值的数量以及整个段的值总和。我们试图使用我们构建的这些特征多样性来捕获信号中的某些隐藏模式。特别是，这里我们包括信号极端振荡的信息：

 X.loc[seg_id, 'max_to_min'] = xc.max() / np.abs(xc.min())
    X.loc[seg_id, 'max_to_min_diff'] = xc.max() - np.abs(xc.min())
    X.loc[seg_id, 'count_big'] = len(xc[np.abs(xc) > 500])
    X.loc[seg_id, 'sum'] = xc.sum() 

我们继续添加多样化的聚合特征，试图捕获原始信号的各个特征。我们进一步计算每个声学数据段前 10K 和最后 50K 数据点的平均变化率（排除空值）：

 X.loc[seg_id, 'mean_change_rate_first_50000'] = np.mean(nonzero((np.diff(xc[:50000]) / xc[:50000][:-1]))[0])
    X.loc[seg_id, 'mean_change_rate_last_50000'] = np.mean(nonzero((np.diff(xc[-50000:]) / xc[-50000:][:-1]))[0])
    X.loc[seg_id, 'mean_change_rate_first_10000'] = np.mean(nonzero((np.diff(xc[:10000]) / xc[:10000][:-1]))[0])
    X.loc[seg_id, 'mean_change_rate_last_10000'] = np.mean(nonzero((np.diff(xc[-10000:]) / xc[-10000:][:-1]))[0]) 

我们添加的一些特征将排除数据中的0元素，以确保只将非零值包含在聚合函数的计算中。使用nonzero函数的代码如下：

def nonzero(x):
    """
    Utility function to simplify call of numpy `nonzero` function
    """
return np.nonzero(np.atleast_1d(x)) 

一组工程特征涉及分位数，具体是整个声学数据段的 01%，05%，95%，和 99%分位数值。分位数是使用numpy的quantile函数计算的。分位数是一个统计术语，指的是将数据集分成等概率的区间。例如，75%的分位数值是 75%的数据值小于该数值的点。50%的分位数是 50%的数据值小于该数值的点（也称为中位数）。我们还添加了 01%，05%，95%，和 99%分位数的绝对值。以下代码用于计算这些特征：

 X.loc[seg_id, 'q95'] = np.quantile(xc, 0.95)
    X.loc[seg_id, 'q99'] = np.quantile(xc, 0.99)
    X.loc[seg_id, 'q05'] = np.quantile(xc, 0.05)
    X.loc[seg_id, 'q01'] = np.quantile(xc, 0.01)

    X.loc[seg_id, 'abs_q95'] = np.quantile(np.abs(xc), 0.95)
    X.loc[seg_id, 'abs_q99'] = np.quantile(np.abs(xc), 0.99)
    X.loc[seg_id, 'abs_q05'] = np.quantile(np.abs(xc), 0.05)
    X.loc[seg_id, 'abs_q01'] = np.quantile (np.abs(xc), 0.01) 

另一种引入的工程特征类型是趋势值（使用不带绝对标志的add_trend_values函数计算）。趋势值将捕获声学数据信号变化的一般方向。对于显示围绕 0 高频振荡的信号，趋势将捕获实际信号平均值的变化。

我们还添加了绝对趋势值（使用带有绝对标志的add_trend_values函数计算）。我们包括这种工程特征来捕获信号绝对值中出现的模式。在这种情况下，对于趋势的计算，我们使用原始信号的绝对值。因此，这种趋势将捕获信号绝对值变化的方向。相应的代码如下：

 X.loc[seg_id, 'trend'] = add_trend_feature(xc)
    X.loc[seg_id, 'abs_trend'] = add_trend_feature(xc, abs_values=True) 

接下来，我们包括绝对值的平均值和绝对值的标准差。中位数绝对偏差（mad）、峰度、偏度（偏度）和中位数值也被计算。这些函数使用numpy实现。中位数绝对偏差是定量数据单变量样本变异性的稳健度量。峰度是分布尾部的综合权重相对于分布中心的度量。偏度（来自偏度）是对称分布的不对称或扭曲的度量。中位数，正如我们之前观察到的，是分隔数据集上半部分和下半部分值的数值。所有这些聚合函数都捕捉到关于信号的有益信息。这些聚合函数的计算代码如下所示：

 X.loc[seg_id, 'abs_mean'] = np.abs(xc).mean()
    X.loc[seg_id, 'abs_std'] = np.abs(xc).std()

    X.loc[seg_id, 'mad'] = xc.mad()
    X.loc[seg_id, 'kurt'] = xc.kurtosis()
    X.loc[seg_id, 'skew'] = xc.skew()
    X.loc[seg_id, 'med'] = xc.median() 

接下来，我们包括几个通过使用信号处理特定的变换函数计算的特征。

使用希尔伯特变换和汉宁窗口派生的特征

我们还计算希尔伯特均值。我们使用scipy.signal.hilbert函数对声学信号段应用希尔伯特变换。这通过希尔伯特变换计算解析信号，然后，我们计算变换数据的绝对值的平均值。希尔伯特变换在信号处理中经常被使用，并捕捉到关于信号的重要信息。因为我们使用聚合函数从我们的时间数据中生成特征，我们希望包括大量、多样化的现有信号处理技术，在训练模型时添加信号的重要补充元素：

 X.loc[seg_id, 'Hilbert_mean'] = np.abs(hilbert(xc)).mean() 

接下来，我们包括一个由汉宁窗口均值派生的特征。我们使用这个由汉宁窗口派生的特征来减少信号边缘的突然不连续性。汉宁窗口均值是通过将原始信号与汉宁窗口的结果进行卷积，然后除以汉宁窗口中所有值的总和来计算的：

 X.loc[seg_id, 'Hann_window_mean'] = (convolve(xc, hann(150), mode='same') / sum(hann(150))).mean() 

我们之前介绍了经典 STA/LTA 的定义。我们计算了 500-10K、5K-100K、3,333-6,666 和 10K-25K STA/LTA 窗口的经典 STA/LTA 均值等几个特征。这些是通过之前介绍的 STA/LTA 函数计算的。我们在聚合工程特征中包括各种变换，试图捕捉不同的信号特性：

 X.loc[seg_id, 'Hilbert_mean'] = np.abs(hilbert(xc)).mean()
    X.loc[seg_id, 'Hann_window_mean'] = (convolve(xc, hann(150), mode='same') / sum(hann(150))).mean()
    X.loc[seg_id, 'classic_sta_lta1_mean'] = classic_sta_lta(xc, 500, 10000).mean()
    X.loc[seg_id, 'classic_sta_lta2_mean'] = classic_sta_lta(xc, 5000, 100000).mean()
    X.loc[seg_id, 'classic_sta_lta3_mean'] = classic_sta_lta(xc, 3333, 6666).mean()
    X.loc[seg_id, 'classic_sta_lta4_mean'] = classic_sta_lta(xc, 10000, 25000).mean() 

最后，我们还将计算基于移动平均的特征。

基于移动平均的特征

接下来，我们计算几个移动平均值，如下所示：

• 700、1.5K、3K 和 6K 窗口的移动平均值均值（排除 NaNs）

• 指数加权移动平均，跨度为 300、3K 和 6K

• 700 和 400 窗口的平均标准差移动平均值

• 700 大小窗口加减 2 倍平均标准差的移动平均值

移动平均线帮助我们辨别模式，减少噪声，并清晰地展示数据中潜在趋势的图像。相应的代码如下：

 X.loc[seg_id, 'Moving_average_700_mean'] = xc.rolling(window=700).mean().mean(skipna=True)
    X.loc[seg_id, 'Moving_average_1500_mean'] = xc.rolling(window=1500).mean().mean(skipna=True)
    X.loc[seg_id, 'Moving_average_3000_mean'] = xc.rolling(window=3000).mean().mean(skipna=True)
    X.loc[seg_id, 'Moving_average_6000_mean'] = xc.rolling(window=6000).mean().mean(skipna=True)
    ewma = pd.Series.ewm
    X.loc[seg_id, 'exp_Moving_average_300_mean'] = (ewma(xc, span=300).mean()).mean(skipna=True)
    X.loc[seg_id, 'exp_Moving_average_3000_mean'] = ewma(xc, span=3000).mean().mean(skipna=True)
    X.loc[seg_id, 'exp_Moving_average_30000_mean'] = ewma(xc, span=6000).mean().mean(skipna=True)
    no_of_std = 2
    X.loc[seg_id, 'MA_700MA_std_mean'] = xc.rolling(window=700).std().mean()
    X.loc[seg_id,'MA_700MA_BB_high_mean'] = (X.loc[seg_id, 'Moving_average_700_mean'] + no_of_std * X.loc[seg_id, 'MA_700MA_std_mean']).mean()
    X.loc[seg_id,'MA_700MA_BB_low_mean'] = (X.loc[seg_id, 'Moving_average_700_mean'] - no_of_std * X.loc[seg_id, 'MA_700MA_std_mean']).mean()
    X.loc[seg_id, 'MA_400MA_std_mean'] = xc.rolling(window=400).std().mean()
    X.loc[seg_id,'MA_400MA_BB_high_mean'] = (X.loc[seg_id, 'Moving_average_700_mean'] + no_of_std * X.loc[seg_id, 'MA_400MA_std_mean']).mean()
    X.loc[seg_id,'MA_400MA_BB_low_mean'] = (X.loc[seg_id, 'Moving_average_700_mean'] - no_of_std * X.loc[seg_id, 'MA_400MA_std_mean']).mean()
    X.loc[seg_id, 'MA_1000MA_std_mean'] = xc.rolling(window=1000).std().mean() 

我们还计算四分位距（IQR）、0.01%和 99.9%的分位数。四分位距（IQR，即四分位数范围）是通过从 75%分位数减去 25%分位数（使用numpy函数）来计算的。四分位距是数据中 50%所在的范围。0.01%和 99.9%的分位数也使用numpy函数进行计算。IQR 和我们所包含的其他分位数是有用的，因为它们提供了关于信号集中趋势和分散的重要信息：

 X.loc[seg_id, 'iqr'] = np.subtract(*np.percentile(xc, [75, 25]))
    X.loc[seg_id, 'q999'] = np.quantile(xc,0.999)
    X.loc[seg_id, 'q001'] = np.quantile(xc,0.001)
    X.loc[seg_id, 'ave10'] = stats.trim_mean(xc, 0.1) 

对于 10、100 和 1,000 个窗口，我们计算移动平均和移动标准差。有了这些值，我们随后计算最小值、最大值、平均值、标准差、平均绝对变化和相对变化、1%、5%、95%和 99%的分位数，以及绝对最大滚动值。我们包括这些特征，因为它们揭示了在指定窗口内信号局部特性的信息。随后，对于从 10、100 和 1,000 个窗口计算得到的移动平均标准差派生出的特征，代码如下：

for windows in [10, 100, 1000]:
        x_roll_std = xc.rolling(windows).std().dropna().values

        X.loc[seg_id, 'ave_roll_std_' + str(windows)] = x_roll_std.mean()
        X.loc[seg_id, 'std_roll_std_' + str(windows)] = x_roll_std.std()
        X.loc[seg_id, 'max_roll_std_' + str(windows)] = x_roll_std.max()
        X.loc[seg_id, 'min_roll_std_' + str(windows)] = x_roll_std.min()
        X.loc[seg_id, 'q01_roll_std_' + str(windows)] = np.quantile(x_roll_std, 0.01)
        X.loc[seg_id, 'q05_roll_std_' + str(windows)] = np.quantile(x_roll_std, 0.05)
        X.loc[seg_id, 'q95_roll_std_' + str(windows)] = np.quantile(x_roll_std, 0.95)
        X.loc[seg_id, 'q99_roll_std_' + str(windows)] = np.quantile(x_roll_std, 0.99)
        X.loc[seg_id, 'av_change_abs_roll_std_' + str(windows)] = np.mean(np.diff(x_roll_std))
        X.loc[seg_id, 'av_change_rate_roll_std_' + str(windows)] = np.mean(nonzero((np.diff(x_roll_std) / x_roll_std[:-1]))[0])
        X.loc[seg_id, 'abs_max_roll_std_' + str(windows)] = np.abs(x_roll_std).max() 

对于从 10、100 和 1,000 个窗口计算得到的移动平均平均值派生出的特征，代码如下：

 for windows in [10, 100, 1000]:
        x_roll_mean = xc.rolling(windows).mean().dropna().values

        X.loc[seg_id, 'ave_roll_mean_' + str(windows)] = x_roll_mean.mean()
        X.loc[seg_id, 'std_roll_mean_' + str(windows)] = x_roll_mean.std()
        X.loc[seg_id, 'max_roll_mean_' + str(windows)] = x_roll_mean.max()
        X.loc[seg_id, 'min_roll_mean_' + str(windows)] = x_roll_mean.min()
        X.loc[seg_id, 'q01_roll_mean_' + str(windows)] = np.quantile(x_roll_mean, 0.01)
        X.loc[seg_id, 'q05_roll_mean_' + str(windows)] = np.quantile(x_roll_mean, 0.05)
        X.loc[seg_id, 'q95_roll_mean_' + str(windows)] = np.quantile(x_roll_mean, 0.95)
        X.loc[seg_id, 'q99_roll_mean_' + str(windows)] = np.quantile(x_roll_mean, 0.99)
        X.loc[seg_id, 'av_change_abs_roll_mean_' + str(windows)] = np.mean(np.diff(x_roll_mean))
        X.loc[seg_id, 'av_change_rate_roll_mean_' + str(windows)] = np.mean(nonzero((np.diff(x_roll_mean) / x_roll_mean[:-1]))[0])
        X.loc[seg_id, 'abs_max_roll_mean_' + str(windows)] = np.abs(x_roll_mean).max() 

对于从训练数据生成的每个 150K 行段，我们计算这些特征。然后，选择当前段最后一行的值作为故障时间：

# iterate over all segments
for seg_id in tqdm_notebook(range(segments)):
    seg = train_df.iloc[seg_id*rows:seg_id*rows+rows]
    create_features(seg_id, seg, train_X)
    train_y.loc[seg_id, 'time_to_failure'] = seg['time_to_failure'].values[-1] 

接下来，我们使用StandardScaler对所有特征进行缩放。如果我们使用基于决策树（如随机森林或 XGBoost）的模型，这一步不是必须的。我们包括这一步是为了考虑到我们可能想使用其他模型的情况，例如基于神经网络的模型，在这种情况下，对特征进行归一化将是一个必要的步骤：

scaler = StandardScaler()
scaler.fit(train_X)
scaled_train_X = pd.DataFrame(scaler.transform(train_X), columns=train_X.columns) 

我们对测试数据段重复同样的过程：

for seg_id in tqdm_notebook(test_X.index):
    seg = pd.read_csv('../input/LANL-Earthquake-Prediction/test/' + seg_id + '.csv')
    create_features(seg_id, seg, test_X)
scaled_test_X = pd.DataFrame(scaler.transform(test_X), columns=test_X.columns) 

在我们分析数据后，我们生成了一组工程特征。我们打算使用这些特征来构建一个基线模型。然后，基于模型评估，我们可以进一步选择保留哪些特征，并最终创建新的特征。

构建基线模型

从原始的时间数据中，通过特征工程，我们为训练数据中的每个时间段生成了时间聚合特征，其持续时间与一个测试集相等。对于本次比赛中展示的基线模型，我们选择了LGBMRegressor，这是当时表现最好的算法之一，在很多情况下，其性能与XGBoost相似。训练数据使用KFold分割成五个部分，我们对每个部分进行训练和验证，直到达到最终的迭代次数或者验证误差在指定步骤数后不再改善（由patience参数给出）。对于每个分割，我们还会对测试集进行预测，使用的是当前折叠的最佳模型——即使用当前折叠的训练分割进行训练的模型，也就是使用训练集的 4/5。最后，我们将计算每个折叠获得的预测的平均值。我们可以使用这种交叉验证方法，因为我们的数据不再是时间序列数据。我们将数据分割成 150K 行的段，与测试数据长度相同，然后从这些数据中创建了聚合特征：

n_fold = 5
folds = KFold(n_splits=n_fold, shuffle=True, random_state=42)
train_columns = scaled_train_X.columns.values 

模型参数如下代码片段所示。我们为 LightGBM 设置了以下一些常用参数：

• 叶子节点数：此参数控制每个树中的叶子节点（或终端节点）的数量。增加叶子节点的数量允许模型捕捉更复杂的模式，但也增加了过拟合的风险。

• 叶子节点中的最小数据量：如果一个叶子节点的样本数低于此阈值，则该节点不会分裂。此参数有助于控制过拟合。

• 目标：这是我们的模型的回归。

• 学习率：这将控制模型学习的速度。

• 提升方法：我们可以选择梯度提升决策树（gbdt）、dart 梯度提升（dgb）和基于梯度的单侧采样（goss）。在这里我们使用gbdt。

• 特征分数：这是在算法使用的树集成中，向树展示的数据子集中的特征百分比。

• 袋装频率、袋装分数和袋装种子：这些控制我们在对算法进行子采样以展示给不同的树时如何划分样本集。

• 度量标准：在这种情况下，mae，即平均绝对误差。

• 正则化因子 lambda_l1

• 详细程度：此参数控制算法在训练过程中打印到控制台的信息量。详细程度为 0 表示静默模式（无信息），而详细程度为 1 会打印有关训练进度的消息。

• 并行处理线程数：此参数控制算法在训练期间使用的并行线程数。

• 随机化因子 random_state：此参数是用于初始化算法各种参数的随机种子。

params = {'num_leaves': 51,
         'min_data_in_leaf': 10, 
         'objective':'regression',
         'max_depth': -1,
         'learning_rate': 0.001,
         "boosting": "gbdt",
         "feature_fraction": 0.91,
         "bagging_freq": 1,
         "bagging_fraction": 0.91,
         "bagging_seed": 42,
         "metric": 'mae',
         "lambda_l1": 0.1,
         "verbosity": -1,
         "nthread": -1,
         "random_state": 42} 

训练、验证和测试（每个折叠）的代码如下：

oof = np.zeros(len(scaled_train_X))
predictions = np.zeros(len(scaled_test_X))
feature_importance_df = pd.DataFrame()
#run model
for fold_, (trn_idx, val_idx) in enumerate(folds.split(scaled_train_X,train_y.values)):
    strLog = "fold {}".format(fold_)
    print(strLog)

    X_tr, X_val = scaled_train_X.iloc[trn_idx], scaled_train_X.iloc[val_idx]
    y_tr, y_val = train_y.iloc[trn_idx], train_y.iloc[val_idx]
    model = lgb.LGBMRegressor(**params, n_estimators = 20000, n_jobs = -1)
    model.fit(X_tr, 
              y_tr, 
              eval_set=[(X_tr, y_tr), (X_val, y_val)], 
              eval_metric='mae',
              verbose=1000, 
              early_stopping_rounds=500)
    oof[val_idx] = model.predict(X_val, num_iteration=model.best_iteration_)
    #feature importance
    fold_importance_df = pd.DataFrame()
    fold_importance_df["Feature"] = train_columns
    fold_importance_df["importance"] = model.feature_importances_[:len(train_columns)]
    fold_importance_df["fold"] = fold_ + 1
    feature_importance_df = pd.concat([feature_importance_df, fold_importance_df], axis=0)
    #predictions
    predictions += model.predict(scaled_test_X, num_iteration=model.best_iteration_) / folds.n_splits 

我们将预测向量（其维度与提交文件的维度相同，即每个测试段一个条目）初始化为零。我们还初始化了一个超出折叠的向量（训练数据的长度，即训练段的数量）。

对于每个折叠，我们为训练集和验证集采样数据子集和目标值。然后，我们将它们用作模型的输入，LGBMRegressor，使用之前定义的模型参数初始化（我们还将估计器的数量和工人的数量添加进去）。我们使用当前折叠对应的训练集子集来拟合模型，并使用相应的训练数据验证子集进行验证。

我们还设置了评估指标：mae——表示平均绝对误差——打印评估错误的频率，以及早期停止的迭代次数。此参数（早期停止的迭代次数）控制算法在验证错误在训练期间没有改进时停止等待的步数。我们在oof（超出折叠）向量中累积验证结果，并将当前折叠的特征重要性向量连接到特征重要性 DataFrame 中。早期停止用于在训练期间保持最佳模型进行预测。

特征重要性将用于在特征工程、特征选择和模型训练的迭代过程中观察——如果对于使用某些特征训练的模型，特征重要性在折叠间没有高变化。在每个折叠中，我们也在运行整个测试集的预测，使用每个折叠训练和验证的模型。然后，我们将每个折叠的值加到预测向量中，除以折叠数。这相当于一个模型集成，其中每个模型都使用对应于每个折叠分割的不同数据子集进行训练。

当评估当前模型时，我们将检查三块信息：训练和验证错误、这些错误在折叠间的变化，以及特征重要性在折叠间的变化。理想情况下，这些变化——训练和验证错误在折叠间的变化，以及特征重要性在折叠间的变化——应该更小。图 8.8显示了基线模型训练的评估图。

如您所见，我们在每 1,000 步（详细程度设置为1000）处绘制训练进度，并在 500 步时实施早期停止（如果验证错误在最后 500 次迭代中没有改进，则停止训练）。保留最佳模型（就验证错误而言）用于测试预测。测试预测是五个分割的平均值。

图 8.8：模型训练评估输出 – 训练和验证错误

图 8.9显示了特征重要性图，包括每个折的平均值和标准差：

图 8.9：特征重要性：平均和标准差值（来自五个折的值） – 前 20 个特征

在我们进行数据分析之后，我们继续构建时间聚合特征，这些特征是在与测试集相同持续时间的训练集子集上形成的。使用具有工程特征的新的数据集，我们训练了一个基线模型。对于基线模型，我们使用了五折交叉验证，并使用每个折的模型来预测测试集。最终预测是通过平均每个折的预测来形成的。

摘要

在本章中，我们深入探讨了处理信号数据的方法，特别是音频信号。我们探讨了此类数据的各种存储格式，并检查了用于加载、转换和可视化此类数据类型的库。为了开发强大的特征，我们应用了一系列信号处理技术。我们的特征工程努力将每个训练段的时间序列数据转换为特征，并为每个测试集汇总特征。

我们将所有特征工程过程合并为一个单一函数，适用于所有训练段和测试集。转换后的特征经过了缩放。然后我们使用这些准备好的数据，利用 LGBMRegressor 算法训练了一个基线模型。该模型采用了交叉验证，我们使用每个折训练的模型来生成测试集的预测。随后，我们将这些预测汇总以创建提交文件。此外，我们还捕获并可视化了每个折的特征重要性。

参考文献

1. LANL 地震预测，你能预测即将到来的实验室地震吗？Kaggle 竞赛：www.kaggle.com/competitions/LANL-Earthquake-Prediction

2. Gabriel Preda, LANL 地震 EDA 和预测：www.kaggle.com/code/gpreda/lanl-earthquake-eda-and-prediction

3. LANL 地震预测，数据集描述：www.kaggle.com/competitions/LANL-Earthquake-Prediction/data

4. BirdCLEF 2021 - 鸟鸣识别，在声音景观录音中识别鸟鸣，Kaggle 竞赛：www.kaggle.com/competitions/birdclef-2021

5. McFee, Brian, Colin Raffel, Dawen Liang, Daniel PW Ellis, Matt McVicar, Eric Battenberg, 和 Oriol Nieto. “librosa: Python 中的音频和音乐信号分析.” 在第 14 届 Python 科学会议论文集中，第 18-25 页。2015

6. librosa 加载函数：librosa.org/doc/main/generated/librosa.load.html

7. 康奈尔鸟鸣识别，Kaggle 比赛：www.kaggle.com/competitions/birdsong-recognition

8. EarthData MERRA2 CO，地球数据 NASA 卫星测量，Kaggle 数据集：www.kaggle.com/datasets/gpreda/earthdata-merra2-co

9. 加布里埃尔·普雷达，EARTHDATA-MERRA2 数据探索，Kaggle 笔记本：www.kaggle.com/code/gpreda/earthdata-merra2-data-exploration

加入我们书籍的 Discord 空间

加入我们的 Discord 社区，与志同道合的人相聚，并在以下地点与超过 5000 名成员一起学习：

packt.link/kaggle