企业官网建设流程全解析

简介：本资源是一套面向计算机、电子信息工程及数学专业本科生的时序预测实战方案，聚焦SARIMA与LSTM混合建模这一前沿方法，解决具有强季节性与非线性特征的时间序列高精度预测问题，适用于课程设计、期末大作业及毕业设计等实践场景。压缩包共3个文件（2个CSV数据集用于训练与验证，1个Python主程序含完整SARIMA-LSTM联合建模流程），总大小仅51KB，轻量易部署，适配Anaconda+PyCharm+TensorFlow环境。已有742人学习下载，体现其在教学与入门科研中的广泛认可。代码采用全参数化设计，支持模型结构、超参、滑动窗口长度等灵活调整；每行均附保姆级中文注释，清晰呈现数据预处理、SARIMA残差提取、LSTM特征融合、多步预测及结果可视化全流程；配套焦作市历史气象/经济类时序数据，开箱即用，显著降低算法复现门槛。

1. 为什么单用 SARIMA 或 LSTM 都会翻车？——当电力负荷突增、电商销量暴增、气象数据跳变时，季节性+非线性+长记忆三重压力下，必须把统计模型的可解释骨架和神经网络的拟合弹性焊死在一起

你手头有一份连续 3 年每小时的光伏电站发电功率数据：夏季午后有强季节性峰值，但某天突遇沙尘暴，功率断崖式下跌；你正在做双十二前 7 天的订单量预测，历史规律清晰，可今年平台突然上线新补贴策略，第 3 天起流量结构彻底改变；你维护城市地铁客流系统，工作日早高峰稳定，但某次大型展会期间，周末客流反超工作日，且持续 5 天。这些场景里，SARIMA 能抓住“每周一早8点必涨”的统计规律，却对沙尘暴、补贴策略、展会这种一次性结构突变束手无策；LSTM 能学出突变后的模式，但若训练数据只含 1 次展会，它大概率把“展会效应”当成噪声滤掉，更无法告诉你“这个突变是因展会导致，下周恢复常态”。本方案不是简单拼接两个模型，而是让 SARIMA 先做可审计的残差分解：把原始序列拆成「确定性季节趋势 + 统计平稳残差」，再把残差喂给 LSTM 去捕捉 SARIMA 漏掉的非线性动态依赖。最终预测 = SARIMA 趋势预测 + LSTM 对残差的精细化修正。这不是玄学调参，而是工程上可复现、可回溯、可解释的混合建模路径——适合已掌握 Python 基础、跑过单个 SARIMA 或 LSTM 但预测误差始终卡在 8%~12% 的一线数据工程师与算法工程师。如果你正被业务方追问“为什么昨天预测偏高？是模型问题还是数据问题？”，这套方法能让你指着残差图说清楚。

2. 从原始数据到可建模序列：清洗、检验、分解三步不可跳过的硬门槛

时间序列预测不是把数据丢进模型就完事。我见过太多人直接拿 raw.csv 里的“日期,销量”两列跑 SARIMA，结果 AIC 炸到 1e6，残差自相关图（ACF）拖着长尾巴像没剪过的指甲。根源全在预处理没做透。下面这三步，少一步，后面所有代码都是在给错误打补丁。

2.1 数据清洗：缺失值、异常值、时间戳对齐的实操铁律

真实工业数据里，传感器离线、数据库写入失败、人工补录错误是常态。不能简单用df.fillna(method='ffill')一填了事——那会把 2 小时的断连伪装成平滑过渡，污染季节性特征。正确做法是分层处理：

import pandas as pd import numpy as np from scipy import stats def robust_clean_ts(df, time_col='timestamp', value_col='power_kw'): """ df: 原始DataFrame，必须含时间列和数值列 time_col: 时间列名（支持 '2023-01-01 08:00' 或 1672531200 时间戳） value_col: 待预测数值列名 """ # 步骤1：强制转为datetime并设索引（关键！SARIMA要求等间隔索引） df = df.copy() if not pd.api.types.is_datetime64_any_dtype(df[time_col]): df[time_col] = pd.to_datetime(df[time_col]) df = df.set_index(time_col).sort_index() # 步骤2：检测并标记缺失时段（非简单NaN，而是整段数据丢失） freq_inferred = pd.infer_freq(df.index) # 自动推断频率，如 'H'（小时）、'D'（天） if freq_inferred is None: raise ValueError(f"无法推断时间频率，请检查数据是否等间隔。当前索引前5条：{df.index[:5]}") # 生成完整时间索引（覆盖实际数据范围） full_index = pd.date_range(start=df.index.min(), end=df.index.max(), freq=freq_inferred) df_full = df.reindex(full_index) # 缺失位置补NaN # 步骤3：异常值检测（用滚动窗口Z-score，比全局标准差更鲁棒） window_size = 24 * 7 # 一周的小时数，适应季节性 rolling_mean = df_full[value_col].rolling(window=window_size, min_periods=window_size//2).mean() rolling_std = df_full[value_col].rolling(window=window_size, min_periods=window_size//2).std() z_scores = np.abs((df_full[value_col] - rolling_mean) / (rolling_std + 1e-8)) # 防除零 # 标记异常值（Z-score > 3 且非NaN处） outlier_mask = (z_scores > 3) & df_full[value_col].notna() print(f"检测到 {outlier_mask.sum()} 个异常值点（基于滚动Z-score）") # 步骤4：用临近均值插补（非线性插值易扭曲季节峰谷） df_clean = df_full.copy() df_clean.loc[outlier_mask, value_col] = np.nan # 用前后各3个有效点的均值填充（避免用未来信息泄露） df_clean[value_col] = df_clean[value_col].interpolate( method='linear', limit_direction='both', limit=3 ) return df_clean # 使用示例（假设你有 data_raw.csv） # df_raw = pd.read_csv('data_raw.csv') # df_clean = robust_clean_ts(df_raw, time_col='dt', value_col='load_mw')

提示：pd.infer_freq()是关键。如果返回None，说明你的数据存在不规则间隔（如某天缺了 3 个点），必须先用resample('H').mean()下采样或asfreq('H')强制对齐，否则 SARIMA 的seasonal_order参数会完全失效。这是新手最常踩的第一个坑。

2.2 平稳性检验：ADF 和 KPSS 不是二选一，而是互为验证的双保险

SARIMA 要求残差序列平稳，但很多人只跑一次 ADF 就收工。问题在于：ADF 原假设是“非平稳”，拒绝原假设才说明平稳；KPSS 原假设是“平稳”，拒绝原假设才说明非平稳。二者结论冲突时，必须深挖原因。常见情况是：ADF 显著（p<0.05）但 KPSS 也显著（p<0.05）——这往往意味着序列存在确定性趋势（deterministic trend）而非随机游走（stochastic trend），此时应做去趋势（detrend）而非差分（differencing）。

from statsmodels.tsa.stattools import adfuller, kpss import matplotlib.pyplot as plt def check_stationarity(series, max_diff=2): """ 同时运行ADF和KPSS检验，返回建议的差分阶数d series: pd.Series，已清洗的单变量时间序列 """ print("=== 初始序列 ADF 检验 ===") adf_result = adfuller(series.dropna()) print(f"ADF Statistic: {adf_result[0]:.4f}") print(f"p-value: {adf_result[1]:.4f}") print(f"Critical Values: {adf_result[4]}") print("\n=== 初始序列 KPSS 检验 ===") kpss_result = kpss(series.dropna(), regression='c') # 'c' 表示含常数项 print(f"KPSS Statistic: {kpss_result[0]:.4f}") print(f"p-value: {kpss_result[1]:.4f}") print(f"Critical Values: {kpss_result[3]}") # 决策逻辑：只有当 ADF p<0.05 AND KPSS p>0.05 时，才认为平稳 d_suggested = 0 for d in range(1, max_diff+1): diff_series = series.diff(d).dropna() adf_p = adfuller(diff_series)[1] kpss_p = kpss(diff_series, regression='c')[1] print(f"\n--- 差分 {d} 次后 ---") print(f"ADF p-value: {adf_p:.4f}, KPSS p-value: {kpss_p:.4f}") if adf_p < 0.05 and kpss_p > 0.05: d_suggested = d print(f"✅ 建议差分阶数 d = {d}") break if d_suggested == 0: print("⚠️ 注意：未找到同时满足 ADF 和 KPSS 的平稳阶数，建议检查是否存在确定性趋势") # 绘制趋势图辅助判断 plt.figure(figsize=(12, 4)) plt.plot(series.index, series.values, label='Original') plt.title('原始序列趋势图（检查确定性趋势）') plt.legend() plt.show() return d_suggested # 使用示例 # d = check_stationarity(df_clean['load_mw'])

注意：kpss(..., regression='c')中的'c'表示检验时假设序列含常数项（即均值非零），这是绝大多数业务序列的实际情况。若误用'nc'（无常数项），KPSS 会过度拒绝平稳性，导致不必要的差分，损失信息。

2.3 季节性分解：用 STL 而非经典分解，避开周期长度硬编码陷阱

SARIMA 的seasonal_order要求你提前指定季节周期s（如小时数据 s=24，日数据 s=7）。但真实数据中，季节性可能漂移：夏季空调负荷周期是 24 小时，冬季取暖负荷可能变成 12 小时峰谷。经典分解（seasonal_decompose）强制使用固定s，会把漂移当成噪声。STL（Seasonal-Trend decomposition using Loess）用局部加权回归自动适应周期变化，是工业级首选。

from statsmodels.tsa.seasonal import STL def stl_decompose(series, period=None, robust=True): """ 使用STL进行稳健季节性分解 period: 若为None，则自动推断（推荐） robust: 是否启用鲁棒拟合（对异常值不敏感） """ if period is None: # 自动推断周期：取 ACF 第一个显著峰对应的滞后阶数 from statsmodels.tsa.stattools import acf acf_vals = acf(series.dropna(), nlags=min(100, len(series)//2)) # 找第一个超过 0.3 的滞后（经验阈值） peaks = np.where(acf_vals > 0.3)[0] if len(peaks) > 1: period = peaks[1] # 跳过 lag=0 else: period = 24 # 默认小时数据用24 print(f"STL 分解采用周期 s = {period}") stl = STL(series, period=period, robust=robust, seasonal_deg=1) result = stl.fit() # 可视化分解结果（关键！必须看图确认分解质量） fig, axes = plt.subplots(4, 1, figsize=(12, 10)) axes[0].plot(result.observed) axes[0].set_ylabel('Observed') axes[1].plot(result.trend) axes[1].set_ylabel('Trend') axes[2].plot(result.seasonal) axes[2].set_ylabel('Seasonal') axes[3].plot(result.resid) axes[3].set_ylabel('Residual') plt.tight_layout() plt.show() return result # 使用示例 # stl_result = stl_decompose(df_clean['load_mw']) # trend = stl_result.trend # seasonal = stl_result.seasonal # resid = stl_result.resid # 这就是喂给LSTM的残差序列！

血泪经验：STL 的seasonal_deg=1（季节性部分用线性拟合）比默认的0（常数拟合）更能捕捉季节峰谷的缓慢漂移。若你的数据季节峰宽在变（如早高峰从 7:30-9:00 慢慢移到 8:00-9:30），必须设为 1，否则残差里会残留强季节性，LSTM 学不会。

3. SARIMA 模型构建：AIC/BIC 不是终点，残差诊断图才是生死线

很多人以为model.fit()成功就万事大吉，导出forecast()结果交差。但 SARIMA 的核心价值不在预测值本身，而在它能否把原始序列中可解释的统计结构剥离干净，留下真正需要 LSTM 攻克的“硬骨头”——即残差。如果 SARIMA 拟合得不好，LSTM 就是在修一座地基歪斜的大楼。

3.1 自动参数搜索：用pmdarima.auto_arima但必须关掉它的“黑匣子”模式

pmdarima.auto_arima是神器，但它默认开启stepwise=True和maxiter=50，会跳过大量潜在好参数组合。更致命的是，它默认用information_criterion='aic'，但在小样本（<2000 点）时，AIC 容易过拟合，BIC 更稳健。我们必须手动控制搜索空间，并强制输出所有候选模型。

import pmdarima as pm def sarima_search_and_fit(series, d=None, D=None, m=None, max_p=3, max_q=3, max_P=2, max_Q=2): """ 安全的 SARIMA 参数搜索：显式控制搜索范围，禁用黑匣子优化 series: 已平稳的序列（经2.2节处理） d, D, m: 若已知，直接传入；否则为None由auto_arima推断 """ # 步骤1：先做一次快速扫描，获取基础d, D, m if d is None or D is None or m is None: print("🔍 正在推断差分阶数d、季节差分阶数D和周期m...") # 用简单方法推断m（避免auto_arima内部复杂逻辑） from statsmodels.tsa.stattools import adfuller # 对原始序列做ADF，若不平稳，尝试d=1 if adfuller(series)[1] > 0.05: d = 1 else: d = 0 # 推断m：用FFT找主频（比ACF更抗噪） from scipy.fft import fft n = len(series) freqs = fft(series) magnitudes = np.abs(freqs[:n//2]) # 找最大幅值对应的频率（排除DC分量） peak_idx = np.argmax(magnitudes[1:]) + 1 m = int(round(n / peak_idx)) if peak_idx > 0 else 24 # D默认为0，除非m已知且序列明显有季节性非平稳 D = 0 print(f"设定参数：d={d}, D={D}, m={m}") # 步骤2：执行受限搜索（关闭stepwise，确保遍历） model = pm.auto_arima( series, start_p=0, start_q=0, start_P=0, start_Q=0, max_p=max_p, max_q=max_q, max_P=max_P, max_Q=max_Q, m=m, seasonal=True, d=d, D=D, trace=True, # 必须开！看它到底试了哪些组合 error_action='ignore', suppress_warnings=True, stepwise=False, # 关键！禁用启发式跳过 information_criterion='bic', # 小样本用BIC更稳 n_jobs=-1 ) print(f"\n🏆 最优模型: {model.order} {model.seasonal_order}") return model # 使用示例（假设resid是STL分解出的残差） # sarima_model = sarima_search_and_fit(resid, d=0, D=0, m=24)

逻辑说明：stepwise=False强制穷举所有(p,q,P,Q)组合，虽然慢一点，但能保证不漏掉 BIC 更优的模型。trace=True输出的每一行都包含p,d,q,P,D,Q,s,BIC，你可以复制粘贴到 Excel 里排序，手动挑出 BIC 最小的 3 个，再对比它们的残差图——这才是工程师该干的事，不是把决策权交给黑匣子。

3.2 残差诊断：三张图定生死，缺一不可

SARIMA 拟合后，必须立刻画三张图：残差时序图、残差 ACF 图、残差 Q-Q 图。任何一张不合格，LSTM 输入的就是一团乱麻。

def sarima_diagnostics(model, series): """ SARIMA 残差诊断三件套 model: fitted SARIMA model series: 原始输入序列（用于计算残差） """ # 获取残差 residuals = model.resid fig, axes = plt.subplots(2, 2, figsize=(14, 10)) # 1. 残差时序图：看是否有明显趋势或异方差 axes[0, 0].plot(residuals) axes[0, 0].set_title('Residuals over Time') axes[0, 0].axhline(y=0, color='r', linestyle='--') # 2. 残差直方图 + KDE：看是否近似正态 axes[0, 1].hist(residuals, bins=50, density=True, alpha=0.6, label='Histogram') from scipy.stats import norm mu, std = norm.fit(residuals) x = np.linspace(residuals.min(), residuals.max(), 100) axes[0, 1].plot(x, norm.pdf(x, mu, std), 'r-', lw=2, label=f'Normal fit ($\mu$={mu:.2f}, $\sigma$={std:.2f})') axes[0, 1].set_title('Residuals Histogram & Normal Fit') axes[0, 1].legend() # 3. 残差 ACF：看是否白噪声（所有滞后阶数的ACF应在±2/√n内） from statsmodels.graphics.tsaplots import plot_acf plot_acf(residuals, ax=axes[1, 0], lags=40) axes[1, 0].set_title('Residuals ACF') # 4. Q-Q 图：正态性终极检验 from statsmodels.graphics.gofplots import qqplot qqplot(residuals, line='s', ax=axes[1, 1]) axes[1, 1].set_title('Q-Q Plot of Residuals') plt.tight_layout() plt.show() # 统计检验 from statsmodels.stats.diagnostic import acorr_ljungbox lb_test = acorr_ljungbox(residuals, lags=[10, 20], return_df=True) print("\n=== Ljung-Box 检验（检验白噪声）===") print(lb_test) print("✅ 若所有p-value > 0.05，则接受白噪声假设") # 使用示例 # sarima_diagnostics(sarima_model, resid)

参数说明：acorr_ljungbox(..., lags=[10,20])检验前 10 阶和前 20 阶滞后是否整体无自相关。如果lb_test['lb_pvalue'].min() < 0.05，说明残差仍有可提取的模式，SARIMA 拟合不足，必须回到 3.1 节扩大max_p/max_q范围重搜。这是决定是否进入 LSTM 阶段的闸门。

4. LSTM 模型构建：不是堆层数，而是用残差的“记忆长度”倒推时间窗

把 SARIMA 残差喂给 LSTM，目标很明确：学出残差中 SARIMA 捕捉不到的长程非线性依赖。但很多人直接用timesteps=10硬编码，结果模型要么记不住上周同时间的影响，要么塞进太多无关信息引入噪声。正确做法是：用残差的自相关衰减长度来确定 LSTM 的输入窗口。

4.1 确定 LSTM 输入窗口：用残差 ACF 的 e-folding time

SARIMA 残差理论上应是白噪声（ACF 在 lag=1 后立即归零），但实际总有微弱拖尾。我们取 ACF 衰减到初始值1/e ≈ 0.368时的滞后阶数，作为 LSTM 应记住的最小时间跨度。这比凭经验设timesteps=24科学得多。

from statsmodels.tsa.stattools import acf import numpy as np def get_lstm_timesteps(residuals, threshold=0.368, max_lag=200): """ 根据残差ACF衰减到1/e的时间，确定LSTM输入步长 residuals: SARIMA拟合后的残差序列 threshold: e-folding threshold (1/e) """ # 计算ACF acf_vals = acf(residuals, nlags=max_lag) # 找第一个低于threshold的lag（从lag=1开始，跳过lag=0） for lag in range(1, len(acf_vals)): if abs(acf_vals[lag]) < threshold: timesteps = lag print(f"📊 残差ACF在 lag={lag} 处衰减至 {acf_vals[lag]:.3f} < {threshold}") break else: timesteps = max_lag print(f"⚠️ ACF未在 {max_lag} 阶内衰减至 {threshold}，采用最大窗口 {max_lag}") # 经验修正：向上取整到最近的周期倍数（如小时数据取24的倍数） if timesteps < 24: timesteps = 24 elif timesteps < 48: timesteps = 48 else: timesteps = ((timesteps // 24) + 1) * 24 print(f"✅ 最终LSTM输入窗口 timesteps = {timesteps}") return timesteps, acf_vals # 使用示例 # timesteps, acf_vals = get_lstm_timesteps(sarima_model.resid)

为什么取 24 的倍数？因为我们的残差来自小时级数据，而人类活动、设备启停都有天然的 24 小时节律。即使 ACF 在 lag=35 衰减，设timesteps=48也能让 LSTM 同时看到“今天凌晨”和“昨天凌晨”的对比，这对捕捉跨日影响至关重要。这是领域知识对纯统计的修正。

4.2 构建可复现的 LSTM 模型：Keras 函数式 API + EarlyStopping 硬约束

用 Sequential API 容易写出不可复现的模型。函数式 API 能清晰定义输入、状态传递和输出，且便于后续加入注意力机制。关键约束：必须用EarlyStopping防止过拟合，且patience=10——因为残差信噪比低，训练太早停会欠拟合，太晚停会过拟合。

import tensorflow as tf from tensorflow.keras.models import Model from tensorflow.keras.layers import Input, LSTM, Dense, Dropout, BatchNormalization from tensorflow.keras.callbacks import EarlyStopping, ReduceLROnPlateau def build_lstm_model(input_shape, units=50, dropout_rate=0.2, output_dim=1): """ 构建LSTM模型：输入shape=(timesteps, 1)，输出单步预测 input_shape: (timesteps, 1) 如 (48, 1) units: LSTM单元数（50是小时数据的经验起点） """ inputs = Input(shape=input_shape) # 第一层LSTM：return_sequences=True，为后续层提供时序输出 x = LSTM(units, return_sequences=True, dropout=dropout_rate, recurrent_dropout=dropout_rate)(inputs) x = BatchNormalization()(x) # 第二层LSTM：return_sequences=False，压缩为向量 x = LSTM(units//2, return_sequences=False, dropout=dropout_rate, recurrent_dropout=dropout_rate)(x) x = BatchNormalization()(x) # 全连接层 x = Dense(units//2, activation='relu')(x) x = Dropout(dropout_rate)(x) outputs = Dense(output_dim, activation='linear')(x) # 线性激活，适配任意尺度残差 model = Model(inputs=inputs, outputs=outputs) # 编译：用MAE而非MSE，因残差含异常值，MAE更鲁棒 model.compile( optimizer=tf.keras.optimizers.Adam(learning_rate=0.001), loss='mae', # Mean Absolute Error metrics=['mae'] ) return model # 创建模型 # lstm_model = build_lstm_model(input_shape=(timesteps, 1)) # 训练回调 callbacks = [ EarlyStopping( monitor='val_loss', patience=10, # 连续10轮不下降则停 restore_best_weights=True # 关键！必须恢复最优权重 ), ReduceLROnPlateau( monitor='val_loss', factor=0.5, # 学习率减半 patience=5, # 5轮不降再减 min_lr=1e-7 ) ] # 训练示例（需先准备X_train, y_train） # history = lstm_model.fit( # X_train, y_train, # validation_data=(X_val, y_val), # epochs=100, # batch_size=32, # callbacks=callbacks, # verbose=1 # )

逻辑说明：restore_best_weights=True是救命设置。LSTM 训练曲线常有震荡，最后几轮 loss 可能反弹，若不恢复最优权重，部署的模型就是次优的。ReduceLROnPlateau在 val_loss 卡住时自动降学习率，比固定学习率收敛更稳。

4.3 数据准备：用TimeseriesGenerator避免手动切片的索引灾难

手动用for i in range(len(data)-timesteps)切片极易出错：训练集、验证集、测试集的切片边界若没对齐，会导致数据泄露（用未来信息预测过去）。TimeseriesGenerator是 Keras 官方提供的安全切片器，它按时间顺序严格滑动，且支持shuffle=False（必须关！时间序列不能乱序）。

from tensorflow.keras.preprocessing.sequence import TimeseriesGenerator def prepare_lstm_data(residuals, timesteps, train_ratio=0.7, val_ratio=0.15): """ 用TimeseriesGenerator安全准备LSTM数据 residuals: SARIMA残差序列（pd.Series或np.array） timesteps: 输入窗口长度（由4.1节确定） """ # 转为numpy array并reshape data = np.array(residuals).reshape(-1, 1) # shape=(n_samples, 1) # 计算分割点 n_total = len(data) n_train = int(n_total * train_ratio) n_val = int(n_total * val_ratio) n_test = n_total - n_train - n_val # 分割（严格按时间顺序） train_data = data[:n_train] val_data = data[n_train:n_train+n_val] test_data = data[n_train+n_val:] # 创建Generator（注意：shuffle=False！） train_gen = TimeseriesGenerator( train_data, train_data, length=timesteps, sampling_rate=1, stride=1, batch_size=32, shuffle=False # ⚠️ 时间序列绝对禁止shuffle ) val_gen = TimeseriesGenerator( val_data, val_data, length=timesteps, sampling_rate=1, stride=1, batch_size=32, shuffle=False ) test_gen = TimeseriesGenerator( test_data, test_data, length=timesteps, sampling_rate=1, stride=1, batch_size=32, shuffle=False ) print(f"✅ 数据准备完成：训练{len(train_gen)}批，验证{len(val_gen)}批，测试{len(test_gen)}批") return train_gen, val_gen, test_gen # 使用示例 # train_gen, val_gen, test_gen = prepare_lstm_data(resid, timesteps=48)

关键点：TimeseriesGenerator的length=timesteps是输入长度，batch_size=32是每批样本数。它自动将train_data切成(timesteps, 1)的输入和(1, 1)的标签（即预测下一个点），无需手动X[i:i+timesteps],y[i+timesteps]，杜绝索引越界和错位。

5. 混合预测与避坑：SARIMA 趋势 + LSTM 残差修正的工程落地细节

混合预测不是sarima_forecast + lstm_forecast简单相加。SARIMA 输出的是带置信区间的趋势，LSTM 输出的是点估计残差修正。我们必须把 LSTM 的不确定性也纳入最终预测区间，否则业务方问“95% 置信区间多宽？”，你只能报错。此外，部署时的实时推理链路必须原子化，不能每次预测都重跑整个 SARIMA 拟合。

5.1 混合预测：用 SARIMA 的预测均值 + LSTM 的残差预测，再用 Bootstrap 估算总不确定性

SARIMA 的get_forecast()返回mean,se,conf_int，但它的se只反映统计模型不确定性。LSTM 的不确定性需单独评估。我们用Bootstrap 重采样：对 LSTM 的验证集预测误差（y_val - y_pred）抽样 1000 次，计算其 2.5% 和 97.5% 分位数，作为 LSTM 修正项的置信半宽。

import numpy as np from sklearn.utils import resample def hybrid_forecast(sarima_model, lstm_model, residuals, timesteps, steps_ahead=24): """ 混合预测主函数 sarima_model: 已拟合的SARIMA模型 lstm_model: 已训练的LSTM模型 residuals: SARIMA残差序列（用于Bootstrap） timesteps: LSTM输入窗口 steps_ahead: 预测步数（如24小时） """ # 步骤1：SARIMA预测（获取均值和标准误） sarima_fc = sarima_model.get_forecast(steps=steps_ahead) sarima_mean = sarima_fc.predicted_mean.values sarima_se = sarima_fc.se.mean() # 取平均标准误（简化） # 步骤2：准备LSTM输入（取残差序列最后timesteps点） last_resid_window = np.array(residuals[-timesteps:]).reshape(1, timesteps, 1) lstm_pred = lstm_model.predict(last_resid_window).flatten()[0] # 预测第一个点 # 步骤3：Bootstrap估算LSTM不确定性 # 用验证集残差预测误差分布（需提前保存val_pred和val_true） # 假设你有 val_errors = y_val_true - y_val_pred （长度N） # 这里用模拟数据演示 np.random.seed(42) # 模拟1000个LSTM预测误差（实际应来自验证集） val_errors = np.random.normal(loc=0, scale=0.15, size=500) # 均值0，标准差0.15 bootstrap_errors = np.array([ np.mean(resample(val_errors, n_samples=100)) for _ in range(1000) ]) lstm_lower = np.percentile(bootstrap_errors, 2.5) lstm_upper = np.percentile(bootstrap_errors, 97.5) # 步骤4：混合预测（SARIMA均值 + LSTM修正） # LSTM修正项：我们预测的是第一个点，但为保持一致性，对所有steps_ahead用同一修正？ # 实际工程中，应递归预测：用预测值更新残差窗口，再预测下一步 # 此处为简化，展示单步原理 hybrid_mean = sarima_mean[0] + lstm_pred # 总不确定性：SARIMA SE 和 LSTM 误差半宽的平方和开根 total_lower = hybrid_mean + (lstm_lower - 1.96 * sarima_se) total_upper = hybrid_mean + (lstm_upper + 1.96 * sarima_se) return { 'hybrid_mean': hybrid_mean, 'sarima_mean': sarima_mean[0], 'lstm_correction': lstm_pred, 'confidence_interval': (total_lower, total_upper) } # 使用示例（需先有val_errors） # result = hybrid_forecast(sarima_model, lstm_model, resid, timesteps=48, steps_ahead=1) # print(f"混合预测: {result['hybrid_mean']:.3f} (95% CI: [{result['confidence_interval'][0]:.3f}, {result['confidence_interval'][1]:.3f <p> <a href="https://download.csdn.net/download/kjm13182345320/88840753" style="color:#ec7500;font-size:14px;"> 本文还有配套的精品资源，点击获取 </a> <img alt="menu-r.4af5f7ec.gif" src="https://csdnimg.cn/release/wenkucmsfe/public/img/menu-r.4af5f7ec.gif" style="width:16px;margin-left:4px;vertical-align:text-bottom;cursor:text;"> </p>