Croston优化模型


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

import warnings
warnings.filterwarnings("ignore")

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

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

使用Statsforecast的Croston优化模型的逐步指南。

在本教程中，我们将熟悉主要的StatsForecast类及一些相关方法，如StatsForecast.plot、StatsForecast.forecast和StatsForecast.cross_validation等。

让我们开始吧！！！

目录

介绍

Croston优化模型是一种为间歇性需求时间序列数据设计的预测方法。它是Croston方法的扩展，该方法最初是为预测偶发需求模式而开发的。

间歇性需求时间序列的特点是非零需求值的不规则和偶发出现，通常伴随着长时间的零需求。传统的预测方法可能难以有效处理这种模式。

Croston优化模型通过结合两个关键组成部分来应对这一挑战：指数平滑和间歇性需求估计。

  1. 指数平滑：Croston优化模型使用指数平滑来捕获间歇性需求数据中的趋势和季节性。这有助于识别潜在模式并进行更准确的预测。

  2. 间歇性需求估计：由于间歇性需求数据通常包括长时间的零需求，Croston优化模型采用单独的估计过程来评估非零需求值的出现和大小。它估计非零需求间隔的发生概率和平均大小，从而能够更好地预测间歇性需求。

Croston优化模型旨在在过度预测和不足预测间歇性需求之间找到平衡，这在传统预测方法中是常见的挑战。通过明确建模间歇性需求模式，它能够为间歇性需求时间序列数据提供更准确的预测。

值得注意的是，Croston优化模型有多种变体和适应性，根据特定的预测场景进行了不同的修改和增强。这些变体可能会结合额外的特性或算法，以进一步提高预测的准确性。

Croston 优化方法

Croston 优化模型可以数学上定义如下：

  1. 初始化：
    • 设 \((y_t)\) 代表时间 \(t\) 的间歇性需求时间序列数据。
    • 初始化两个变量集：\((p_t)\) 用于事件发生的概率，以及 \((q_t)\) 用于非零需求间隔的平均大小。
    • 初始化预测 \((F_t)\) 和预测误差 \((E_t)\) 变量为零。
  2. 计算 \((p_t)\) 和 \((q_t)\)：
    • 使用指数平滑法计算间歇性需求发生概率 \((p_t)\)： \[[p_t = \alpha + (1 - \alpha)(p_{t-1}),]\] 其中 \((\alpha)\) 是平滑参数（通常设定在 0.1 到 0.3 之间）。
    • 使用指数平滑法计算非零需求间隔的平均大小 \((q_t)\)： \[[q_t = \beta \cdot y_t + (1 - \beta)(q_{t-1}),]\] 其中 \((\beta)\) 是平滑参数（通常设定在 0.1 到 0.3 之间）。
  3. 预测：
    • 如果 \((y_t > 0)\) （非零需求发生）：
      • 将预测 \((F_t)\) 计算为前一个预测 \((F_{t-1})\) 除以非零需求间隔的平均大小 \((q_{t-1})\)： \[[F_t = \frac{{F_{t-1}}}{{q_{t-1}}}]\]
      • 将预测误差 \((E_t)\) 计算为实际需求 \((y_t)\) 与预测 \((F_t)\) 之间的差值： \[[E_t = y_t - F_t]\]
    • 如果 \((y_t = 0)\) （零需求发生）：
      • 将预测 \((F_t)\) 和预测误差 \((E_t)\) 设定为零。
  4. 更新模型：
    • 使用步骤 2 中描述的指数平滑法更新间歇性需求发生概率 \((p_t)\) 和非零需求间隔的平均大小 \((q_t)\)。
  5. 针对时间序列中的每个时间点重复步骤 3 和 4。

Croston 优化模型利用指数平滑法捕捉间歇性需求数据中的趋势和季节性，并单独估计事件发生的概率和非零需求间隔的平均大小，以有效处理间歇性需求模式。通过根据观察数据更新模型参数，它为间歇性需求时间序列提供了更好的预测。

优化克罗斯顿模型的一些特性

优化克罗斯顿模型是经典克罗斯顿模型的一个改进，旨在预测间歇性需求。经典克罗斯顿模型使用历史订单的加权平均和订单之间的平均间隔来预测需求。优化克罗斯顿模型则使用概率函数来预测订单之间的平均间隔。

研究表明，优化克罗斯顿模型在处理不规则需求的时间序列时，比经典克罗斯顿模型更为准确。优化克罗斯顿模型对于不同类型的间歇性时间序列也更具适应性。

优化克罗斯顿模型具有以下特性：

  • 对于不规则需求的时间序列也能够精准。
  • 适应不同类型的间歇性时间序列。
  • 易于实现和理解。
  • 对异常值具有鲁棒性。

优化克罗斯顿模型已成功用于预测广泛的间歇性时间序列，包括产品需求、服务需求和资源需求。

以下是优化克罗斯顿模型的一些特性：

  • 精确性： 研究表明，对于不规则需求的时间序列，优化克罗斯顿模型比经典克罗斯顿模型更为准确。这是因为，优化克罗斯顿模型使用概率函数来预测订单之间的平均间隔，其准确性优于历史订单的加权平均。
  • 适应性： 优化克罗斯顿模型对不同类型的间歇性时间序列更具适应性。这是因为，优化克罗斯顿模型使用概率函数来预测订单之间的平均间隔，能够适应不同的需求模式。
  • 易于实现和理解： 优化克罗斯顿模型易于实现和理解。这是因为，优化克罗斯顿模型是经典克罗斯顿模型的改进，经典模型广为人知且易于理解。
  • 鲁棒性： 优化克罗斯顿模型对异常值也表现出良好的鲁棒性。这是由于优化克罗斯顿模型使用概率函数来预测订单之间的平均间隔，从而可以忽略异常值。

加载库和数据

Tip

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

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

import matplotlib.pyplot as plt
import seaborn as sns
from statsmodels.graphics.tsaplots import plot_acf, 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/intermittend_demand2")

df.head()
date sales
0 2022-01-01 00:00:00 0
1 2022-01-01 01:00:00 10
2 2022-01-01 02:00:00 0
3 2022-01-01 03:00:00 0
4 2022-01-01 04:00:00 100

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 2022-01-01 00:00:00 0 1
1 2022-01-01 01:00:00 10 1
2 2022-01-01 02:00:00 0 1
3 2022-01-01 03:00:00 0 1
4 2022-01-01 04:00:00 100 1 
print(df.dtypes)
ds           object
y             int64
unique_id    object
dtype: object

我们可以看到我们的时间变量 (ds) 是以对象格式表示的，我们需要将其转换为日期格式。

df["ds"] = pd.to_datetime(df["ds"])

使用plot方法探索数据

使用StatsForecast类中的plot方法绘制一些系列。该方法从数据集中随机打印一个系列，适合进行基本的探索性数据分析(EDA)。

from statsforecast import StatsForecast

StatsForecast.plot(df)

自相关图

自相关（ACF）和偏自相关（PACF）图是用于分析时间序列的统计工具。ACF图显示时间序列值与其滞后值之间的相关性，而PACF图则显示在去除了之前滞后值影响后的时间序列值与其滞后值之间的相关性。

ACF和PACF图可以用来识别时间序列的结构，这在选择合适的时间序列模型时非常有帮助。例如，如果ACF图显示出重复的峰谷模式，这表明时间序列是平稳的，这意味着它在时间上具有相同的统计特性。如果PACF图显示出快速递减的波动模式，这表明时间序列是可逆的，这意味着它可以被逆转以获得平稳的时间序列。

ACF和PACF图的重要性在于它们可以帮助分析师更好地理解时间序列的结构。这种理解可以帮助选择合适的时间序列模型，从而提高预测未来时间序列值的能力。

分析ACF和PACF图的方法：

  • 查找图表中的模式。常见的模式包括重复的峰谷、锯齿形模式和平台模式。
  • 比较ACF和PACF图。PACF图通常比ACF图有更少的波动。
  • 考虑时间序列的长度。较长时间序列的ACF和PACF图将有更多的波动。
  • 使用置信区间。ACF和PACF图还显示自相关值的置信区间。如果某个自相关值在置信区间之外，则可能是显著的。
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 plotly.subplots import make_subplots
from statsmodels.tsa.seasonal import seasonal_decompose

def plotSeasonalDecompose(
    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
plotSeasonalDecompose(
    df["y"],
    model="additive",
    period=24,
    title="Seasonal Decomposition")

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

让我们将数据划分为几个集合

让我们将数据划分为几个集合 1. 用于训练我们的 Croston 优化模型 的数据。 2. 用于测试我们模型的数据

对于测试数据，我们将使用最近的 500 个小时来测试和评估我们模型的性能。

train = df[df.ds<='2023-01-31 19:00:00'] 
test = df[df.ds>'2023-01-31 19:00:00'] 
train.shape, test.shape
((9500, 3), (500, 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("Store visit");
plt.xlabel("Hours")
plt.show()

CrostonOptimized 的实现与 StatsForecast

要了解 CrostonOptimized Model 函数的参数，以下是它们的列表。有关更多信息，请访问 文档。

alias : str
    模型的自定义名称。

加载库

from statsforecast import StatsForecast
from statsforecast.models import CrostonOptimized

实例化模型

导入并实例化模型。设置参数有时会很棘手。罗伯·海德曼（Rob Hyndmann）撰写的关于季节周期的文章可能对season_length有所帮助。

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

# 我们称之为即将使用的模型
models = [CrostonOptimized()]

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

模型：一个模型列表。从模型中选择所需的模型并导入它们。

  • freq: 一个字符串，指示数据的频率。（请参见 pandas 可用的频率.）

  • n_jobs: n_jobs: int，表示并行处理使用的工作数量，使用 -1 表示所有核心。

  • fallback_model: 如果某个模型失败，将使用的模型。

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

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

拟合模型

# 拟合模型
sf.fit()
StatsForecast(models=[CrostonOptimized])

让我们看看我们 Croston 优化模型 的结果。我们可以通过以下指令进行观察：

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

预测方法

如果您希望在产品环境中提高速度，处理多个序列或模型，我们建议使用 StatsForecast.forecast 方法，而不是 .fit 和 .predict。

主要区别是 .forecast 不会存储拟合值，并且在分布式环境中具有很高的可扩展性。

预测方法接受两个参数：预测未来的 h（时间范围）和 level。

  • h (int): 表示预测未来 h 步。在这种情况下，预测 500 小时后。

这里的预测对象是一个新的数据框，包含一个模型名称的列和 y hat 值，以及不确定性区间的列。根据您的计算机，这一步预计需要大约 1 分钟。

Y_hat = sf.forecast(horizon)
Y_hat
ds CrostonOptimized
unique_id
1 2023-02-21 16:00:00 23.606695
1 2023-02-21 17:00:00 23.606695
1 2023-02-21 18:00:00 23.606695
... ... ...
1 2023-03-14 09:00:00 23.606695
1 2023-03-14 10:00:00 23.606695
1 2023-03-14 11:00:00 23.606695

500 rows × 2 columns

Y_hat=Y_hat.reset_index()
Y_hat
unique_id ds CrostonOptimized
0 1 2023-02-21 16:00:00 23.606695
1 1 2023-02-21 17:00:00 23.606695
2 1 2023-02-21 18:00:00 23.606695
... ... ... ...
497 1 2023-03-14 09:00:00 23.606695
498 1 2023-03-14 10:00:00 23.606695
499 1 2023-03-14 11:00:00 23.606695

500 rows × 3 columns

Y_hat1 = pd.concat([df,Y_hat])
Y_hat1
ds y unique_id CrostonOptimized
0 2022-01-01 00:00:00 0.0 1 NaN
1 2022-01-01 01:00:00 10.0 1 NaN
2 2022-01-01 02:00:00 0.0 1 NaN
... ... ... ... ...
497 2023-03-14 09:00:00 NaN 1 23.606695
498 2023-03-14 10:00:00 NaN 1 23.606695
499 2023-03-14 11:00:00 NaN 1 23.606695

10500 rows × 4 columns

fig, ax = plt.subplots(1, 1)
plot_df = pd.concat([df, Y_hat1]).set_index('ds')
plot_df['y'].plot(ax=ax, linewidth=2)
plot_df[ "CrostonOptimized"].plot(ax=ax, linewidth=2, color="yellow")
ax.set_title(' Forecast', fontsize=22)
ax.set_ylabel("Store visit (Hourly data)", fontsize=20)
ax.set_xlabel('Hours', fontsize=20)
ax.legend(prop={'size': 15})
ax.grid(True)

带置信区间的预测方法

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

预测方法接收两个参数：预测下一个 h（表示时间跨度）和 level。

  • h (int): 表示预测未来 h 步。在这种情况下，预测500小时后。

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

此步骤应少于1秒。

forecast_df = sf.predict(h=horizon) 
forecast_df
ds CrostonOptimized
unique_id
1 2023-02-21 16:00:00 23.606695
1 2023-02-21 17:00:00 23.606695
1 2023-02-21 18:00:00 23.606695
... ... ...
1 2023-03-14 09:00:00 23.606695
1 2023-03-14 10:00:00 23.606695
1 2023-03-14 11:00:00 23.606695

500 rows × 2 columns

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

pd.concat([df, forecast_df]).set_index('ds')
y unique_id CrostonOptimized
ds
2022-01-01 00:00:00 0.0 1 NaN
2022-01-01 01:00:00 10.0 1 NaN
2022-01-01 02:00:00 0.0 1 NaN
... ... ... ...
2023-03-14 09:00:00 NaN NaN 23.606695
2023-03-14 10:00:00 NaN NaN 23.606695
2023-03-14 11:00:00 NaN NaN 23.606695

10500 rows × 3 columns

df_plot= pd.concat([df, forecast_df]).set_index('ds').tail(5000)
df_plot
y unique_id CrostonOptimized
ds
2022-08-18 04:00:00 0.0 1 NaN
2022-08-18 05:00:00 80.0 1 NaN
2022-08-18 06:00:00 0.0 1 NaN
... ... ... ...
2023-03-14 09:00:00 NaN NaN 23.606695
2023-03-14 10:00:00 NaN NaN 23.606695
2023-03-14 11:00:00 NaN NaN 23.606695

5000 rows × 3 columns

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

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

plt.title("Store visit (Hourly data)");
plt.xlabel("Hourly")
plt.ylabel("Store visit")
plt.legend()
plt.show();

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

sf.plot(df, forecast_df)

交叉验证

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

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

下图展示了这种交叉验证策略：

执行时间序列交叉验证

时间序列模型的交叉验证被认为是最佳实践，但大多数实现非常缓慢。statsforecast库将交叉验证实现为分布式操作，从而使该过程的执行不那么耗时。如果您拥有大型数据集，您也可以使用Ray、Dask或Spark在分布式集群中执行交叉验证。

在这种情况下，我们希望评估每个模型在过去5个月的表现 (n_windows=)，每隔两个月进行一次预测 (step_size=50)。根据您的计算机，这一步应大约需要1分钟。

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

  • df: 训练数据框

  • h (int): 表示未来h步的预测。在这种情况下，是500小时。

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

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

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

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

  • unique_id: 索引。如果您不喜欢使用索引，只需运行 crossvalidation_df.resetindex()。
  • ds: 日期戳或时间索引
  • cutoff: n_windows 的最后日期戳或时间索引。
  • y: 真实值
  • model: 包含模型名称和拟合值的列。
crossvalidation_df
ds cutoff y CrostonOptimized
unique_id
1 2023-01-23 12:00:00 2023-01-23 11:00:00 0.0 23.655830
1 2023-01-23 13:00:00 2023-01-23 11:00:00 0.0 23.655830
1 2023-01-23 14:00:00 2023-01-23 11:00:00 0.0 23.655830
... ... ... ... ...
1 2023-02-21 13:00:00 2023-01-31 19:00:00 60.0 27.418417
1 2023-02-21 14:00:00 2023-01-31 19:00:00 20.0 27.418417
1 2023-02-21 15:00:00 2023-01-31 19:00:00 20.0 27.418417

2500 rows × 4 columns

模型评估

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

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

计算RMSE的函数接受两个参数：

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

致谢

我们要感谢 Naren Castellon 编写本教程。

参考文献

  1. Changquan Huang • Alla Petukhina. 施普林格系列（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