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

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

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

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

在时间序列预测与分析的领域中，识别出最适合特定系列组的模型是一项相对复杂的任务。通常，这一选择过程在很大程度上依赖于直觉，而这些直觉可能并不一定与我们数据集的实证现实相符。

在本教程中，我们旨在为不同系列组的模型选择提供一种更加结构化、以数据驱动的方法，基于M5基准数据集。这个在预测领域广为人知的数据集使我们能够展示我们方法论的多样性和力量。

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

StatsForecast

• 基线模型：这些模型简单但通常在提供对预测问题的初步视角时非常有效。我们将使用SeasonalNaive和HistoricAverage模型来代表这一类别。
• 间歇性模型：针对具有零星、不连续需求的系列，我们将采用如CrostonOptimized、IMAPA和ADIDA等模型。这些模型特别适用于处理零膨胀系列。
• 状态空间模型：这些是利用系统的数学描述进行预测的统计模型。来自statsforecast库的AutoETS模型属于这一类别。

MLForecast

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

NeuralForecast

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

使用Nixtla的库套件，我们将能够通过数据推动我们的模型选择过程，确保利用最适合的数据集中各特定系列组的模型。

大纲：

• 数据读取：在这个初始步骤中，我们将数据集加载到内存中，使其可供后续的分析和预测使用。在这个阶段，理解数据集的结构和细微差别是非常重要的。

• 使用统计和深度学习方法进行预测：我们应用从基本统计技术到先进深度学习模型的多种预测方法。目标是基于我们的数据集生成未来28天的预测。

• 在不同窗口上评估模型性能：我们评估模型在不同窗口上的表现。

• 为一组时间序列选择最佳模型：通过性能评估，我们确定每组时间序列的最佳模型。此步骤确保所选择的模型适合每组的独特特性。

• 筛选最佳可能预测：最后，我们筛选所选模型生成的预测，以获得最有前景的预测。这是我们的最终输出，并代表根据我们的模型为每个序列提供的最佳可能预测。

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 类中的绘图方法绘制一些序列。此方法会从数据集中随机打印 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、DilatedRNN、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)

使用plot函数来探索模型和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（整数）：预测范围，以我们希望预测的未来步数表示。例如，如果我们正在预测每小时的数据，h=24将表示24小时的预测。
• step_size（整数）：每个交叉验证窗口之间的步长。该参数决定了我们希望多频繁地运行预测过程。
• n_windows（整数）：用于交叉验证的窗口数量。该参数定义了我们希望评估的过去预测过程的数量。

这些参数允许我们控制交叉验证过程的范围和粒度。通过调整这些设置，我们可以在计算成本和交叉验证的彻底性之间找到平衡。

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 （int）：代表向未来 h 步的预测。在这种情况下，表示提前 24 小时。
• step_size （int）：每个窗口之间的步长。换句话说：您希望多频繁地运行预测过程。
• n_windows （int）：用于交叉验证的窗口数量。换句话说：您希望评估过去的多少次预测过程。
• 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

前3个模型：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')

选择系列组的模型

特性：

• 一个统一的数据框，包含所有不同模型的预测
• 简单的集成
• 例如：平均预测
• 或最小值最大值（选择就是集成）

# 为每个时间序列、指标和交叉验证窗口选择最佳模型
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等。

进一步的材料

• 可用模型 StatsForecast
• 可用模型 NeuralForecast
• 缩放器和损失函数
• 开始使用 NeuralForecast
• 层次性调和
• 分布式机器学习预测 (树)
• 使用 StatsForecast 训练数百万个时间序列
• 在 Databricks 上使用 Nixtla 进行间歇性需求预测

Give us a ⭐ on Github