简介面向时间序列区间预测的完整程序包基于双向长短期记忆网络并结合分位数回归构建QRBiLSTM模型适用于风速、负荷、功率等单变量序列数据。与传统点预测不同该模型能够给出指定置信水平下的预测上限与下限并以区间形式刻画未来值的不确定性。压缩包内共5个文件包括4个m源文件与1个Excel数据表源文件分别实现主程序、分位数回归层定义以及区间覆盖率、区间宽度等评价指标数据表提供示例序列可直接替换为自有数据。整体仅20KB运行环境为MATLAB 2018及以上输入输出均为单变量目前已有318人学习下载。读者通过该程序可以理解BiLSTM如何融合前向与后向信息、捕捉长期与短期依赖并掌握分位数回归输出预测区间的工作原理将数据替换后还可迁移至股票、气象、交通等非线性时序场景便于科研与工程人员快速实验、对比和二次开发。1. 区间预测不是把点预测加个误差带QRBiLSTM 在 MATLAB 里到底解决了什么问题做电力负荷预测那阵子我先拿到的是点预测模型领导看完补了一句误差你能给到几兆瓦我信但我更想知道什么范围外要启动备用机组。点预测给不出这个答案因为它只输出一个期望值把不确定性全吞掉了。时间序列区间预测要做的是输出一个条件分位数区间——比如未来 24 小时每小时的 5% 到 95% 分位带让调度员知道什么时刻波动大、什么时刻可以放心。QRBiLSTM 就是把双向长短期记忆网络和分位数回归组合起来BiLSTM 负责从时序上下文里提取特征分位数回归负责把特征映射成多个分位数输出最后在 MATLAB 里用自定义损失函数一次训练出整个区间。这篇文章写给要做负荷、风电、金融波动、设备退化这类区间预测的工程师从分位数损失的定义讲到滑窗构建、自定义回归层、训练调参和评估指标最后给出我踩过的四个真实坑。2. 分位数回归与 BiLSTM先搞清楚区间预测到底要让网络输出什么2.1 点预测的局限与分位数回归原理任何回归模型输出的都是一个条件期望 E[y|x]也就是在给定输入特征 x 时目标变量的平均值。时间序列里的真实过程往往带有异方差性负荷在早晚高峰波动大、凌晨波动小残差的方差是时变的。点预测把这种时变方差压缩成一个数字连“这个预测有多可信”都没有表达。分位数回归的思路完全不同它估计的是条件分位数 Q_τ(y|x)其中 τ ∈ (0,1)。τ0.5 时得到的是条件中位数τ0.05 和 τ0.95 构成的区间就是 90% 预测区间。它不对残差做正态性假设能捕捉非对称分布和厚尾这在电力负荷、风速、金融收益率上比均值回归稳得多。分位数回归的损失函数是 pinball loss也叫分位数损失。给定预测值 ŷ 和真实值 y单个样本在第 τ 分位上的损失为L_τ(y, ŷ) max(τ·(y−ŷ), (τ−1)·(y−ŷ))当 ŷ 低于 y 时损失权重是 τ当 ŷ 高于 y 时损失权重是 1−τ。这个不对称性很关键τ0.95 时模型把值预测低了要受重罚于是它会更倾向于给出略高的估计这个“略高”的程度由数据分布自适应决定。把这个损失在训练集上最小化得到的 ŷ 就是条件分位数的估计。τ0.5 时 pinball loss 退化为 MAE所以分位数回归是 MAE 回归的自然推广。在 MATLAB 里实现这个损失不难难的是 Deep Learning Toolbox 的默认回归层用的是均方误差不支持分位数权重。所以要跑 QRBiLSTM必须自定义一个回归层把 pinball loss 写进去。这个层的输入是网络输出预测分位数输出是损失值训练过程对它求梯度。2.2 BiLSTM 的双向读取机制与适用边界LSTM 靠门控机制解决长序列梯度消失问题但单向 LSTM 只能从过去读到现在信息流是单向的。BiLSTM 在结构上多了一条反向路径序列先正向过一遍 LSTM再反向过一遍然后把两个方向的隐藏状态拼接起来作为每个时间步的输出。对离线建模来说这等于让每个时刻的表示既包含历史上下文也包含“未来”上下文。序列内部的模式往往是对称的——负荷早高峰的形状、风速的爬坡段、传感器故障前的退化纹路——双向读取能把这种形状特征学得更完整。但这个结构有个边界严格意义的在线预测拿不到未来数据。训练时用 BiLSTM 没问题推理时如果还是把完整序列喂进去等于用未来信息预测当下这就是数据泄漏。我一般采取两个方案之一要么训练用 BiLSTM推理时只取正向路径的输出相当于把反向路径丢掉模型仍然受益于训练阶段的双向特征要么在滚动预测里给未来若干步填入占位值或前向估计值让反向路径有东西可读但这引入了额外误差。对“区间预测”这个场景多数论文和落地项目其实是离线评估——拿到一整段测试序列预测它的区间覆盖率这时候双向结构是安全的。如果部署要求是真正的逐点在线预测建议用单向 LSTM 或者把推理时的反向输入截断。MATLAB 里 bilstmLayer 封装了完整双向结构输入输出格式和 lstmLayer 基本一致没有额外的使用门槛。需要留意的是隐藏单元数的含义bilstmLayer(50) 表示正反向各 50 个单元拼接后送往全连接层的特征维度是 100比同参数的 LSTM 多一倍参数训练时间也相应上升。2.3 MATLAB 实现分位数损失的两条路线单模型多输出与三个独立模型把 pinball loss 接到 BiLSTM 上有两种常见做法。第一种是训练三个独立模型分别预测 5%、50%、95% 分位。每个模型结构一样只是损失函数里的 τ 不同。优点是不存在分位数交叉问题——三个模型互不干扰训练稳定排查起来很直观。缺点是训练成本乘以三而且三个模型各自收敛分位数之间的一致性全靠数据分布撑着。第二种是单模型多输出网络最后一层全连接输出三个节点分别对应 0.05、0.5、0.95 分位损失函数把三个 pinball loss 求和取平均。训练一次出三个分位数参数共享让特征提取部分只需要学一遍效率高输出还天然共享同一套隐含表示。我一般倾向第二种路线但需要处理分位数交叉问题模型可能在某些样本上输出 Q_0.95 Q_0.5这在业务上不成立。解决办法是在自定义损失里加一个排序惩罚项或者预测完之后做一次单调化处理。后面第 5 章会详细说。MATLAB 里第二种路线要做的就是自定义回归层属性里存一个分位数向量forwardLoss 里按 τ 逐列计算 pinball loss。这个层的写法直接决定训练能不能收敛、梯度能不能回传所以第 3 章把它单独拆出来逐行讲。3. 从原始序列到可训练模型滑窗构建、归一化与自定义损失层3.1 构建滑窗数据集用过去 L 步预测未来 h 步区间预测的监督学习样本需要把原始时间序列切片。核心思路是给定历史窗口长度 lag 和预测步长 horizon用 t−lag1 到 t 时刻的观测值预测 thorizon 时刻的值或值区间。窗口滑动一步就得到一个样本。lag 的选择直接影响特征质量我一般先用自相关图看数据的记忆长度负荷数据 lag 取 24 或 4824 小时周期金融日频数据取 5 到 20。horizon 决定预测目标的远近单步预测 horizon1多步可以先训单步模型再滚动也可以一次输出 h 个目标。QRBiLSTM 的常见做法是输出未来多个时刻的区间但为了控制复杂度本文先讲单步区间预测第 6 章再给滚动预测的扩展。下面这段 MATLAB 代码把一维序列转成训练用的输入输出。输入 X 是 cell 数组每个元素是 lag×1 的列向量输出 Y 是样本数×分位数个数 的矩阵。目标值是同一个复制成三列分别对应三个分位数。function [X, Y] makeSlidingWindow(data, lag, horizon, numQuantiles) % data: 一维列向量原始或归一化后的值 % lag: 历史窗口长度horizon: 预测未来第几步 % numQuantiles: 分位数个数如 3 n numel(data); nSamples n - lag - horizon 1; X cell(nSamples, 1); tgt zeros(nSamples, 1); for k 1:nSamples X{k} data(k : k lag - 1); % 过去 lag 步 tgt(k) data(k lag horizon - 1); % 未来第 horizon 步 end Y repmat(tgt, 1, numQuantiles); % 分位数回归的 target 复制多份 end参数说明lag 和 horizon 是超参数。lag 太小模型看不到周期性lag 太大样本数骤减且训练变慢。horizon 增大时预测难度非线性上升区间会明显变宽这是正常现象。numQuantiles 通常取 3低分位、中位数、高分位或 5额外加 0.1 和 0.9。repmat 把同一个 target 复制成多列配合自定义损失层逐列计算 pinball loss。这里不把 X 组装成矩阵而用 cell 数组是因为 sequenceInputLayer 要求每个样本是独立序列cell 数组每个元素对应一个序列样本。实际使用时要先做归一化再做滑窗还是先滑窗再归一化我建议先滑窗再归一化但要保证统计量只在训练集上计算。下面 3.2 节细说。3.2 归一化与数据划分统计量只在训练集上拟合测试集不能提前参与时间序列的自然顺序不能打乱这是和图像分类最大的一点不同。把序列随机 shuffle 会破坏时间依赖结构模型在测试集上表现会虚高因为测试样本里混着未来的上下文。我一般按时间顺序划分前 70% 训练、中间 15% 验证、最后 15% 测试。验证集用来做早停和调超参测试集只做最终评估。注意不能用 K 折交叉验证的默认随机拆分要用时间序列交叉验证前向链式划分否则同样是数据泄漏。归一化用 zscore先算训练集的均值和标准差再用同一组统计量去转换验证集和测试集。如果先在全序列上算均值再划分测试集的信息会通过统计量流进训练过程相当于把未来的尺度信息提前给了模型。这个错误很隐蔽因为训练误差不会变大测试误差反而显得更好一上真实部署就原形毕露。% 假设 rawData 是完整时间序列先按时间顺序切分 trainRaw rawData(1 : round(0.7*numel(rawData))); valRaw rawData(round(0.7*numel(rawData))1 : round(0.85*numel(rawData))); testRaw rawData(round(0.85*numel(rawData))1 : end); % 统计量只由训练集算出 mu mean(trainRaw); sigma std(trainRaw); trainNorm (trainRaw - mu) / sigma; valNorm (valRaw - mu) / sigma; testNorm (testRaw - mu) / sigma; % 滑窗并准备训练数据 lag 24; horizon 1; numQ 3; [XTr, YTr] makeSlidingWindow(trainNorm, lag, horizon, numQ); [XVa, YVa] makeSlidingWindow(valNorm, lag, horizon, numQ); [XTe, YTe] makeSlidingWindow(testNorm, lag, horizon, numQ);这里有个细节滑窗后训练集样本数不等于原始长度的 70%因为窗口会消耗 lag horizon − 1 个样本作为前缀。数据少时要留意这个边界损耗。推理阶段做反归一化时用 mu 和 sigma 把预测分位数还原到原始量纲。保存模型时把 mu 和 sigma 一起存成 .mat 文件否则换了机器推理会对不上。3.3 自定义分位数损失层把 pinball loss 接进 BiLSTM 的最后一环MATLAB 自定义回归层要继承 nnet.layer.RegressionLayer实现两个关键方法constructor 和 forwardLoss。forwardLoss 接收网络输出 Y样本数×输出节点数和目标 T样本数×输出节点数返回标量损失。下面的类定义把分位数向量作为属性传入逐列计算 pinball loss 后求平均。Y 和 T 的第 i 列对应第 i 个分位数因为前面 makeSlidingWindow 里用 repmat 把 target 复制到了每一列。classdef quantileLossLayer nnet.layer.RegressionLayer properties Taus % 分位数向量例如 [0.05 0.5 0.95] end methods function layer quantileLossLayer(taus) layer.Taus taus; layer.Name quantileLoss; layer.Description pinball loss for quantile regression; end function loss forwardLoss(layer, Y, T) % Y: numObservations x numQuantiles % T: numObservations x numQuantilestarget 各列相同 res T - Y; % 残差 q reshape(layer.Taus, 1, []); lossPerCol mean( max(q .* res, (q - 1) .* res), 1 ); loss mean(lossPerCol); % 各分位损失取平均 end end endforwardLoss 里的 max 实现了 pinball loss 的分段表达当残差为正预测偏低权重是 q残差为负预测偏高权重是 q−1。因为 q 是行向量res 是矩阵MATLAB 的隐式扩展会自动把 q 广播到每一列。lossPerCol 先对每个分位在所有样本上求平均再对所有分位求平均保证三个分位的损失在同等权重下参与梯度回传。这个层的梯度由 MATLAB 自动微分引擎计算不需要手写 backward。网络组装如下。sequenceInputLayer 接收 cell 数组序列bilstmLayer 后接 fullyConnectedLayer 输出三个节点最后接分位数损失层。这里输出节点数必须等于分位数个数不然 forwardLoss 里列对不上。numFeatures 1; % 单变量序列 numHidden 64; % 正反向各 64 numQuantiles 3; layers [ sequenceInputLayer(numFeatures, Name, in) bilstmLayer(numHidden, Name, bilstm) fullyConnectedLayer(numQuantiles, Name, fc) quantileLossLayer([0.05 0.5 0.95]) ]; % 训练选项 options trainingOptions(adam, ... MaxEpochs, 150, ... MiniBatchSize, 64, ... InitialLearnRate, 0.001, ... GradientThreshold, 5, ... ValidationData, {XVa, YVa}, ... ValidationFrequency, 20, ... Shuffle, never, ... % 时序数据不要打乱 Plots, training-progress); net trainNetwork(XTr, YTr, layers, options);Shuffle 参数设成 never 是我踩坑换来的。默认 trainingOptions 会在每个 epoch 打乱样本顺序对时序模型来说这等于破坏了序列间的先后关系——虽然每个样本内部的时间顺序没变但训练过程依赖样本间的分布一致性打乱后收敛会变得不稳定。BiLSTM 的双向路径如果训练时看到的是乱序样本反向路径学到的上下文是错位的。ValidationFrequency 设小一点能早点发现过拟合。GradientThreshold 设 5 是防止 LSTM 梯度爆炸的第一道防线后面第 5 章会讲 NaN 损失的处理。4. 训练与评估让 QRBiLSTM 真正收敛并用区间指标衡量好坏4.1 训练选项与超参数从 Adam、学习率到早停策略分位数损失的曲面不像 MSE 那样光滑它在残差符号切换处有个折点梯度方向会突然变化。用动量 SGD 在这个损失上训练容易震荡Adam 的自适应学习率对这类损失更友好。InitialLearnRate 我一般从 0.001 起步太高0.01 以上容易出现前几个 epoch 损失直接变 NaN太低则收敛过慢。MaxEpochs 设 150 到 300配合早停。MATLAB 内置的 trainingOptions 没有直接暴露早停耐心但可以用 ValidationData 配合 OutputFcn 或 Plots 观察曲线或者在代码里手动按 epoch 保存验证损失最小的模型。隐藏单元数这个超参数值得多说。bilstmLayer(64) 实际隐含维度是 128正反向拼接全连接层再映射到 3 个分位。隐藏单元太少比如 8 或 16序列的周期性和长程依赖学不进去区间会整体偏窄且覆盖率差太多比如 256 以上在样本量几千到几万时过拟合很快训练时间线性上涨。我的经验是先用 64 跑通再在 32 到 128 之间网格搜索用验证集 PICP 和 PINAW 的折中选最优。样本量只有几百时 32 足够样本量上万可以试 128。初始遗忘门偏置这项MATLAB 的 bilstmLayer 没有直接参数暴露但可以通过设置 InputWeightsInitializer 等初始化选项间接控制普通项目不必过度纠结用默认初始化就行。训练完成后模型保存用 save 命令把 net、mu、sigma、lag、horizon、分位数向量一起存进 .mat 文件。注意 MATLAB 训练好的网络对象在预测阶段要把输入同样做成 cell 数组格式和训练时一致否则会报维数不匹配。% 预测阶段拿测试集最后一个窗口举例 lastWin testNorm(end - lag 1 : end); % 列向量 pred predict(net, {lastWin}); % 输出 1x3对应三个分位数 pred mu sigma * pred; % 反归一化回原始量纲这里 predict 返回的是 cell 数组需要取 {1} 转成行向量。如果是多步滚动预测每滚动一步要把预测中位数接回窗口末尾再丢掉窗口第一个值。4.2 区间评估三板斧PICP、PINAW 与 CWC分位数回归模型的评估不能用 MSE 或 MAE因为优化目标是分位数损失用点预测指标选模型会选偏。评估区间预测有三个指标配合使用。PICPPrediction Interval Coverage Probability衡量区间覆盖率真实值落在预测区间内的比例越接近名义置信水平越好。但只看 PICP 不够——把区间拉到无穷宽覆盖率为 100%没有意义所以必须有 PINAWPrediction Interval Normalized Average Width控制区间宽度区间平均宽度除以测试集目标值的变化范围越窄越好。PICP 和 PINAW 是一对矛盾最终用 CWCCoverage Width Criterion做综合打分覆盖率低于目标时加指数惩罚高于目标时只看平均宽度。% 假设 Ytrue 是测试集真实值列向量Ypred 是 Nx3 预测矩阵已反归一化 % 列顺序与 Taus 一致[低分位, 中位数, 高分位] lower Ypred(:, 1); upper Ypred(:, 3); PICP mean(Ytrue lower Ytrue upper); PINAW mean(upper - lower) / (max(Ytrue) - min(Ytrue)); % 名义置信水平 0.9 mu_target 0.9; eta 5; % 覆盖率不足时的惩罚系数常用 5 或 10 CWC PINAW eta * exp(-(PICP - mu_target)) * (PICP mu_target);eta 的取值代表了业务上对覆盖率的重视程度。电力调度里漏报一个极值可能造成切负荷eta 取 10 以上如果只是趋势参考eta 取 2 到 5 即可。PICP 小于目标时 CWC 是指数惩罚一旦达到目标exp 项为零CWCPINAW这时模型之间比的就是谁区间更窄。这套评估逻辑应该直接写成一个脚本每次调参后自动输出三个指标不要只看训练曲线的损失值下降就收工。另外要画一张预测区间时序图横轴是时间纵轴是真实值和三条分位曲线。区间收窄的地方和真实值波动小的区域对应说明模型学到了异方差结构而不只是输出一个固定带宽的区间。如果区间宽度基本恒定说明模型没有捕捉到时变不确定性需要检查输入特征和信息流。5. QRBiLSTM 落地避坑分位数交叉、数据泄漏、NaN 与评估错用5.1 分位数交叉低分位预测值跑到高分位上面去了现象。训练完在测试集上预测某几个时间点的 5% 分位数比 95% 分位数还大区间看起来像麻花。这在单模型多输出结构中不少见尤其是训练数据噪声大、分位数值比较接近时。原因。三个分位数的输出是独立的全连接节点损失函数里没有约束 Q_0.05 ≤ Q_0.5 ≤ Q_0.95 的单调关系。网络拟合数据分布时如果某些样本的条件分布不对称梯度可能把两个节点的输出推成交叉状态。解决。三个手段按成本递增排序。第一预测后做单调化排序把三个输出按大小排序后直接当成单调分位数。这个方法简单但会改变原始预测值对覆盖率指标有一定影响。第二在自定义损失里加排序惩罚项算出相邻分位数的负差给负差加一个较大权重的平方惩罚强制网络学出单调输出。第三如果交叉频繁出现回到三个独立模型的方案彻底避免多输出之间的互相干扰。% 方式一事后单调化推理阶段 Ypred_sorted sort(Ypred, 2, ascend); % 此时列 1、2、3 分别对应低、中、高分位 % 方式二损失层里加单调惩罚 % 在 forwardLoss 中追加如下代码 diff12 Y(:, 2) - Y(:, 1); % 中位 - 低分位 diff23 Y(:, 3) - Y(:, 2); % 高分位 - 中位 penalty mean(max(0, -diff12).^2 max(0, -diff23).^2); loss loss 0.5 * penalty;penalty 的权重我一般先设 0.5观察训练损失和 PICP 的变化。权重大了会让网络牺牲原始分位数准确性来硬凑单调性权重小了压制不住交叉。交叉不频繁时事后排序最省事交叉超过 10% 的样本时建议直接换独立模型。5.2 数据泄漏归一化、序列乱序与窗口越界的三种翻车方式现象。训练过程的损失曲线正常下降验证集 PICP 高达 98%远超名义水平 90%但模型一部署到新数据上区间完全失灵覆盖率掉到 60%。原因。第一个泄露源是归一化用全序列统计量测试集的均值方差通过标准化流进了训练数据。第二个是训练时打开了 Shuffle验证集样本被混进训练批次或者验证集本身参与过训练。第三个是滑窗边界没算对测试集最后一个窗口的输入里包含测试集后面的时间点预测目标却更靠前序列出现“未来信息穿越”。解决。严格按第 3.2 节的三步走先切分再算统计量统计量只依赖训练集trainingOptions 里 Shuffle 设 never滑窗代码里确保 X 的最后一个元素时间下标小于等于目标值下标再往前一个 horizon。排查方法很简单把训练集和测试集的最大最小值打印出来如果测试集的归一化结果明显超出 [−3, 3]说明统计量很可能污染了。另一个检查手段是把训练好的网络在测试集上逐条预测时间对齐画出错误率随时间的分布如果某个时间段误差骤降大概率那个窗口里混了未来信息。5.3 训练不收敛与 NaN先查梯度再查输入现象。训练第几个 epoch 后损失变成 NaN训练曲线中断或者损失一直在高位震荡怎么也不降。原因。最常见的是学习率过大pinball loss 在残差符号切换处的梯度跳变被放大参数更新直接越过稳定区。其次是输入未归一化原始量纲数值动辄几百上千LSTM 内部的门控激活饱和梯度爆炸。第三个原因是样本量太少而隐藏单元太多网络在几十个样本上过拟合验证集损失不降反升看起来像不收敛。解决。GradientThreshold 设 5 到 10 是最有效的防线。学习率降到 0.0005 重新跑。检查训练集输入是否有 Inf 或 NaN这类脏数据会让梯度直接失效。如果隐藏单元 128 时不收敛换成 32 先验证数据和损失层正确再逐步加大。还有一个容易被忽视的细节自定义损失层的 forwardLoss 返回的是标量如果误返回成向量训练引擎会报维度错误或计算出 NaN。5.4 用 MSE 选模型评估目标与优化目标错位现象。用 MAE 或 RMSE 在验证集上选超参选出来的模型测试集 PICP 只有 70%而 PINAW 小得离谱——区间太窄覆盖面不够。原因。分位数回归优化的是分位数损失不是均方损失。MSE 对中间区域的样本给权重大对尾部样本权重小选出的模型会倾向于压缩区间宽度来降低平方误差牺牲掉极端值的覆盖。区间预测的评估逻辑必须是覆盖率优先、宽度其次。解决。所有超参选择都用 PICP、PINAW、CWC 三个指标综合判断。调参时先跑一组基线参数记录 PICP 和 PINAW然后每次只动一个超参对比 CWC 的变化。如果名义置信水平是 90%PICP 长期低于 80%优先增大隐藏单元数或把 lag 窗口加长而不是盲目调学习率。模型最终汇报时也只用这三个指标不带 MSE因为读者和业务方关心的就是“区间有没有罩住真实值、罩得紧不紧”。6. 进阶思路多步滚动预测与区间校准多步预测不能简单地把 horizon 设成 24 让网络一次吐 24 个分位数误差会累积且输出维度过大。我一般用滚动方式训练一个 horizon1 的 QRBiLSTM 模型预测完第一步把中位数预测值拼到窗口末尾删掉窗口第一个值继续预测第二步。滚动预测最重要的点是窗口里混入的预测值会引入误差传播所以每滚动一步要记录区间宽度变化——区间迅速变宽说明误差累积明显这时候应该考虑直接多步输出或者引入外部特征。区间校准是另一个值得投入的方向。QRBiLSTM 输出的区间本质上是条件分位数估计覆盖率可能系统性偏离名义水平。常见做法是用一个简单的 conformal 校准在验证集上计算残差真实值落在预测区间外的那些点的距离取经验分布的分位数作为校准偏移量预测时整体上下界向外扩或向内收这个量。这个后处理步骤十几行 MATLAB 代码就能完成能显著提升覆盖率稳定性。我今天能给你的最后一个建议是做区间预测不要一上来就追求复杂的网络结构。先用线性分位数回归跑通数据流和评估脚本再用 QRBiLSTM 替换特征提取部分你就能清楚地看到 BiLSTM 带来的增量到底在哪里。如果训练多次都过拟合回到 lag 窗口和隐藏单元这两个最敏感的超参上。这套流程我做了不下十个项目踩过的坑基本都写在这篇里了希望帮到你。本文还有配套的精品资源点击获取