统计、机器学习和神经网络预测方法

在这个笔记本中，您将为M5数据集进行预测，使用交叉验证选择每个时间序列的最佳模型。

统计、机器学习和神经网络预测方法
在本教程中，我们将通过利用最适合每个时间序列的模型来探索对M5数据集进行预测的过程。我们将通过一种称为交叉验证的基本技术来实现这一目标。这种方法帮助我们估计模型的预测性能，并选择对每个时间序列最佳的模型。

M5数据集包含来自沃尔玛的五年期层级销售数据。目标是预测未来28天的每日销售额。该数据集划分为美国的50个州，每个州有10家商店。

在时间序列预测和分析领域，识别最适合特定系列组的模型是一个较为复杂的任务。通常，这一选择过程在很大程度上依赖于直觉，而这种直觉不一定与我们的数据集的实际情况相符。

在本教程中，我们旨在为M5基准数据集中不同系列组的模型选择提供一种更结构化的数据驱动方法。该数据集在预测领域颇为知名，使我们能够展示我们方法论的多样性和强大功能。

我们将训练来自各种预测范式的多种模型：

StatsForecast

MLForecast

机器学习：利用 LightGBM、XGBoost 和 LinearRegression 等机器学习模型可以带来优势，因为它们能够揭示数据中的复杂模式。我们将为此目的使用MLForecast库。

NeuralForecast

深度学习：深度学习模型，如变换器（AutoTFT）和神经网络（AutoNHITS），使我们能够处理时间序列数据中的复杂非线性依赖关系。我们将利用NeuralForecast库来处理这些模型。

使用Nixtla库套件，我们将能够以数据驱动的方式进行模型选择，确保我们为数据集中特定系列组使用最合适的模型。

大纲：

Warning

本教程最初是在一个 c5d.24xlarge EC2 实例上执行的。

安装库

%%capture
!pip install statsforecast mlforecast neuralforecast datasetforecast s3fs pyarrow

下载和准备数据

这个示例使用了M5数据集。它由30,490条底层时间序列组成。

import pandas as pd
# 从提供的URL加载训练目标数据集
Y_df = pd.read_parquet('https://m5-benchmarks.s3.amazonaws.com/data/train/target.parquet')

# Rename columns to match the Nixtlaverse's expectations
# The 'item_id' becomes 'unique_id' representing the unique identifier of the time series
# The 'timestamp' becomes 'ds' representing the time stamp of the data points
# The 'demand' becomes 'y' representing the target variable we want to forecast
Y_df = Y_df.rename(columns={
    'item_id': 'unique_id', 
    'timestamp': 'ds', 
    'demand': 'y'
})

# Convert the 'ds' column to datetime format to ensure proper handling of date-related operations in subsequent steps
Y_df['ds'] = pd.to_datetime(Y_df['ds'])

为了简单起见，我们只保留一个类别。

Y_df = Y_df.query('unique_id.str.startswith("FOODS_3")').reset_index(drop=True)

Y_df['unique_id'] = Y_df['unique_id'].astype(str)

基本绘图

使用 StatsForecast 类中的 plot 方法绘制一些序列。此方法从数据集中随机打印 8 个序列，非常适合基本的 EDA。

from statsforecast import StatsForecast
/home/ubuntu/statsforecast/statsforecast/core.py:23: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from tqdm.autonotebook import tqdm
# 功能：为探索性数据分析绘制随机序列
StatsForecast.plot(Y_df)
# 功能：绘制一系列组的图表以进行探索性数据分析
StatsForecast.plot(Y_df, unique_ids=["FOODS_3_432_TX_2"])
# 功能：绘制系列组以进行探索性数据分析
StatsForecast.plot(Y_df, unique_ids=["FOODS_3_432_TX_2"], engine ='matplotlib')

使用统计、机器学习和神经网络方法创建预测。

StatsForecast

StatsForecast 是一个全面的库，提供了一套流行的单变量时间序列预测模型，所有模型都专注于高性能和可扩展性。

以下是 StatsForecast 成为时间序列预测强大工具的原因：

  • 本地模型集合：StatsForecast 提供了多样的本地模型，可以单独应用于每个时间序列，从而捕捉每个序列中的独特模式。

  • 简单性：使用 StatsForecast，训练、预测和回测多个模型变得非常简单，仅需几行代码即可完成。这种简单性使其成为初学者和经验丰富的从业者的便捷工具。

  • 速度优化：StatsForecast 中模型的实现经过速度优化，确保大规模计算高效进行，从而减少模型训练和预测的整体时间。

  • 水平可扩展性：StatsForecast 的一个显著特点是其水平可扩展性。它兼容分布式计算框架，如 Spark、Dask 和 Ray。此特性使其能够通过在集群中的多个节点间分配计算来处理大规模数据集，因此成为大规模时间序列预测任务的首选解决方案。

StatsForecast 接收一个模型列表来拟合每个时间序列。由于我们处理的是日数据，使用7作为季节性会比较有利。

# 从 statsforecast 库中导入必要的模型
from statsforecast.models import (
    # 季节性朴素法：一种使用上一季节数据作为预测的模型。
    SeasonalNaive,
    # 朴素法：一种简单的模型，使用最近观测到的值作为预测值。
    Naive,
    # 历史平均值：该模型使用所有历史数据的平均值作为预测值。
    HistoricAverage,
    # Croston优化：一种专为间歇性需求预测设计的模型
    CrostonOptimized,
    # ADIDA：间歇性需求方法的自适应组合，一种专为间歇性需求设计的模型
    ADIDA,
    # IMAPA：间歇性乘法自回归平均模型，一种用于间歇性序列的模型，该模型结合了自相关性。
    IMAPA,
    # AutoETS：自动指数平滑模型，该模型基于AIC自动选择最佳的指数平滑模型。
    AutoETS
)

我们通过实例化一个新的 StatsForecast 对象来拟合模型，使用以下参数：

  • models：一个模型列表。从模型中选择所需的模型并导入它们。
  • freq：一个字符串，指示数据的频率。（请参见 pandas 的可用频率。）
  • n_jobs：整数，用于并行处理的作业数量，使用 -1 表示所有核心。
  • fallback_model：在模型失败时要使用的模型。 任何设置都会传递到构造函数中。然后调用其 fit 方法并传入历史数据框。
horizon = 28
models = [
    SeasonalNaive(season_length=7),
    Naive(),
    HistoricAverage(),
    CrostonOptimized(),
    ADIDA(),
    IMAPA(),
    AutoETS(season_length=7)
]
# 实例化 StatsForecast 类
sf = StatsForecast(
    models=models,  # 用于预测的模型列表
    freq='D',  # The frequency of the time series data (in this case, 'D' stands for daily frequency)
    n_jobs=-1,  # 用于并行执行的CPU核心数量（-1表示使用所有可用核心）
)

预测方法接受两个参数：预测下一个 h（时间范围）和水平。

  • h（整数）：表示预测未来 h 步的值。在这种情况下，是 12 个月。
  • level（浮点数列表）：此可选参数用于概率预测。设置预测区间的水平（或置信百分位）。例如，level=[90] 意味着模型预计真实值在该区间内的概率为 90%。

这里的预测对象是一个新的数据框，包含一个列用于模型名称和 y hat 值，以及不确定性区间的列。

这段代码用于计算运行 StatsForecast 类的预测函数所需的时间，该函数预测未来 28 天（h=28）。置信水平设置为 [90]，这意味着它将计算 90% 的预测区间。时间以分钟为单位计算，并在最后打印出来。

from time import time

# 在预测开始前获取当前时间，这将用于测量执行时间。
init = time()

# 调用StatsForecast实例的forecast方法，预测未来28天（h=28） 
# 置信水平设为[90]，这意味着将计算90%的预测区间。
fcst_df = sf.forecast(df=Y_df, h=28, level=[90])

# 预测结束后获取当前时间
end = time()

# 计算并打印预测所花费的总时间（以分钟为单位）
print(f'Forecast Minutes: {(end - init) / 60}')

Forecast Minutes: 2.270755163828532
fcst_df.head()
ds SeasonalNaive SeasonalNaive-lo-90 SeasonalNaive-hi-90 Naive Naive-lo-90 Naive-hi-90 HistoricAverage HistoricAverage-lo-90 HistoricAverage-hi-90 CrostonOptimized ADIDA IMAPA AutoETS AutoETS-lo-90 AutoETS-hi-90
unique_id
FOODS_3_001_CA_1 2016-05-23 1.0 -2.847174 4.847174 2.0 0.098363 3.901637 0.448738 -1.009579 1.907055 0.345192 0.345477 0.347249 0.381414 -1.028122 1.790950
FOODS_3_001_CA_1 2016-05-24 0.0 -3.847174 3.847174 2.0 -0.689321 4.689321 0.448738 -1.009579 1.907055 0.345192 0.345477 0.347249 0.286933 -1.124136 1.698003
FOODS_3_001_CA_1 2016-05-25 0.0 -3.847174 3.847174 2.0 -1.293732 5.293732 0.448738 -1.009579 1.907055 0.345192 0.345477 0.347249 0.334987 -1.077614 1.747588
FOODS_3_001_CA_1 2016-05-26 1.0 -2.847174 4.847174 2.0 -1.803274 5.803274 0.448738 -1.009579 1.907055 0.345192 0.345477 0.347249 0.186851 -1.227280 1.600982
FOODS_3_001_CA_1 2016-05-27 0.0 -3.847174 3.847174 2.0 -2.252190 6.252190 0.448738 -1.009579 1.907055 0.345192 0.345477 0.347249 0.308112 -1.107548 1.723771

MLForecast

MLForecast 是一个强大的库，提供了时间序列预测的自动特征创建，促进了全局机器学习模型的使用。它专为高性能和可扩展性而设计。

MLForecast 的主要功能包括：

  • 支持 sklearn 模型：MLForecast 兼容遵循 scikit-learn API 的模型。这使得它具有高度的灵活性，能够与各种机器学习算法无缝集成。

  • 简单性：使用 MLForecast，可以用简短的代码完成模型的训练、预测和回测任务。这种简化的操作使其对各个级别的从业者都非常友好。

  • 优化速度：MLForecast 旨在快速执行任务，这在处理大数据集和复杂模型时至关重要。

  • 水平可扩展性：MLForecast 能够通过使用分布式计算框架如 Spark、Dask 和 Ray 实现水平扩展。此功能使其能够通过将计算分配到集群中的多个节点上高效处理大规模数据集，非常适合大规模时间序列预测任务。

from mlforecast import MLForecast
from mlforecast.target_transforms import Differences
from mlforecast.utils import PredictionIntervals
from window_ops.expanding import expanding_mean
%%capture
!pip install lightgbm xgboost
# 从各个库中导入必要的模型

# LGBMRegressor：一种使用LightGBM库中基于树的学习算法的梯度提升框架
from lightgbm import LGBMRegressor

# XGBRegressor：来自XGBoost库的梯度提升回归模型
from xgboost import XGBRegressor

# 线性回归：来自scikit-learn库的一个简单线性回归模型
from sklearn.linear_model import LinearRegression

要使用 MLForecast 进行时间序列预测，我们实例化一个新的 MLForecast 对象，并为其提供各种参数，以便将建模过程调整为满足我们的特定需求：

  • models：该参数接受一个您希望用于预测的机器学习模型列表。您可以从 scikit-learn、lightgbm 和 xgboost 导入您喜欢的模型。

  • freq：这是一个字符串，表示您的数据的频率（每小时、每日、每周等）。该字符串的特定格式应与 pandas 识别的频率字符串对齐。

  • target_transforms：这些是应用于目标变量的转换，发生在模型训练之前和模型预测之后。当处理可能受益于转换的数据时，例如对高度偏斜数据进行对数转化，这可能会很有用。

  • lags：该参数接受特定的滞后值，用作回归变量。滞后表示您在为模型创建特征时希望回顾多久的时间。例如，如果您希望使用前一天的数据作为预测今天值的特征，则您应该指定滞后为 1。

  • lags_transforms：这些是针对每个滞后的特定转换。这允许您对滞后特征应用转换。

  • date_features：该参数指定与日期相关的特征，用作回归变量。例如，您可能希望在模型中包括星期几或月份作为特征。

  • num_threads：该参数控制用于并行化特征创建的线程数量，帮助加快处理大数据集时的速度。

所有这些设置都传递给 MLForecast 构造函数。一旦 MLForecast 对象使用这些设置初始化，我们调用其 fit 方法，并将历史数据框作为参数传递。fit 方法根据提供的历史数据训练模型，为未来的预测任务做好准备。

# 实例化 MLForecast 对象
mlf = MLForecast(
    models=[LGBMRegressor(), XGBRegressor(), LinearRegression()],  # 预测模型列表：LightGBM、XGBoost 和线性回归
    freq='D',  # Frequency of the data - 'D' for daily frequency
    lags=list(range(1, 7)),  # 作为回归变量的特定滞后期：1至6天
    lag_transforms = {
        1:  [expanding_mean],  # 对滞后1天的数据应用扩展均值变换
    },
    date_features=['year', 'month', 'day', 'dayofweek', 'quarter', 'week'],  # 用作回归变量的日期特征
)

只需调用 fit 模型以训练选择的模型。在这种情况下，我们正在生成保形预测区间。

# 启动计时器以计算模型拟合所需的时间
init = time()

# 将MLForecast模型拟合到数据中，使用28天的窗口大小设置预测区间
mlf.fit(Y_df, prediction_intervals=PredictionIntervals(window_size=28))

# 计算模型拟合后的结束时间
end = time()

# 打印拟合MLForecast模型所花费的时间，以分钟为单位
print(f'MLForecast Minutes: {(end - init) / 60}')

MLForecast Minutes: 2.2809854547182717

之后，只需调用 predict 来生成预测结果。

fcst_mlf_df = mlf.predict(28, level=[90])
fcst_mlf_df.head()
unique_id ds LGBMRegressor XGBRegressor LinearRegression LGBMRegressor-lo-90 LGBMRegressor-hi-90 XGBRegressor-lo-90 XGBRegressor-hi-90 LinearRegression-lo-90 LinearRegression-hi-90
0 FOODS_3_001_CA_1 2016-05-23 0.549520 0.598431 0.359638 -0.213915 1.312955 -0.020050 1.216912 0.030000 0.689277
1 FOODS_3_001_CA_1 2016-05-24 0.553196 0.337268 0.100361 -0.251383 1.357775 -0.201449 0.875985 -0.216195 0.416917
2 FOODS_3_001_CA_1 2016-05-25 0.599668 0.349604 0.175840 -0.203974 1.403309 -0.284416 0.983624 -0.150593 0.502273
3 FOODS_3_001_CA_1 2016-05-26 0.638097 0.322144 0.156460 0.118688 1.157506 -0.085872 0.730160 -0.273851 0.586771
4 FOODS_3_001_CA_1 2016-05-27 0.763305 0.300362 0.328194 -0.313091 1.839701 -0.296636 0.897360 -0.657089 1.313476

神经预测

NeuralForecast是一个强大的神经预测模型集合，专注于可用性和性能。它包括多种模型架构，从经典网络如多层感知机（MLP）和循环神经网络（RNN）到新颖的贡献，如N-BEATS、N-HITS、时序融合变压器（TFT）等。

NeuralForecast的主要特点包括：

  • 广泛的全球模型集合。开箱即用的MLP、LSTM、RNN、TCN、扩张RNN、NBEATS、NHITS、ESRNN、TFT、Informer、PatchTST和HINT的实现。
  • 简单直观的接口，允许在几行代码中训练、预测和回测各种模型。
  • 支持GPU加速，以提高计算速度。

这台机器没有GPU，但Google Colab提供了一些免费的GPU。

使用Colab的GPU来训练NeuralForecast。

# 从Colab读取结果
fcst_nf_df = pd.read_parquet('https://m5-benchmarks.s3.amazonaws.com/data/forecast-nf.parquet')
fcst_nf_df.head()
unique_id ds AutoNHITS AutoNHITS-lo-90 AutoNHITS-hi-90 AutoTFT AutoTFT-lo-90 AutoTFT-hi-90
0 FOODS_3_001_CA_1 2016-05-23 0.0 0.0 2.0 0.0 0.0 2.0
1 FOODS_3_001_CA_1 2016-05-24 0.0 0.0 2.0 0.0 0.0 2.0
2 FOODS_3_001_CA_1 2016-05-25 0.0 0.0 2.0 0.0 0.0 1.0
3 FOODS_3_001_CA_1 2016-05-26 0.0 0.0 2.0 0.0 0.0 2.0
4 FOODS_3_001_CA_1 2016-05-27 0.0 0.0 2.0 0.0 0.0 2.0 
# 合并来自StatsForecast和NeuralForecast的预测
fcst_df = fcst_df.merge(fcst_nf_df, how='left', on=['unique_id', 'ds'])

# 将MLForecast的预测结果合并到综合预测数据框中
fcst_df = fcst_df.merge(fcst_mlf_df, how='left', on=['unique_id', 'ds'])

fcst_df.head()
unique_id ds SeasonalNaive SeasonalNaive-lo-90 SeasonalNaive-hi-90 Naive Naive-lo-90 Naive-hi-90 HistoricAverage HistoricAverage-lo-90 ... AutoTFT-hi-90 LGBMRegressor XGBRegressor LinearRegression LGBMRegressor-lo-90 LGBMRegressor-hi-90 XGBRegressor-lo-90 XGBRegressor-hi-90 LinearRegression-lo-90 LinearRegression-hi-90
0 FOODS_3_001_CA_1 2016-05-23 1.0 -2.847174 4.847174 2.0 0.098363 3.901637 0.448738 -1.009579 ... 2.0 0.549520 0.598431 0.359638 -0.213915 1.312955 -0.020050 1.216912 0.030000 0.689277
1 FOODS_3_001_CA_1 2016-05-24 0.0 -3.847174 3.847174 2.0 -0.689321 4.689321 0.448738 -1.009579 ... 2.0 0.553196 0.337268 0.100361 -0.251383 1.357775 -0.201449 0.875985 -0.216195 0.416917
2 FOODS_3_001_CA_1 2016-05-25 0.0 -3.847174 3.847174 2.0 -1.293732 5.293732 0.448738 -1.009579 ... 1.0 0.599668 0.349604 0.175840 -0.203974 1.403309 -0.284416 0.983624 -0.150593 0.502273
3 FOODS_3_001_CA_1 2016-05-26 1.0 -2.847174 4.847174 2.0 -1.803274 5.803274 0.448738 -1.009579 ... 2.0 0.638097 0.322144 0.156460 0.118688 1.157506 -0.085872 0.730160 -0.273851 0.586771
4 FOODS_3_001_CA_1 2016-05-27 0.0 -3.847174 3.847174 2.0 -2.252190 6.252190 0.448738 -1.009579 ... 2.0 0.763305 0.300362 0.328194 -0.313091 1.839701 -0.296636 0.897360 -0.657089 1.313476

5 rows × 32 columns

预测图表

sf.plot(Y_df, fcst_df, max_insample_length=28 * 3)

使用绘图函数探索模型和ID。

sf.plot(Y_df, fcst_df, max_insample_length=28 * 3, 
        models=['CrostonOptimized', 'AutoNHITS', 'SeasonalNaive', 'LGBMRegressor'])

验证模型的性能

这三个库 - StatsForecast、MLForecast 和 NeuralForecast - 提供了专门为时间序列设计的开箱即用的交叉验证功能。这使我们能够使用历史数据来评估模型的性能，从而获得对每个模型在未见数据上表现的无偏评估。

来自现代预测实践课程

StatsForecast中的交叉验证

StatsForecast类中的cross_validation方法接受以下参数：

  • df：一个表示训练数据的DataFrame。
  • h（int）：预测范围，表示我们希望预测的未来步数。例如，如果我们预测的是每小时数据，则h=24将表示24小时的预测。
  • step_size（int）：每个交叉验证窗口之间的步长。此参数确定我们希望多频繁地运行预测过程。
  • n_windows（int）：用于交叉验证的窗口数量。此参数定义了我们希望评估多少个过去的预测过程。

这些参数使我们能够控制交叉验证过程的范围和粒度。通过调整这些设置，我们能够在计算成本和交叉验证的彻底性之间取得平衡。

init = time()
cv_df = sf.cross_validation(df=Y_df, h=horizon, n_windows=3, step_size=horizon, level=[90])
end = time()
print(f'CV Minutes: {(end - init) / 60}')
/home/ubuntu/statsforecast/statsforecast/ets.py:1041: RuntimeWarning:

divide by zero encountered in double_scalars

CV Minutes: 5.206169327100118

crossvaldation_df对象是一个新的数据框，包含以下列：

  • unique_id索引: （如果您不喜欢使用索引，您可以运行forecasts_cv_df.resetindex()）
  • ds: 日期戳或时间索引
  • cutoff: n_windows的最后一个日期戳或时间索引。如果n_windows=1，则只有一个唯一的cutoff值；如果n_windows=2，则有两个唯一的cutoff值。
  • y: 真值
  • "model": 包含模型名称和拟合值的列。
cv_df.head()
ds cutoff y SeasonalNaive SeasonalNaive-lo-90 SeasonalNaive-hi-90 Naive Naive-lo-90 Naive-hi-90 HistoricAverage HistoricAverage-lo-90 HistoricAverage-hi-90 CrostonOptimized ADIDA IMAPA AutoETS AutoETS-lo-90 AutoETS-hi-90
unique_id
FOODS_3_001_CA_1 2016-02-29 2016-02-28 0.0 2.0 -1.878885 5.878885 0.0 -1.917011 1.917011 0.449111 -1.021813 1.920036 0.618472 0.618375 0.617998 0.655286 -0.765731 2.076302
FOODS_3_001_CA_1 2016-03-01 2016-02-28 1.0 0.0 -3.878885 3.878885 0.0 -2.711064 2.711064 0.449111 -1.021813 1.920036 0.618472 0.618375 0.617998 0.568595 -0.853966 1.991155
FOODS_3_001_CA_1 2016-03-02 2016-02-28 1.0 0.0 -3.878885 3.878885 0.0 -3.320361 3.320361 0.449111 -1.021813 1.920036 0.618472 0.618375 0.617998 0.618805 -0.805298 2.042908
FOODS_3_001_CA_1 2016-03-03 2016-02-28 0.0 1.0 -2.878885 4.878885 0.0 -3.834023 3.834023 0.449111 -1.021813 1.920036 0.618472 0.618375 0.617998 0.455891 -0.969753 1.881534
FOODS_3_001_CA_1 2016-03-04 2016-02-28 0.0 1.0 -2.878885 4.878885 0.0 -4.286568 4.286568 0.449111 -1.021813 1.920036 0.618472 0.618375 0.617998 0.591197 -0.835987 2.018380

MLForecast

MLForecast类中的cross_validation方法接受以下参数：

  • data：训练数据框
  • window_size（整数）：表示预测的未来h步。在这种情况下，是向前24小时。
  • step_size（整数）：每个窗口之间的步长。换句话说：您希望多频繁地运行预测过程。
  • n_windows（整数）：用于交叉验证的窗口数量。换句话说：您希望评估过去多少个预测过程。
  • prediction_intervals：用于计算合规区间的类。
init = time()
cv_mlf_df = mlf.cross_validation(
    data=Y_df, 
    window_size=horizon, 
    n_windows=3, 
    step_size=horizon, 
    level=[90],
)
end = time()
print(f'CV Minutes: {(end - init) / 60}')
/home/ubuntu/miniconda/envs/statsforecast/lib/python3.10/site-packages/mlforecast/forecast.py:576: UserWarning:

Excuting `cross_validation` after `fit` can produce unexpected errors

/home/ubuntu/miniconda/envs/statsforecast/lib/python3.10/site-packages/mlforecast/forecast.py:468: UserWarning:

Please rerun the `fit` method passing a proper value to prediction intervals to compute them.

/home/ubuntu/miniconda/envs/statsforecast/lib/python3.10/site-packages/mlforecast/forecast.py:468: UserWarning:

Please rerun the `fit` method passing a proper value to prediction intervals to compute them.

/home/ubuntu/miniconda/envs/statsforecast/lib/python3.10/site-packages/mlforecast/forecast.py:468: UserWarning:

Please rerun the `fit` method passing a proper value to prediction intervals to compute them.

CV Minutes: 2.961174162228902

crossvaldation_df 对象是一个新的数据框，其中包括以下列：

  • unique_id 索引： （如果您不喜欢使用索引，可以运行 forecasts_cv_df.resetindex()）
  • ds: 日期戳或时间索引
  • cutoff: n_windows 的最后一个日期戳或时间索引。如果 n_windows=1，则为一个唯一的 cutoff 值；如果 n_windows=2，则有两个唯一的 cutoff 值。
  • y: 真实值
  • "model": 包含模型名称和拟合值的列。
cv_mlf_df.head()
unique_id ds cutoff y LGBMRegressor XGBRegressor LinearRegression
0 FOODS_3_001_CA_1 2016-02-29 2016-02-28 0.0 0.435674 0.556261 -0.312492
1 FOODS_3_001_CA_1 2016-03-01 2016-02-28 1.0 0.639676 0.625806 -0.041924
2 FOODS_3_001_CA_1 2016-03-02 2016-02-28 1.0 0.792989 0.659650 0.263699
3 FOODS_3_001_CA_1 2016-03-03 2016-02-28 0.0 0.806868 0.535121 0.482491
4 FOODS_3_001_CA_1 2016-03-04 2016-02-28 0.0 0.829106 0.313353 0.677326

神经预测

此机器没有GPU，但Google Colab提供了一些免费的GPU。

使用Colab的GPU来训练NeuralForecast。

cv_nf_df = pd.read_parquet('https://m5-benchmarks.s3.amazonaws.com/data/cross-validation-nf.parquet')
cv_nf_df.head()
unique_id ds cutoff AutoNHITS AutoNHITS-lo-90 AutoNHITS-hi-90 AutoTFT AutoTFT-lo-90 AutoTFT-hi-90 y
0 FOODS_3_001_CA_1 2016-02-29 2016-02-28 0.0 0.0 2.0 1.0 0.0 2.0 0.0
1 FOODS_3_001_CA_1 2016-03-01 2016-02-28 0.0 0.0 2.0 1.0 0.0 2.0 1.0
2 FOODS_3_001_CA_1 2016-03-02 2016-02-28 0.0 0.0 2.0 1.0 0.0 2.0 1.0
3 FOODS_3_001_CA_1 2016-03-03 2016-02-28 0.0 0.0 2.0 1.0 0.0 2.0 0.0
4 FOODS_3_001_CA_1 2016-03-04 2016-02-28 0.0 0.0 2.0 1.0 0.0 2.0 0.0

合并交叉验证预测结果

cv_df = cv_df.merge(cv_nf_df.drop(columns=['y']), how='left', on=['unique_id', 'ds', 'cutoff'])
cv_df = cv_df.merge(cv_mlf_df.drop(columns=['y']), how='left', on=['unique_id', 'ds', 'cutoff'])

绘制交叉验证图表

cutoffs = cv_df['cutoff'].unique()
for cutoff in cutoffs:
    img = sf.plot(
        Y_df, 
        cv_df.query('cutoff == @cutoff').drop(columns=['y', 'cutoff']), 
        max_insample_length=28 * 5, 
        unique_ids=['FOODS_3_001_CA_1'],
    )
    img.show()

总需求

agg_cv_df = cv_df.loc[:,~cv_df.columns.str.contains('hi|lo')].groupby(['ds', 'cutoff']).sum(numeric_only=True).reset_index()
agg_cv_df.insert(0, 'unique_id', 'agg_demand')
agg_Y_df = Y_df.groupby(['ds']).sum(numeric_only=True).reset_index()
agg_Y_df.insert(0, 'unique_id', 'agg_demand')
for cutoff in cutoffs:
    img = sf.plot(
        agg_Y_df, 
        agg_cv_df.query('cutoff == @cutoff').drop(columns=['y', 'cutoff']),
        max_insample_length=28 * 5,
    )
    img.show()

每个系列和交叉验证窗口的评估

在这一部分，我们将评估每个模型在每个时间序列和每个交叉验证窗口中的表现。由于我们有许多组合，我们将使用 dask 来并行化评估。并行化将使用 fugue 来完成。

from typing import List, Callable

from distributed import Client
from fugue import transform
from fugue_dask import DaskExecutionEngine
from datasetsforecast.losses import mse, mae, smape

evaluate 函数接收一个时间序列和一个窗口的唯一组合，并为 df 中的每个模型计算不同的 metrics。

def evaluate(df: pd.DataFrame, metrics: List[Callable]) -> pd.DataFrame:
    eval_ = {}
    models = df.loc[:, ~df.columns.str.contains('unique_id|y|ds|cutoff|lo|hi')].columns
    for model in models:
        eval_[model] = {}
        for metric in metrics:
            eval_[model][metric.__name__] = metric(df['y'], df[model])
    eval_df = pd.DataFrame(eval_).rename_axis('metric').reset_index()
    eval_df.insert(0, 'cutoff', df['cutoff'].iloc[0])
    eval_df.insert(0, 'unique_id', df['unique_id'].iloc[0])
    return eval_df
str_models = cv_df.loc[:, ~cv_df.columns.str.contains('unique_id|y|ds|cutoff|lo|hi')].columns
str_models = ','.join([f"{model}:float" for model in str_models])
cv_df['cutoff'] = cv_df['cutoff'].astype(str)
cv_df['unique_id'] = cv_df['unique_id'].astype(str)

让我们创建一个 dask 客户端。

client = Client() # 没有这个，Dask 不在分布式模式下运行。
# fugue.dask.dataframe.default.partitions 决定了新 DaskDataFrame 的默认分区数。
engine = DaskExecutionEngine({"fugue.dask.dataframe.default.partitions": 96})

transform 函数接受 evaluate 函数，并将其应用于每个时间序列 (unique_id) 和交叉验证窗口 (cutoff) 的组合，使用我们之前创建的 dask 客户端。

evaluation_df = transform(
    cv_df.loc[:, ~cv_df.columns.str.contains('lo|hi')], 
    evaluate, 
    engine="dask",
    params={'metrics': [mse, mae, smape]}, 
    schema=f"unique_id:str,cutoff:str,metric:str, {str_models}", 
    as_local=True,
    partition={'by': ['unique_id', 'cutoff']}
)
/home/ubuntu/miniconda/envs/statsforecast/lib/python3.10/site-packages/distributed/client.py:3109: UserWarning:

Sending large graph of size 49.63 MiB.
This may cause some slowdown.
Consider scattering data ahead of time and using futures.

evaluation_df.head()
unique_id cutoff metric SeasonalNaive Naive HistoricAverage CrostonOptimized ADIDA IMAPA AutoETS AutoNHITS AutoTFT LGBMRegressor XGBRegressor LinearRegression
0 FOODS_3_003_WI_3 2016-02-28 mse 1.142857 1.142857 0.816646 0.816471 1.142857 1.142857 1.142857 1.142857 1.142857 0.832010 1.020361 0.887121
1 FOODS_3_003_WI_3 2016-02-28 mae 0.571429 0.571429 0.729592 0.731261 0.571429 0.571429 0.571429 0.571429 0.571429 0.772788 0.619949 0.685413
2 FOODS_3_003_WI_3 2016-02-28 smape 71.428574 71.428574 158.813507 158.516235 200.000000 200.000000 200.000000 71.428574 71.428574 145.901947 188.159164 178.883743
3 FOODS_3_013_CA_3 2016-04-24 mse 4.000000 6.214286 2.406764 3.561202 2.267853 2.267600 2.268677 2.750000 2.125000 2.160508 2.370228 2.289606
4 FOODS_3_013_CA_3 2016-04-24 mae 1.500000 2.142857 1.214286 1.340446 1.214286 1.214286 1.214286 1.107143 1.142857 1.140084 1.157548 1.148813 
# 计算每个交叉验证窗口的平均指标
evaluation_df.groupby(['cutoff', 'metric']).mean(numeric_only=True)
SeasonalNaive Naive HistoricAverage CrostonOptimized ADIDA IMAPA AutoETS AutoNHITS AutoTFT LGBMRegressor XGBRegressor LinearRegression
cutoff metric
2016-02-28 mae 1.744289 2.040496 1.730704 1.633017 1.527965 1.528772 1.497553 1.434938 1.485419 1.688403 1.514102 1.576320
mse 14.510710 19.080585 12.858994 11.785032 11.114497 11.100909 10.347847 10.010982 10.964664 10.436206 10.968788 10.792831
smape 85.202042 87.719086 125.418488 124.749908 127.591858 127.704102 127.790672 79.132614 80.983368 118.489983 140.420578 127.043137
2016-03-27 mae 1.795973 2.106449 1.754029 1.662087 1.570701 1.572741 1.535301 1.432412 1.502393 1.712493 1.600193 1.601612
mse 14.810259 26.044472 12.804104 12.020620 12.083861 12.120033 11.315013 9.445867 10.762877 10.723589 12.924312 10.943772
smape 87.407471 89.453247 123.587196 123.460030 123.428459 123.538521 123.612991 79.926781 82.013168 116.089699 138.885941 127.304871
2016-04-24 mae 1.785983 1.990774 1.762506 1.609268 1.527627 1.529721 1.501820 1.447401 1.505127 1.692946 1.541845 1.590985
mse 13.476350 16.234917 13.151311 10.647048 10.072225 10.062395 9.393439 9.363891 10.436214 10.347073 10.774202 10.608137
smape 89.238815 90.685867 121.124947 119.721245 120.325401 120.345284 120.649582 81.402748 83.614029 113.334198 136.755234 124.618622

结果在之前的实验中显示。

模型 均方误差 (MSE)
MQCNN 10.09
DeepAR-student_t 10.11
DeepAR-lognormal 30.20
DeepAR 9.13
NPTS 11.53

顶级模型：DeepAR、AutoNHITS、AutoETS。

错误分布

%%capture
!pip install seaborn
import matplotlib.pyplot as plt
import seaborn as sns
evaluation_df_melted = pd.melt(evaluation_df, id_vars=['unique_id', 'cutoff', 'metric'], var_name='model', value_name='error')

SMAPE（对称平均绝对百分比误差）

sns.violinplot(evaluation_df_melted.query('metric=="smape"'), x='error', y='model')

为系列组选择模型

特征：

  • 具有所有不同模型预测的统一数据框
  • 简单的集成
  • 例如：平均预测
  • 或 MinMax（选择即是集成）
# 为每个时间序列、指标和交叉验证窗口选择最佳模型
evaluation_df['best_model'] = evaluation_df.idxmin(axis=1, numeric_only=True)
# 统计每个模型在每个指标和交叉验证窗口中获胜的次数
count_best_model = evaluation_df.groupby(['cutoff', 'metric', 'best_model']).size().rename('n').to_frame().reset_index()
# 绘制结果
sns.barplot(count_best_model, x='n', y='best_model', hue='metric')

统一中的多样性：一种包容性的预测饼图。

# 对于均方误差（mse），计算模型获胜的次数
eval_series_df = evaluation_df.query('metric == "mse"').groupby(['unique_id']).mean(numeric_only=True)
eval_series_df['best_model'] = eval_series_df.idxmin(axis=1)
counts_series = eval_series_df.value_counts('best_model')
plt.pie(counts_series, labels=counts_series.index, autopct='%.0f%%')
plt.show()

sf.plot(Y_df, cv_df.drop(columns=['cutoff', 'y']), 
        max_insample_length=28 * 6, 
        models=['AutoNHITS'],
        unique_ids=eval_series_df.query('best_model == "AutoNHITS"').index[:8])

为不同组系列选择预测方法

# 合并每个时间序列数据框的最佳模型
# 并根据该数据框过滤预测结果
# 对于每个时间序列
fcst_df = pd.melt(fcst_df.set_index('unique_id'), id_vars=['ds'], var_name='model', value_name='forecast', ignore_index=False)
fcst_df = fcst_df.join(eval_series_df[['best_model']])
fcst_df[['model', 'pred-interval']] = fcst_df['model'].str.split('-', expand=True, n=1)
fcst_df = fcst_df.query('model == best_model')
fcst_df['name'] = [f'forecast-{x}' if x is not None else 'forecast' for x in fcst_df['pred-interval']]
fcst_df = pd.pivot_table(fcst_df, index=['unique_id', 'ds'], values=['forecast'], columns=['name']).droplevel(0, axis=1).reset_index()
sf.plot(Y_df, fcst_df, max_insample_length=28 * 3)

技术债务

  • 在完整数据集中训练统计模型。
  • 增加神经自编码模型中的 num_samples 数量。
  • 包括其他模型，如 Theta、ARIMA、RNN、LSTM 等。

进一步的材料

Give us a ⭐ on Github