如何开发用于空气污染的多变量多步时间序列预测模型

现实世界的时间序列预测之所以具有挑战性，原因有很多，不仅限于问题特征，例如存在多个输入变量、需要预测多个时间步长以及需要针对多个物理站点执行相同类型的预测。

EMC 数据科学全球黑客马拉松数据集，简称“空气质量预测”数据集，描述了多个地点的天气状况，并需要预测未来三天内的空气质量测量值。

机器学习算法可以应用于时间序列预测问题，并提供处理具有噪声复杂依赖关系的多个输入变量的能力等优势。

在本教程中，您将学习如何开发机器学习模型，用于空气污染的多步时间序列预测。

完成本教程后，您将了解：

• 如何估算缺失值并转换时间序列数据，以便可以通过监督学习算法对其进行建模。
• 如何开发和评估一套用于多步时间序列预测的线性算法。
• 如何开发和评估一套用于多步时间序列预测的非线性算法。

通过我的新书《时间序列预测深度学习》启动您的项目，其中包括分步教程和所有示例的 Python 源代码文件。

让我们开始吧。

• 2019 年 6 月更新：更新了 numpy.load()，设置 allow_pickle=True。

如何开发用于多变量多步空气污染时间序列预测的机器学习模型
图片由 Eric Schmuttenmaer 拍摄，保留部分权利。

教程概述

本教程分为九个部分；它们是

1. 问题描述
2. 模型评估
3. 机器学习建模
4. 机器学习数据准备
5. 模型评估测试工具
6. 评估线性算法
7. 评估非线性算法
8. 调整滞后大小

问题描述

空气质量预测数据集描述了多个站点的天气状况，需要预测未来三天内的空气质量测量值。

具体来说，提供了多个站点八天内每小时的温度、压力、风速和风向等天气观测数据。目标是预测多个站点未来 3 天的空气质量测量值。预测提前期不连续；相反，必须在 72 小时预测期内预测特定的提前期。它们是

此外，数据集被划分为不相交但连续的数据块，其中包含八天的数据，然后是需要预测的三天。

并非所有站点或数据块都提供所有观测数据，也并非所有站点或数据块都提供所有输出变量。存在大量缺失数据，必须加以解决。

该数据集在 2012 年 Kaggle 网站上作为短期机器学习竞赛（或黑客马拉松）的基础。

竞赛提交作品根据从参与者那里保留的真实观测数据进行评估，并使用平均绝对误差 (MAE) 进行评分。提交作品要求在由于数据缺失而无法进行预测的情况下指定值 -1,000,000。事实上，提供了一个用于插入缺失值的模板，并且要求所有提交作品都必须采用该模板（真麻烦）。

一位获胜者在保留的测试集（私有排行榜）上使用滞后观测值的随机森林实现了 0.21058 的 MAE。该解决方案的详细说明可在以下帖子中找到

• 将所有内容扔进随机森林：本·哈姆纳谈赢得空气质量预测黑客马拉松, 2012.

在本教程中，我们将探讨如何为该问题开发朴素预测，这些预测可用作基线，以确定模型是否在该问题上具有技能。

模型评估

在评估朴素预测方法之前，我们必须开发一个测试工具。

这至少包括如何准备数据以及如何评估预测。

加载数据集

第一步是下载数据集并将其加载到内存中。

可以从 Kaggle 网站免费下载数据集。您可能需要创建帐户并登录才能下载数据集。

将整个数据集（例如“下载所有”）下载到您的工作站，并在当前工作目录中解压缩存档到名为“AirQualityPrediction”的文件夹中。

• EMC 数据科学全球黑客马拉松（空气质量预测）数据

我们的重点将放在“TrainingData.csv”文件，该文件包含训练数据集，特别是以数据块形式存在的数据，其中每个数据块是八个连续的观测日和目标变量。

我们可以使用 Pandas 的 read_csv() 函数将数据文件加载到内存中，并指定第 0 行作为标题行。

我们可以按“chunkID”变量（列索引 1）对数据进行分组。

首先，让我们获取唯一块标识符的列表。

然后，我们可以收集每个块标识符的所有行，并将它们存储在字典中以便于访问。

下面定义了一个名为 to_chunks() 的函数，该函数接受加载数据的 NumPy 数组，并返回一个从 chunk_id 到块行的字典。

加载数据集并将其拆分为块的完整示例如下所示。

运行示例会打印数据集中块的数量。

数据准备

现在我们知道如何加载数据并将其拆分为块，我们可以将其分成训练和测试数据集。

每个块包含八天每小时的观测数据，尽管每个块中实际观测的数量可能差异很大。

我们可以将每个块的前五天的观测数据用于训练，后三天用于测试。

每个观测值都有一个名为“position_within_chunk”的行，其值从 1 到 192（8 天 * 24 小时）不等。因此，我们可以将此列中值小于或等于 120 (5 * 24) 的所有行作为训练数据，将大于 120 的任何值作为测试数据。

此外，任何在训练或测试拆分中没有观测值的块都可以被视为不可用而丢弃。

在使用朴素模型时，我们只对目标变量感兴趣，而不是任何输入气象变量。因此，我们可以删除输入数据，使训练和测试数据仅由每个块的 39 个目标变量以及块内位置和观测小时组成。

下面的 split_train_test() 函数实现了此行为；给定一个块字典，它将每个块拆分为训练和测试块数据的列表。

我们不需要整个测试数据集；相反，我们只需要三天内特定提前期处的观测数据，特别是以下提前期

其中，每个提前期都相对于训练期结束。

首先，我们可以将这些提前期放入一个函数中，以便于引用

接下来，我们可以将测试数据集缩减为仅包含首选提前期的数据。

我们可以通过查看“position_within_chunk”列，并使用提前期作为训练数据集结束的偏移量来实现，例如 120 + 1、120 + 2 等。

如果在测试集中找到匹配的行，则保存该行，否则生成一行 NaN 观测值。

下面的 to_forecasts() 函数实现了此功能，并返回一个 NumPy 数组，其中每个块的每个预测提前期都有一行。

我们可以将所有这些结合起来，将数据集拆分为训练集和测试集，并将结果保存到新文件中。

完整的代码示例如下所示。

运行示例首先会提示块 69 因数据不足而被从数据集中删除。

然后我们可以看到，训练集和测试集中各有 42 列，其中一列用于块 ID，一列用于块内位置，一列用于一天中的小时，以及 39 个训练变量。

我们还可以看到，测试数据集的版本明显较小，只包含预测提前期处的行。

新的训练和测试数据集分别保存在“naive_train.csv”和“naive_test.csv”文件中。

预测评估

一旦做出预测，就需要对其进行评估。

评估预测时，采用更简单的格式会很有帮助。例如，我们将使用 [块][变量][时间] 的三维结构，其中变量是目标变量编号（从 0 到 38），时间是提前期索引（从 0 到 9）。

模型预计将以这种格式进行预测。

我们还可以重新构造测试数据集，使其具有此数据集以进行比较。下面的 prepare_test_forecasts() 函数实现了此功能。

我们将使用平均绝对误差 (MAE) 来评估模型。这是竞赛中使用的指标，考虑到目标变量的非高斯分布，这是一个明智的选择。

如果某个提前期在测试集中不包含数据（例如 NaN），则不会计算该预测的误差。如果该提前期在测试集中有数据但在预测中没有数据，则观测值的全部大小将被视为误差。最后，如果测试集有观测值并且进行了预测，则绝对差将被记录为误差。

calculate_error() 函数实现了这些规则并返回给定预测的误差。

误差在所有块和所有提前期上求和，然后求平均值。

将计算总体 MAE，但我们还将计算每个预测提前期的 MAE。这通常有助于模型选择，因为某些模型在不同的提前期可能表现不同。

下面的 evaluate_forecasts() 函数实现了此功能，它计算给定预测和预期值（以 [块][变量][时间] 格式）的 MAE 和每个提前期 MAE。

一旦我们评估了模型，我们就可以展示它。

下面的 summarize_error() 函数首先打印模型性能的一行摘要，然后创建每个预测提前期的 MAE 图。

现在我们准备开始探索朴素预测方法的性能。

机器学习建模

该问题可以通过机器学习进行建模。

大多数机器学习模型不直接支持时间上的观测概念。相反，滞后观测必须被视为输入特征才能进行预测。

这是机器学习算法用于时间序列预测的优势。具体来说，它们能够支持大量输入特征。这些可以是单个或多个输入时间序列的滞后观测值。

机器学习算法用于时间序列预测相对于经典方法的其他一般优势包括：

• 能够支持嘈杂的特征和变量之间的关系中的噪声。
• 能够处理不相关的特征。
• 能够支持变量之间复杂的关系。

该数据集的一个挑战是需要进行多步预测。机器学习方法可以用于进行多步预测的两种主要方法；它们是：

• 直接。为预测每个预测提前期开发一个单独的模型。
• 递归。开发一个单一模型进行一步预测，并递归使用该模型，其中先前的预测用作输入以预测后续的提前期。

当预测一个短的连续提前期块时，递归方法可能更合理，而当预测不连续的提前期时，直接方法可能更有意义。考虑到我们有兴趣预测三天内 10 个连续和不连续提前期的混合，直接方法可能更适合空气污染预测问题。

数据集有 39 个目标变量，我们为每个目标变量、每个预测提前期开发一个模型。这意味着我们需要 (39 * 10) 390 个机器学习模型。

机器学习算法用于时间序列预测的关键是输入数据的选择。我们可以考虑三种主要的数据来源，它们可以用作输入并映射到目标变量的每个预测提前期；它们是：

• 单变量数据，例如来自正在预测的目标变量的滞后观测值。
• 多变量数据，例如来自其他变量（天气和目标）的滞后观测值。
• 元数据，例如关于预测日期或时间的数据。

数据可以从所有块中提取，为学习从输入到目标预测提前期的映射提供丰富的数据集。

39 个目标变量实际上由 14 个站点上的 12 个变量组成。

由于数据提供的方式，默认的建模方法是将每个变量-站点视为独立的。有可能按变量折叠数据，并在多个站点对同一变量使用相同的模型。

一些变量被故意错误标记（例如，不同的数据使用具有相同标识符的变量）。尽管如此，也许可以识别这些错误标记的变量并将其排除在多站点模型之外。

机器学习数据准备

在探索该数据集的机器学习模型之前，我们必须以可以拟合模型的方式准备数据。

这需要两个数据准备步骤：

• 处理缺失数据。
• 准备输入-输出模式。

目前，我们将专注于 39 个目标变量，忽略气象和元数据。

处理缺失数据

块由五天或更少的每小时观测数据组成，针对 39 个目标变量。

许多块没有所有五天的数据，也没有任何块包含所有 39 个目标变量的数据。

在某个块缺少目标变量数据的情况下，不需要进行预测。

在某个块确实包含某个目标变量的一些数据，但没有全部五天数据的情况下，序列中将存在间隙。这些间隙可能是几个小时到一天多的观测值，有时甚至更长。

处理这些间隙的三个候选策略如下：

• 忽略间隙。
• 使用无间隙数据。
• 填补间隙。

我们可以忽略这些空白。这样做的问题是，当将数据拆分为输入和输出时，数据将不连续。当训练模型时，输入将不一致，但可能意味着最近 n 小时的数据，或分散在最近 n 天的数据。这种不一致性将使从输入到输出的学习映射非常嘈杂，并且可能比模型所需的时间更困难。

我们只能使用没有间隙的数据。这是一个不错的选择。风险在于我们可能没有足够的数据来拟合模型。

最后，我们可以填补空白。这称为数据插补，有许多策略可以用来填补空白。三种可能表现良好的方法包括：

• 将最后观察到的值向前延续（线性）。
• 使用该块内一天中该小时的中位数。
• 使用跨块一天中该小时的中位数。

在本教程中，我们将使用后一种方法，通过使用跨块的一天中该小时的中位数来填补空白。经过一些测试后，这种方法似乎能带来更多的训练样本和更好的模型性能。

对于给定的变量，可能存在由缺失行定义的缺失观测值。具体来说，每个观测值都有一个“position_within_chunk”。我们期望训练数据集中的每个块都有 120 个观测值，其中“positions_within_chunk”从 1 到 120（包括）。

因此，我们可以为每个变量创建一个包含 120 个 NaN 值的数组，使用“positions_within_chunk”值标记块中的所有观测值，任何剩余的都将标记为 NaN。然后我们可以绘制每个变量并查找间隙。

下面的 variable_to_series() 函数将获取一个块的行和目标变量的给定列索引，并将返回一个包含 120 个时间步长的变量序列，其中所有可用数据都使用块中的值进行标记。

我们需要为每个块计算一个并行的日期小时序列，我们可以用它来为块中的每个变量估算特定小时的数据。

给定一个部分填充的日期小时序列，下面的 interpolate_hours() 函数将填充缺失的日期小时。它通过找到第一个标记的小时，然后向前计数，填充日期小时，然后反向执行相同的操作来实现。

我们可以调用相同的 variable_to_series() 函数（上方）来创建带有缺失值的小时序列（列索引 2），然后调用 interpolate_hours() 来填充缺失值。

然后我们可以将小时传递给任何可能使用它的插补函数。

现在我们可以尝试使用同一序列中相同小时的值来填充块中的缺失值。具体来说，我们将找到序列中所有具有相同小时的行，并计算中位数。

下面的 impute_missing() 函数接受块中的所有行、为块准备好的日期小时序列、具有缺失值的变量序列以及变量的列索引。

它首先检查序列是否完全缺失数据，如果是，则立即返回，因为无法执行插补。然后，它遍历序列的时间步长，当检测到没有数据的时间步长时，它收集序列中具有相同小时数据的所有行并计算中位数。

监督表示

我们需要将每个目标变量的序列转换为包含输入和输出的行，以便我们可以拟合监督机器学习算法。

具体来说，我们有一个序列，例如：

当使用 2 个滞后变量预测提前期 +1 时，我们将序列拆分为输入 (X) 和输出 (y) 模式，如下所示：

这首先要求我们选择作为输入使用的滞后观测数量。没有正确答案；相反，最好测试不同的数量并查看哪种有效。

然后，我们必须将序列拆分为监督学习格式，用于 10 个预测提前期中的每一个。例如，使用 2 个滞后观测值预测 +24 可能如下所示：

然后为 39 个目标变量中的每一个重复此过程。

然后可以将为每个提前期、每个目标变量准备的模式跨块聚合，以提供模型的训练数据集。

我们还必须准备一个测试数据集。也就是说，每个目标变量的每个块的输入数据 (X)，以便我们可以将其用作输入来预测测试数据集中的提前期。如果我们将滞后设置为 2，那么测试数据集将由每个目标变量的每个块的最后两个观测值组成。非常简单。

我们可以从定义一个函数开始，该函数将为给定的完整（插补）序列创建输入-输出模式。

下面的 supervised_for_lead_time() 函数将接受一个序列，一个用作输入的滞后观测数量，以及一个要预测的预测提前期，然后将返回从序列中提取的输入/输出行列表。

理解这一部分很重要。

我们可以在小型测试数据集上测试此函数，并探索不同数量的滞后变量和预测提前期。

以下是一个完整的示例，它生成一个包含20个整数的序列，并创建了一个具有两个输入滞后时间的序列，用于预测+6的提前期。

运行示例将打印出生成的模式，显示滞后观测值及其相关的预测提前期。

请尝试此示例，以熟悉此数据转换，因为它对于使用机器学习算法建模时间序列至关重要。

我们现在可以为给定目标变量序列的每个预测提前期调用 supervised_for_lead_time() 函数。

下面的 target_to_supervised() 函数实现了这一点。首先，目标变量被转换为序列，并使用上一节中开发的函数进行插补。然后为每个目标提前期创建训练样本。还为目标变量创建了一个测试样本。

然后返回此目标变量的每个预测提前期的训练数据和测试输入数据。

我们已经有了这些部分；现在我们需要定义一个函数来驱动数据准备过程。

此函数构建训练和测试数据集。

方法是枚举每个目标变量，并从所有块中收集每个提前期的训练数据。同时，我们收集在对测试数据集进行预测时作为输入所需的样本。

结果是一个训练数据集，其维度为 [变量][提前期][样本]，其中最后一个维度是目标变量的预测提前期的训练样本行。该函数还返回维度为 [块][变量][样本] 的测试数据集，其中最后一个维度是用于对块的目标变量进行预测的输入数据。

下面的 data_prep() 函数实现了此行为，它以块格式获取数据和指定数量的滞后观测值作为输入。

完整示例

我们可以将所有内容连接起来，并准备一个具有监督学习格式的训练和测试数据集，用于机器学习算法。

我们将使用前12小时的滞后观测值作为输入来预测每个预测提前期。

然后将生成的训练和测试数据集保存为二进制NumPy数组。

完整的示例如下所示。

运行示例可能需要一分钟。

结果是两个二进制文件，包含训练和测试数据集，我们可以在以下部分加载它们以训练和评估机器学习算法。

模型评估测试工具

在开始评估算法之前，我们需要测试工具的更多元素。

首先，我们需要能够根据训练数据拟合scikit-learn模型。下面的 fit_model() 函数将克隆模型配置并根据提供的训练数据进行拟合。我们将需要拟合许多（360个）每个配置模型的版本，因此此函数将被大量调用。

接下来，我们需要为每个变量和预测提前期组合拟合一个模型。

我们可以通过首先按变量，然后按提前期枚举训练数据集来完成此操作。然后，我们可以拟合一个模型并将其存储在一个具有相同结构的列表列表中，具体为： [变量][时间][模型]。

下面的 fit_models() 函数实现了这一点。

模型拟合是缓慢的部分，可以通过并行化（例如使用Joblib库）来提高效率。这留作扩展。

一旦模型拟合好，它们就可以用来对测试数据集进行预测。

准备好的测试数据集首先按块组织，然后按目标变量组织。进行预测很快，首先需要检查是否可以进行预测（我们有输入数据），如果可以，则使用目标变量的相应模型。然后，变量的每个10个预测提前期将使用这些提前期的每个直接模型进行预测。

下面的 make_predictions() 函数实现了这一点，它将模型列表列表和加载的测试数据集作为参数，并返回一个结构为 [块][变量][时间] 的预测数组。

我们需要一个模型列表进行评估。

我们可以定义一个通用的 get_models() 函数，它负责定义一个模型名称映射到配置好的 scikit-learn 模型对象的字典。

最后，我们需要一个函数来驱动模型评估过程。

给定模型字典，枚举模型，首先在训练数据上拟合模型矩阵，对测试数据集进行预测，评估预测结果，并总结结果。

下面的 evaluate_models() 函数实现了这一点。

我们现在拥有评估机器学习模型所需的一切。

评估线性算法

在本节中，我们将抽查一组线性机器学习算法。

线性算法假定输出是输入变量的线性函数。这与ARIMA等经典时间序列预测模型的假设非常相似。

抽查意味着评估一系列模型以粗略了解哪些模型有效。我们对任何优于简单自回归模型AR(2)（其MAE误差约为0.487）的模型感兴趣。

我们将测试八种具有默认配置的线性机器学习算法；具体来说：

• 线性回归
• Lasso 线性回归
• 岭回归
• 弹性网络回归
• Huber 回归
• Lasso Lars 线性回归
• Passive Aggressive 回归
• 随机梯度下降回归

我们可以在 get_models() 函数中定义这些模型。

完整的代码示例如下所示。