ADIDA模型


# 此单元格将不会被渲染,但用于隐藏警告并限制显示的行数。

import warnings
warnings.filterwarnings("ignore")

import logging
logging.getLogger('statsforecast').setLevel(logging.ERROR)

import pandas as pd
pd.set_option('display.max_rows', 6)

使用ADIDA模型Statsforecast的逐步指南。

在本指南中,我们将熟悉主要的StatsForecast类和一些相关的方法,如StatsForecast.plotStatsForecast.forecastStatsForecast.cross_validation

目录

介绍

聚合-拆分间歇性需求方法(ADIDA)是一种用于预测表现出间歇性需求模式的产品需求的预测方法。间歇性需求模式的特点是观察到大量零值,这可能使得预测变得具有挑战性。

ADIDA方法通过时间聚合来减少零观察值的数量,并减轻在时间间隔中观察到的方差影响。该方法使用等大小的时间桶进行不重叠的时间聚合,并预测在预先指定的交货时间内的需求。时间桶被设定为平均需求间隔时间,即两个连续非零观察之间的平均时间。

该方法采用简单指数平滑(SES)技术来获得预测。SES是一种常用的时间序列预测技术,由于其简单性和有效性,通常用于生成准确的预测。

ADIDA方法有几个优点。它易于实施,并且可以用于各种间歇性需求模式。该方法还提供准确的预测,可以用来预测在预先指定的交货时间内的需求。

然而,ADIDA方法也存在一些局限性。该方法假设时间桶大小相等,但并不是所有间歇性需求模式都符合这一假设。此外,该方法可能不适用于具有复杂模式或趋势的时间序列数据。

总体而言,ADIDA方法是适用于间歇性需求模式的有用预测技术,可以帮助减轻零观察值的影响,并生成准确的需求预测。

ADIDA 模型

什么是间歇性需求?

间歇性需求是一种需求模式,其特点是在事件或销售的发生上表现出不规律和偶发性。换句话说,它指的是某种产品或服务的需求在一段时间内间歇性地出现,在这些时间段内,销售或重大事件不会发生。

间歇性需求不同于恒定或规律需求,在后者中,销售在时间上呈现可预测且一致的方式。相反,在间歇性需求中,没有销售的时间段可能会很长,并且事件的发生并不规律。

这种类型的需求可以在不同的行业和情境中出现,例如低消费产品、季节性产品、高变异产品、生命周期短的产品,或在某些需求取决于特定事件或外部因素的情况下。

间歇性需求在预测和库存管理中可能带来挑战,因为很难预测销售何时会发生以及数量。像我之前提到的克罗斯顿模型等方法被用来应对间歇性需求,并为这种类型的需求模式生成更准确和适当的预测。

间歇性需求的问题

间歇性需求在库存管理和需求预测中可能会呈现出多种挑战和问题。一些与间歇性需求相关的常见问题如下:

  1. 不可预测的变异性:间歇性需求可能具有不可预测的变异性,从而使规划和预测变得困难。需求模式可能不规律,并且销售期与非销售期之间的波动可能非常剧烈。

  2. 销售频率低:间歇性需求的特点是长时间没有销售。这可能导致库存管理的困难,因为必须保持足够的库存以满足发生时的需求,同时避免在非销售期积压过剩库存。

  3. 预测误差:预测间歇性需求比预测恒定需求更难以确定。传统的预测模型可能无法充分捕捉间歇性需求中的变异性和缺乏模式,这可能导致对未来需求的估计出现显著误差。

  4. 对供应链的影响:间歇性需求可能影响供应链的效率,并在生产计划、供应商管理和物流方面造成困难。交货时间和库存水平必须调整以满足不可预测的需求。

  5. 运营成本:在间歇性需求的情况下管理库存可能会增加运营成本。维持在非销售期的适当库存和管理库存水平可能需要在存储和物流方面增加额外投入。

为了解决这些问题,采用了针对间歇性需求管理的特定方法,例如专门的预测模型、产品分类技术和定制的库存策略。这些解决方案旨在最小化间歇性需求中变异性和缺乏模式的影响,优化库存管理,改善供应链效率。

ADIDA 模型

ADIDA 模型基于简单指数平滑 (SES) 方法,并使用时间聚合来处理间歇性需求的问题。该模型的数学推导可以总结如下:

\(S_t\) 为时间 \(t\) 的需求,其中 \(t = 1, 2, ..., T\)。平均间需求间隔用 MI 表示,MI 是两个连续非零需求之间的平均时间。时间桶的大小设置为 MI。

然后将需求数据聚合到大小为 MI 的不重叠时间桶中。设 \(B_t\) 为桶 \(t\) 的需求,其中 \(t = 1, 2, ..., T/MI\)。聚合后的需求数据可以表示为:

\[B_t = \sum S_t, \text{ for } (t-1)*MI + 1 ≤ j ≤ t*MI\]

然后将 SES 方法应用于聚合的需求数据,以获得预测。桶 \(t\) 的预测用 \(F_t\) 表示。SES 方法涉及使用以下公式基于时间 \(t\) 的实际需求 \(D_t\) 和前一个时间段的估计水平 \(L_{t-1}\) 来估计水平 \(L_t\)

\[L_t = \alpha * D_t + (1 - α) * L_{t-1}\]

其中 \(\alpha\) 是控制当前需求值权重的平滑参数。

\(t\) 的预测随后通过使用前一个时间段的估计水平 \(L_{t-1}\) 获得,如下所示:

\[F_t = L_{t-1}\]

然后将预测进行解聚,以获得原始时间段的需求预测。设 \(Y_t\) 为时间 \(t\) 的需求预测。解聚可以使用以下公式进行:

\[Y_t = F_t / MI, \text{ for } (t-1)*MI + 1 ≤ j ≤ t*MI\]

如何确定ADIDA模型是否适合特定数据集?

要确定ADIDA模型是否适合特定数据集,可以遵循以下步骤:

  1. 分析需求模式:检查数据的需求模式,以确定它是否符合间歇性模式。间歇性数据的特点是在特定时期内有很高比例的零需求和偶发需求。

  2. 评估季节性:检查数据中是否存在明显的季节性。ADIDA模型假设数据没有季节性,或者可以通过时间聚合来处理。如果数据表现出复杂的季节性或无法通过时间聚合来处理,则ADIDA模型可能不适合。

  3. 数据要求:考虑ADIDA模型的数据要求。该模型需要历史需求数据,并能够计算非零需求之间的平均间隔。确保您有足够的数据以估计参数,并且数据的频率适合时间聚合。

  4. 性能评估:对ADIDA模型在特定数据集上的表现进行性能评估。将模型生成的预测与实际需求值进行比较,并使用均方误差(MAE)或均方根误差(MSE)等评估指标。如果模型在数据集上表现良好并产生准确的预测,这表明该模型适合该数据集。

  5. 与其他模型比较:将ADIDA模型的性能与适合间歇性数据的其他预测模型进行比较。考虑像Croston、Syntetos-Boylan近似(SBA)或专门为间歇性数据开发的基于指数平滑技术的模型。如果ADIDA模型的表现与其他模型相似或更好,则可以认为它适合。

请记住,ADIDA模型的适用性可能取决于数据的具体性质和预测问题的背景。建议进行全面分析,并尝试不同的模型,以确定最适合该数据集的方法。

加载库和数据

Tip

需要使用Statsforecast。有关安装,请参见 说明

接下来,我们导入绘图库并配置绘图样式。

import matplotlib.pyplot as plt
import seaborn as sns
from statsmodels.graphics.tsaplots import plot_acf
from statsmodels.graphics.tsaplots import plot_pacf
import plotly.graph_objects as go
plt.style.use('grayscale') # 五三八 灰度 经典
plt.rcParams['lines.linewidth'] = 1.5
dark_style = {
    'figure.facecolor': '#008080',  # #212946
    'axes.facecolor': '#008080',
    'savefig.facecolor': '#008080',
    'axes.grid': True,
    'axes.grid.which': 'both',
    'axes.spines.left': False,
    'axes.spines.right': False,
    'axes.spines.top': False,
    'axes.spines.bottom': False,
    'grid.color': '#000000',  #2A3459
    'grid.linewidth': '1',
    'text.color': '0.9',
    'axes.labelcolor': '0.9',
    'xtick.color': '0.9',
    'ytick.color': '0.9',
    'font.size': 12 }
plt.rcParams.update(dark_style)

from pylab import rcParams
rcParams['figure.figsize'] = (18,7)
import pandas as pd

df = pd.read_csv("https://raw.githubusercontent.com/Naren8520/Serie-de-tiempo-con-Machine-Learning/main/Data/tipos_malarias_choco_colombia.csv", sep=";", usecols=[0,4])
df = df.dropna()
df.head()
semanas malaria_falciparum
0 2007-12-31 50.0
1 2008-01-07 62.0
2 2008-01-14 76.0
3 2008-01-21 64.0
4 2008-01-28 38.0

StatsForecast的输入始终是一个长格式的数据框,包含三列:unique_id、ds 和 y:

  • unique_id(字符串、整数或类别)代表系列的标识符。

  • ds(日期戳)列应为Pandas期望的格式,理想情况下应为YYYY-MM-DD的日期格式或YYYY-MM-DD HH:MM:SS的时间戳格式。

  • y(数值)代表我们希望预测的测量值。

df["unique_id"]="1"
df.columns=["ds", "y", "unique_id"]
df.head()
ds y unique_id
0 2007-12-31 50.0 1
1 2008-01-07 62.0 1
2 2008-01-14 76.0 1
3 2008-01-21 64.0 1
4 2008-01-28 38.0 1
print(df.dtypes)
ds            object
y            float64
unique_id     object
dtype: object

我们需要将object类型转换为日期时间和数字类型。

df["ds"] = pd.to_datetime(df["ds"])
df["y"] = df["y"].astype(float).astype("int64")

使用 plot 方法探索数据

使用 StatsForecast 类中的 plot 方法绘制系列。该方法从数据集中打印随机系列,对于基本的探索性数据分析非常有用。

from statsforecast import StatsForecast

StatsForecast.plot(df)

自相关图

fig, axs = plt.subplots(nrows=1, ncols=2)

plot_acf(df["y"],  lags=30, ax=axs[0],color="fuchsia")
axs[0].set_title("Autocorrelation");

plot_pacf(df["y"],  lags=30, ax=axs[1],color="lime")
axs[1].set_title('Partial Autocorrelation')

plt.show();

时间序列的分解

如何分解时间序列以及为什么要这样做?

在时间序列分析中,预测新值时,了解过去的数据非常重要。更正式地说,了解值随时间变化的模式是非常重要的。我们的预测值可能朝错误的方向波动的原因有很多。基本上,时间序列由四个组件组成。这些组件的变化导致时间序列模式的改变。这些组件是:

  • 水平(Level): 这是随着时间平均的主要值。
  • 趋势(Trend): 趋势是导致时间序列中增加或减少模式的值。
  • 季节性(Seasonality): 这是时间序列中短时间内发生的循环事件,导致时间序列中短期增加或减少的模式。
  • 残差/噪声(Residual/Noise): 这些是时间序列中的随机变化。

将这些组件随着时间结合在一起,形成了一个时间序列。大多数时间序列由水平和噪声/残差组成,趋势或季节性是可选的值。

如果季节性和趋势是时间序列的一部分,那么将会对预测值产生影响。因为预测的时间序列模式可能与之前的时间序列不同。

时间序列中组件的组合可以分为两种类型: * 加法型(Additive) * 乘法型(Multiplicative)

加法型时间序列

如果时间序列的组件是相加而形成该时间序列的,则该时间序列被称为加法型时间序列。通过可视化,我们可以说,如果时间序列的增加或减少模式在整个序列中是相似的,则该时间序列是加法型的。任何加法型时间序列的数学函数可以表示为: \[y(t) = level + Trend + seasonality + noise\]

乘法型时间序列

如果时间序列的组件是相乘在一起的,则该时间序列被称为乘法型时间序列。通过可视化,如果时间序列随时间呈指数增长或下降,则该时间序列可以视为乘法型时间序列。乘法型时间序列的数学函数可以表示为: \[y(t) = Level * Trend * seasonality * Noise\]

from statsmodels.tsa.seasonal import seasonal_decompose
from plotly.subplots import make_subplots
import plotly.graph_objects as go

def plot_seasonal_decompose(
    x,
    model='additive',
    filt=None,
    period=None,
    two_sided=True,
    extrapolate_trend=0,
    title="Seasonal Decomposition"):

    result = seasonal_decompose(
            x, model=model, filt=filt, period=period,
            two_sided=two_sided, extrapolate_trend=extrapolate_trend)
    fig = make_subplots(
            rows=4, cols=1,
            subplot_titles=["Observed", "Trend", "Seasonal", "Residuals"])
    for idx, col in enumerate(['observed', 'trend', 'seasonal', 'resid']):
        fig.add_trace(
            go.Scatter(x=result.observed.index, y=getattr(result, col), mode='lines'),
                row=idx+1, col=1,
            )
    return fig
plot_seasonal_decompose(
    df["y"],
    model="additive",
    period=52,
    title="Seasonal Decomposition")

将数据分为训练集和测试集

让我们将数据分为两部分 1. 用于训练我们的ADIDA模型的数据。 2. 用于测试我们模型的数据。

对于测试数据,我们将使用最后25周的数据来测试和评估我们模型的性能。

train = df[df.ds<='2022-07-04'] 
test = df[df.ds>'2022-07-04'] 
train.shape, test.shape
((758, 3), (25, 3))

现在让我们绘制训练数据和测试数据。

sns.lineplot(train,x="ds", y="y", label="Train", linestyle="--",linewidth=2)
sns.lineplot(test, x="ds", y="y", label="Test", linewidth=2, color="yellow")
plt.title("Falciparum Malaria");
plt.show()

ADIDA Model 的实现与 StatsForecast

要了解更多关于 ADIDA Model 函数的参数,它们列在下面。有关更多信息,请访问 文档

alias : str
    模型的自定义名称。
prediction_intervals : Optional[ConformalIntervals]
    用于计算符合预测区间的信息。
    默认情况下,模型将计算原生预测区间。

加载库

from statsforecast import StatsForecast
from statsforecast.models import ADIDA

实例化模型

导入并实例化模型。设置参数有时比较棘手。Rob Hyndmann大师的文章季节周期对于season_length可能很有帮助。

season_length = 52 # 每小时数据 
horizon = len(test) # 预测数量

# 我们将要使用的模型称为
models = [ADIDA()]

我们通过实例化一个新的 StatsForecast 对象并设定以下参数来拟合模型:

模型:模型的列表。选择您想要的模型并导入它们。

  • freq: 一个字符串,表示数据的频率。(请参见 pandas 可用的频率.)

  • n_jobs: n_jobs: 整数,表示在并行处理时使用的作业数量,使用 -1 表示使用所有核心。

  • fallback_model: 如果某个模型失败,则使用的模型。

所有设置都会传递给构造函数。然后您调用其 fit 方法并传入历史数据框。

sf = StatsForecast(df=df,
                   models=models,
                   freq='W', 
                   n_jobs=-1)

拟合模型

在这里,我们调用 fit() 方法来拟合模型。

sf.fit()
StatsForecast(models=[ADIDA])

让我们看看我们的 ADIDA 模型 的结果。我们可以通过以下指令进行观察:

result=sf.fitted_[0,0].model_
result
{'mean': array([132.52887], dtype=float32)}

预测方法

如果您希望在具有多个系列或模型的生产环境中提高速度,我们建议使用 StatsForecast.forecast 方法,而不是 .fit.predict

主要区别在于 forecast() 方法不存储拟合值,并且在分布式环境中具有很高的可扩展性。

预测方法接受两个参数:预测下一个 h(时间范围)和 level

  • h (int): 代表预测未来 h 步。在本例中,表示提前 25 周。

此处的预测对象是一个新的数据框,其中包含一列模型名称和 y hat 值,以及不确定性区间的列。根据您的计算机情况,该步骤通常需要约 1 分钟。

Y_hat = sf.forecast(horizon)
Y_hat
ds ADIDA
unique_id
1 2023-01-01 132.52887
1 2023-01-08 132.52887
1 2023-01-15 132.52887
... ... ...
1 2023-06-04 132.52887
1 2023-06-11 132.52887
1 2023-06-18 132.52887

25 rows × 2 columns

Y_hat=Y_hat.reset_index()
Y_hat
unique_id ds ADIDA
0 1 2023-01-01 132.52887
1 1 2023-01-08 132.52887
2 1 2023-01-15 132.52887
... ... ... ...
22 1 2023-06-04 132.52887
23 1 2023-06-11 132.52887
24 1 2023-06-18 132.52887

25 rows × 3 columns

将预测结果与真实值合并

Y_hat1 = pd.concat([df,Y_hat],  keys=['unique_id', 'ds'])
Y_hat1
ds y unique_id ADIDA
unique_id 0 2007-12-31 50.0 1 NaN
1 2008-01-07 62.0 1 NaN
2 2008-01-14 76.0 1 NaN
... ... ... ... ... ...
ds 22 2023-06-04 NaN 1 132.52887
23 2023-06-11 NaN 1 132.52887
24 2023-06-18 NaN 1 132.52887

808 rows × 4 columns

fig, ax = plt.subplots(1, 1)
plot_df = pd.concat([df, Y_hat1]).set_index('ds').tail(300)
plot_df['y'].plot(ax=ax, linewidth=2)
plot_df[ "ADIDA"].plot(ax=ax, linewidth=2, color="yellow")
ax.set_title(' Forecast', fontsize=22)
ax.set_ylabel('Falciparum Malaria (Weekly data)', fontsize=20)
ax.set_xlabel('Monthly [t]', fontsize=20)
ax.legend(prop={'size': 15})
ax.grid(True)

带置信区间的预测方法

要生成预测,请使用预测方法。

预测方法接受两个参数:预测下一个 h(代表时间范围)和 level

  • h (int): 表示预测向前的 h 步。在这种情况下,向前预测 25 周。

此处的预测对象是一个新的数据框,包含一个模型名称的列和 y 预测值,以及不确定性区间的列。

此步骤应在 1 秒钟内完成。

forecast_df = sf.predict(h=horizon) 

forecast_df.head()
ds ADIDA
unique_id
1 2023-01-01 132.52887
1 2023-01-08 132.52887
1 2023-01-15 132.52887
1 2023-01-22 132.52887
1 2023-01-29 132.52887

我们可以使用pandas函数pd.concat()将预测结果与历史数据连接,然后可以使用这个结果进行绘图。

pd.concat([df, forecast_df]).set_index('ds')
y unique_id ADIDA
ds
2007-12-31 50.0 1 NaN
2008-01-07 62.0 1 NaN
2008-01-14 76.0 1 NaN
... ... ... ...
2023-06-04 NaN NaN 132.52887
2023-06-11 NaN NaN 132.52887
2023-06-18 NaN NaN 132.52887

808 rows × 3 columns

df_plot= pd.concat([df, forecast_df]).set_index('ds').tail(200)
df_plot
y unique_id ADIDA
ds
2019-08-26 191.0 1 NaN
2019-09-02 205.0 1 NaN
2019-09-09 200.0 1 NaN
... ... ... ...
2023-06-04 NaN NaN 132.52887
2023-06-11 NaN NaN 132.52887
2023-06-18 NaN NaN 132.52887

200 rows × 3 columns

现在让我们可视化我们的预测结果以及时间序列的历史数据。

plt.plot(df_plot['y'],label="Actual", linewidth=2.5)
plt.plot(df_plot['ADIDA'], label="ADIDA", color="yellow") # '-', '--', '-.', ':',

plt.title("Falciparum Malaria (Weekly data)");
plt.xlabel("Weekly")
plt.ylabel("")
Text(0, 0.5, '')

让我们使用Statsforecast中的plot函数绘制相同的图表,如下所示。

sf.plot(df, forecast_df)

交叉验证

在之前的步骤中,我们使用历史数据来预测未来。然而,为了评估其准确性,我们还希望了解模型在过去的表现。为了评估您在数据上的模型的准确性和稳健性,请执行交叉验证。

对于时间序列数据,交叉验证是通过在历史数据中定义一个滑动窗口并预测其后的时期来完成的。这种形式的交叉验证使我们能够更好地估计模型在更广泛的时间实例中的预测能力,同时保持训练集中的数据是连续的,这也是我们的模型所要求的。

以下图表描述了这样的交叉验证策略:

img

执行时间序列交叉验证

时间序列模型的交叉验证被认为是一种最佳实践,但大多数实现都非常慢。statsforecast库将交叉验证实现为分布式操作,使得执行过程更加省时。如果您有大型数据集,还可以使用Ray、Dask或Spark在分布式集群中执行交叉验证。

在本例中,我们希望评估每个模型在过去5个月的表现(n_windows=),每隔两个月进行一次预测(step_size=12)。根据您的计算机,此步骤大约需要1分钟。

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

  • df: 训练数据框

  • h (int): 代表预测的未来h步。在本例中,为12个月。

  • step_size (int): 每个窗口之间的步长。换句话说:您希望多长时间运行一次预测过程。

  • n_windows(int): 用于交叉验证的窗口数量。换句话说:您希望评估过去多少个预测过程。

crossvalidation_df = sf.cross_validation(df=df,
                                         h=horizon,
                                         step_size=30,
                                         n_windows=5)

crossvaldation_df 对象是一个新的数据框,包括以下列:

  • unique_id: 索引。如果你不喜欢使用索引,只需运行 crossvalidation_df.resetindex()
  • ds: 日期戳或时间索引
  • cutoff: n_windows 的最后一个日期戳或时间索引。
  • y: 真实值
  • model: 包含模型名称和拟合值的列。
crossvalidation_df
ds cutoff y ADIDA
unique_id
1 2020-03-22 2020-03-15 317.0 251.901505
1 2020-03-29 2020-03-15 332.0 251.901505
1 2020-04-05 2020-03-15 306.0 251.901505
... ... ... ... ...
1 2022-12-11 2022-07-03 151.0 336.747375
1 2022-12-18 2022-07-03 97.0 336.747375
1 2022-12-25 2022-07-03 42.0 336.747375

125 rows × 4 columns

模型评估

我们现在可以使用适当的准确率指标来计算预测的准确性。在这里,我们将使用均方根误差(RMSE)。为此,我们首先需要安装datasetsforecast,这是一个由Nixtla开发的Python库,其中包含计算RMSE的函数。

%%capture
!pip install datasetsforecast
from datasetsforecast.losses import rmse

计算RMSE的函数需要两个参数:

  1. 实际值。
  2. 预测值,在本例中为ADIDA模型
rmse = rmse(crossvalidation_df['y'], crossvalidation_df["ADIDA"])
print("RMSE using cross-validation: ", rmse)
RMSE using cross-validation:  84.45186

致谢

我们要感谢 Naren Castellon 为撰写本教程所做的努力。

参考文献

  1. Changquan Huang • Alla Petukhina. Springer系列 (2022). 使用Python的应用时间序列分析与预测。
  2. Ivan Svetunkov. 增强动态自适应模型(ADAM)的预测与分析
  3. James D. Hamilton. 时间序列分析 普林斯顿大学出版社, 新泽西州普林斯顿, 第1版, 1994.
  4. Nixtla 参数
  5. Pandas可用频率
  6. Rob J. Hyndman和George Athanasopoulos (2018). “预测原则与实践,时间序列交叉验证”。
  7. 季节性周期- Rob J Hyndman

Give us a ⭐ on Github