时间序列预测的不确定性区间估计：基于EnbPI的方法与应用研究|估计器|子集|应用研究|时间序列|算法|观测值

在现代预测分析领域，准确评估预测结果的不确定性已成为一个关键挑战。预测的不确定性量化不仅能够提供更可靠的决策支持，还能深入揭示模型的预测能力边界。本文聚焦于时间序列预测中的不确定性量化问题，重点探讨基于一致性预测理论的集成批量预测区间（Ensemble Batch Prediction Interval, EnbPI）方法。

传统一致性预测方法的核心假设——数据可交换性（exchangeability）——在时间序列分析中面临重大挑战。时间序列数据本质上包含了重要的时序依赖特征，如趋势、周期性和季节性模式，这使得观测值的顺序信息对预测至关重要。因此，在构建时间序列预测的不确定性量化框架时，必须保持数据的时序完整性。

EnbPI的工作原理

2021年提出的EnbPI方法为解决这一技术难题提供了创新解决方案。该方法突破了传统一致性预测对数据可交换性的依赖，能够有效处理非平稳时间序列数据和复杂的时空依赖关系。

EnbPI的核心思想是在数据的不同子集上训练基础模型，并通过集成这些模型来构建预测区间。EnbPI通过自举采样（Bootstrap Sampling）方法创建数据子集，即对时间序列进行有放回随机采样。

作为时间序列领域的一致性预测扩展，EnbPI具有以下技术优势：

构建分布无关的预测区间，实现近似边际覆盖率
支持任意估计器的集成
推理阶段计算效率高，无需模型重训练
适用于小规模数据集，不需要额外的校准集

需要注意的是，由于需要训练多个模型，EnbPI的训练阶段会消耗较多计算资源。

EnbPI的核心技术框架

EnbPI方法的核心思想是通过集成学习范式构建预测区间。该方法通过自举采样技术在时间序列数据上构建多个子集，并在这些子集上训练独立的预测模型。这种方法具有以下技术优势：

方法论优势

构建分布无关的预测区间
实现近似边际覆盖率
保持时序数据的结构特征

实践优势

支持任意基础预测器的集成
推理阶段计算效率高
无需额外的模型重训练
适用于小规模数据集分析

技术局限

训练阶段计算开销较大
需要针对具体应用场景优化集成规模

这种系统性的方法论框架为时间序列预测的不确定性量化提供了坚实的理论基础，同时也为实践应用提供了可操作的技术路径。

EnbPI算法实现流程

EnbPI算法可分为训练和预测两个主要阶段，每个阶段包含两个关键步骤。在开始训练EnbPI集成之前，需要首先选择合适的基础估计器。该估计器可以是任何机器学习模型，包括梯度提升树、神经网络、线性回归或传统统计模型。

EnbPI算法流程图。黄色框表示训练阶段的步骤，蓝色框表示预测阶段的步骤，绿色框表示可选的在线更新步骤。

1、训练阶段详解

训练阶段主要包含两个核心步骤：从训练数据中生成非重叠子集，并在每个子集上训练模型。

第1步：自举子集采样

在生成子集之前，首先需要确定集成模型的规模。每个集成成员都需要一个独立的自举样本。

子集采样的关键要求是保持非重叠性。这意味着每个子集都是独特的，且与其他子集相互独立。每个子集包含原始数据集中的不同数据片段。通过在不同数据子集上训练模型，能够引入模型间的变异性，从而增强集成的多样性并降低过拟合风险。

集成规模的选择需要根据数据集大小进行权衡。数据集较小时，可生成的非重叠子集数量有限，因为每个子集都需要足够的数据量来保证模型训练的有效性。较小的集成规模会导致多样性不足，从而影响EnbPI方法的性能，具体表现为预测区间偏宽。

增加模型数量通常能提升性能表现，产生更窄的预测区间。这是因为更大规模的集成能够捕获数据中更丰富的变异性，同时通过聚合更多预测来提高模型稳健性。这也意味着更高的计算成本。

实践表明，20到50个模型的集成规模通常能够在计算效率和预测准确性之间取得良好的平衡。

确定集成规模后，需要为每个集成成员创建相应的数据子集。子集通过从原始数据集中有放回抽样获得。

值得注意的是，为了保持时间序列数据的时序依赖性，采样单位为数据块而非单个观测值。由于采用有放回采样，某些数据块可能在子集中重复出现，而其他块可能未被选中。

时间序列自举采样流程示意图。原始时间序列被分割为多个数据块，每个子集通过随机有放回采样这些数据块构建而成，数据块可能在子集内重复出现。

完成自举采样后，即可进入模型训练阶段。

第2步：训练自举集成模型

对每个自举子集训练一个独立的模型，从而构建一个具有多样性的模型集成。由于各个模型在不同数据子集上训练，即使面对相同的输入数据也会产生不同的预测结果。这种预测的多样性是获得稳健预测区间估计的关键因素。

2、预测阶段实现机制

第3步：留一法（LOO）估计实现

在完成集成模型训练后，需要评估集成模型在未知数据上的预测方差，用于校准集成并确定预测区间的宽度。在EnbPI框架中，非一致性分数定义为真实值与预测值之间的偏差。这里需要解决的关键问题是校准数据集的选择。

由于集成中的每个估计器都在训练集的不同子集上进行训练，没有任何单个估计器接触过完整的训练数据。这一特性使得我们可以直接利用训练集进行校准，无需额外的校准数据集。

对训练集中的每个观测值，采用集成进行预测，但仅使用训练过程中未接触过该观测值的估计器。随后通过聚合函数（如均值或中位数）对这些模型的预测结果进行汇总。

聚合函数的选择会影响预测的稳健性。均值聚合虽然对异常值较为敏感，但能够有效降低整体误差；中位数聚合则具有较强的异常值抵抗能力，适用于噪声较大的数据场景。因此聚合函数的选择需要根据具体的数据特征和应用需求进行确定。

通过聚合预测可以计算每个观测值的非一致性分数。这种方法实质上利用了样本外误差作为非一致性度量，为后续未知数据的预测校准提供基础。

EnbPI集成示例（包含4个基础模型）的校准过程图解。对于训练集中的每个观测值，仅使用训练过程中未见过该观测值的模型进行预测。例如，对第一个观测值仅使用模型3和4进行预测（因为该观测值出现在模型1和2的训练子集中）。预测结果通过均值聚合（紫色点）。对第二个观测值可使用模型1、3和4，因为仅模型2的训练集包含该观测值。依此类推完成所有训练集观测值的预测。通过计算集成均值与真实值之间的误差得到非一致性分数。最后根据分数分布和预设的显著性水平α确定截断值q。

第4步：预测区间构建方法

对于未知数据的预测，使用所有训练完成的集成估计器进行预测，并采用与第3步相同的聚合函数对单个预测结果进行汇总，得到的值作为预测区间的中心点。

预测区间的构建基于第3步获得的残差分布。通过预先设定的显著性水平确定截断值，将该值在预测中心点的基础上进行加减运算，从而得到预测区间的上下界。

第5步：非一致性分数的动态更新机制（可选）

上述方法使用的非一致性分数仅基于训练数据计算，在接收到新数据后并不进行更新。这导致预测区间的宽度保持不变，无法反映数据分布或模型性能的动态变化。此限制在基础数据特征或模型表现发生显著变化时可能带来问题。

为解决这一问题，可以在获取新观测值时动态更新非一致性分数。这一过程无需重新训练集成估计器，仅需重新计算非一致性分数。通过这种机制，预测区间能够及时反映最新的数据特征和模型动态，实现区间宽度的自适应调整。

3、EnbPI方法的实验验证

为了验证EnbPI方法的有效性，我们这里选择了德国电力价格预测这一实际应用场景。在实验数据中，我们特意在最后两周的数据中引入了100欧元/兆瓦时的人工偏移，用于评估预测区间对突发性变化的适应能力。

通过初步数据分析，我们观察到电力价格数据具有明显的双重周期性特征：

日内周期：由于用电负载的日内分布规律，价格在早晚高峰时段显著上升
周度周期：工作日的电力价格普遍高于周末，这反映了工商业用电需求的周期性变化

基于对数据特征的理解，我们构建了以下特征集：

时间特征：包括日、月、年等基础时间维度，以及工作日/周末标识
历史数据特征：包含过去7天的同时段历史价格，以及基于过去24小时的滑动平均价格（滞后一天）

以下是特征工程的具体实现：

# %% 特征工程实现
forecast_horizon = 24 # 预测时间窗口设置为24小时
# 时间特征提取
price_data["year"] = price_data["dt"].dt.year
price_data["month"] = price_data["dt"].dt.month
price_data["day"] = price_data["dt"].dt.day
# 周末标识特征构建
price_data["is_weekend"] = (price_data["dt"].dt.dayofweek >= 5).astype(int)
# 历史同期价格特征构建
for day in range(forecast_horizon, 7 * forecast_horizon, forecast_horizon):
price_data[f"hour_lag_{day}"] = price_data["y"].shift(day)
# 24小时滑动平均价格特征
price_data["ma_24"] = (
price_data["y"].rolling(window=24, center=False, closed="left").mean().shift(24)
)
price_data.set_index("dt", inplace=True)

为了全面评估模型性能，我们采用了较长的测试周期。虽然实际应用中通常只需要预测未来24小时的价格，但我们选择使用最后30天的数据作为测试集，这使我们能够更好地观察和评估预测区间的动态变化特性。

# %% 数据集划分
test_horizon = forecast_horizon * 30 # 测试集长度设为30天
X_train = price_data.iloc[:-test_horizon].drop(columns=["y"])
y_train = price_data["y"].iloc[:-test_horizon]
X_test = price_data.iloc[-test_horizon:].drop(columns=["y"])
y_test = price_data["y"].iloc[-test_horizon:]
_, ax = plt.subplots(1, figsize=(10, 3))
ax.plot(y_train, label="train")
ax.plot(y_test, label="test")
ax.set_ylabel("price in EUR/MWh")
ax.legend()
ax.grid()

为了便于结果可视化和性能评估，我们实现了一个专门的结果分析函数：

# %% 结果分析与可视化函数
def plot_results(y_test, y_pred, y_pred_int):
fig, ax = plt.subplots(figsize=(10, 3))
ax.plot(y_test, label="Actual")
ax.plot(y_test.index, y_pred, label="Predicted", ls="--")
ax.fill_between(
y_test.index, y_pred_int[:, 0, 0], y_pred_int[:, 1, 0], color="green", alpha=0.2
)
ax.set_ylabel("price in EUR/MWh")
ax.legend(loc="best")
ax.grid()
fig.tight_layout()
# 计算覆盖率和区间宽度指标
coverage = regression_coverage_score(y_test, y_pred_int[:, 0, 0], y_pred_int[:, 1, 0])
width_interval = regression_mean_width_score(y_pred_int[:, 0, 0], y_pred_int[:, 1, 0])
print(f"Coverage: {coverage:}")
print(f"Width of the interval: {width_interval}")

基础模型选择与集成实现

为了简单起见，我们选择LightGBM作为基础预测模型。需要说明的是，虽然本文重点不在于优化基础模型性能，但在实际应用中，基础模型的预测准确性直接影响预测区间的质量。

EnbPI集成模型的构建与训练

EnbPI的实现采用了模块化的方式。首先，使用Mapie的BlockBootsrap类生成自举样本。我们设置了20个不重叠的数据块，每个块的长度与预测时域（24小时）相匹配。

EnbPI模型的具体实现

以下代码展示了EnbPI模型的完整实现过程。使用Mapie库的MapieTimeSeriesRegressor类来封装模型，并通过参数配置实现预测功能：

# %% EnbPI模型实现
cv_mapie_ts = BlockBootstrap(n_resamplings=20, length=24, overlapping=False, random_state=42)
params = {"objective": "rmse"}
model = lgb.LGBMRegressor(**params)
mapie_enbpi = MapieTimeSeriesRegressor(
model,
method="enbpi",
cv=cv_mapie_ts,
agg_function="mean",
n_jobs=-1,
)
mapie_enbpi = mapie_enbpi.fit(X_train, y_train)
# 设置显著性水平为0.05，对应95%的置信水平
alpha = 0.05
y_pred, y_pred_int = mapie_enbpi.predict(X_test, alpha=alpha, ensemble=True, optimize_beta=True)
plot_results(y_test, y_pred, y_pred_int)

4、实验结果分析

实验结果显示，模型达到了91%的覆盖率，略低于预期的95%目标值，预测区间的平均宽度为142.62欧元/兆瓦时。覆盖率未达预期的主要原因可能是测试期间出现的价格突变。预测区间宽度在整个测试期间保持相对稳定。

5、非一致性分数的动态更新实现

为了提高模型对数据变化的适应能力，还可以实现了基于partial_fit()方法的动态更新机制。这种实现模拟了实际场景中逐步获取新数据的情况：

# %% 实现带动态更新的EnbPI模型
cv_mapie_ts = BlockBootstrap(n_resamplings=20, length=24, overlapping=False, random_state=42)
model = lgb.LGBMRegressor(**params)
mapie_enbpi = MapieTimeSeriesRegressor(
model,
method="enbpi",
cv=cv_mapie_ts,
agg_function="mean",
n_jobs=-1,
)
mapie_enpbi = mapie_enbpi.fit(X_train, y_train)
# 初始化预测结果存储数组
y_pred_pfit = np.zeros(y_pred.shape)
y_pred_int_pfit = np.zeros(y_pred_int.shape)
# 首次预测
y_pred_pfit[:forecast_horizon], y_pred_int_pfit[:forecast_horizon, :, :] = mapie_enbpi.predict(
X_test.iloc[:forecast_horizon, :], alpha=alpha, ensemble=True, optimize_beta=True
)
# 逐步更新预测
for step in range(forecast_horizon, len(X_test), forecast_horizon):
# 使用新数据更新模型
mapie_enbpi.partial_fit(
X_test.iloc[(step - forecast_horizon) : step, :],
y_test.iloc[(step - forecast_horizon) : step],
)
# 生成下一时段预测
(
y_pred_pfit[step : step + forecast_horizon],
y_pred_int_pfit[step : step + forecast_horizon, :, :],
) = mapie_enbpi.predict(
X_test.iloc[step : (step + forecast_horizon), :],
alpha=alpha,
ensemble=True,
optimize_beta=True,
)
plot_results(y_test, y_pred, y_pred_int)