简介基于MATLAB与QRLSTM分位数回归长短期记忆网络的时间序列区间预测程序包面向风速、负荷、功率等单变量预测场景适合需要从点预测进阶到区间预测的研究者与工程师。QRLSTM借助LSTM提炼序列长期依赖再利用分位数回归输出不同置信水平下的上下界有效刻画预测不确定性弥补传统点预测信息量不足的局限。资源共5个文件4个M脚本覆盖自定义分位数回归层、网络训练、预测输出与PICP/PIMWP等区间评价指标另附1个Excel数据文件输入输出均为单变量便于替换自有数据后直接运行运行环境为MATLAB 2018及以上版本。压缩包仅20KB轻量紧凑当前已有495人学习下载。读者可结合完整源码与数据系统掌握基于QRLSTM的时间序列区间预测建模流程、分位数损失函数设计及预测区间覆盖率、区间宽度等评估方法减少从零搭建的重复工作可直接用于风速、负荷、功率预测等科研或工程对比实验。1. 区间预测不是点预测加个误差条QRLSTM 到底改了什么做时间序列预测的人迟早会撞上一堵墙LSTM 输出的是一个点而业务方要的是一个区间。电价要上下限风速要置信带股票回撤要最坏情况——点预测给出均值或中位数但无法回答最差会到多少。传统的做法是在残差上假设正态分布加两个标准差当作区间可真实残差偏偏厚尾、偏态、异方差。QR-LSTM分位数回归长短期记忆神经网络换了个思路不预测期望直接预测条件分位数让网络自己学出不同概率水平下的取值边界。这样做出来的区间不依赖残差分布假设覆盖率和区间宽度都更可信。这篇博文从损失函数推导出发给出在 MATLAB 中从数据预处理、网络搭建到训练预测的完整路径代码可以直接复制改写。适合已经跑通 LSTM 点预测、想进一步做区间预测的工程师也适合在学术论文里需要对比区间预测方法的读者。整套方案用 MATLAB 原生深度学习工具箱实现不需要额外安装第三方库。2. 分位数回归与 LSTM 的组合逻辑为什么在 MATLAB 里用 QRLSTM2.1 分位数回归损失函数从对称误差到非对称误差普通回归最小化均方误差得到的是条件均值。分位数回归最小化的是 pinball loss分位数损失核心在于给正负误差分配不对称的权重。对分位点 $\tau \in (0,1)$单个样本的损失定义为L(y, y_pred, tau) (y - y_pred) * tau if y y_pred (y_pred - y) * (1 - tau) if y y_pred理解这个函数的关键在于高估和低估的代价不同。当预测值低于真实值时误差项带着权重 $\tau$反之带着 $1-\tau$。取 $\tau0.5$ 时两边权重相等最小化它得到的是条件中位数比 MSE 对异常值鲁棒得多。取 $\tau0.9$ 时低估的代价远高于高估优化器会偏向输出偏高的预测值从而逼近第 90 百分位。将这个损失函数接在 LSTM 输出层之后就得到 QRLSTM。网络结构本身与普通 LSTM 没有区别变化的只是损失函数和输出层的激活函数选择。输出层通常用线性激活或 identity因为分位数回归不允许输出被压缩到固定范围否则尾部行为会被截断。在 MATLAB 里pinball loss 不需要自己写循环。用自定义损失函数或直接在训练循环里计算均可。下面这段代码定义了一个可用的分位数损失函数function loss pinballLoss(Y, YPred, tau) % Y: 真实值大小 [numResponses, numTimeSteps, numObservations] % YPred: 预测值与 Y 同尺寸 % tau: 分位点标量 0 tau 1 err Y - YPred; loss mean(mean(max(tau * err, (tau - 1) * err), 2), 1); end这段逻辑里max(tau * err, (tau - 1) * err)用矢量化的方式实现了分段定义当err为正时取tau * err为负时取(tau - 1) * err。两个mean依次对时间步和样本求平均输出一个标量损失值。注意err的符号含义Y - YPred 0表示低估正是上文说的 $\tau$ 权重对应的情况。2.2 从三个分位点到完整区间覆盖率和区间宽度怎么权衡单次训练只输出一个分位数。要得到区间需要对多个 $\tau$ 值分别训练模型或者让网络输出多个值。$\tau0.1$ 和 $\tau0.9$ 的预测值围成的区间就是 80% 预测区间。区间质量看两个指标覆盖率真实值落在区间内的比例和平均区间宽度。这两个指标天然互相牵制——区间越宽覆盖率越高但预测也就越没信息量。这里有一个常见误用把多分位模型拆成多个独立 LSTM 训练。每个模型随机初始化不同收敛路径不同会导致 $\tau0.1$ 的输出偶尔高于 $\tau0.9$ 的输出即分位数穿越问题。解决办法有两个推荐第二个训练后排序预测完成后对输出做单调化排序。实现简单但破坏网络输出的连续性。共享隐藏层多头输出在 LSTM 层之上加一个全连接层输出维度设为分位点个数在损失函数中累加多个 $\tau$ 的 pinball loss。参数共享让不同分位数输出自然保持相关性穿越问题大幅减少。MATLAB 中第二种做法的损失计算代码如下function loss multiQuantileLoss(Y, YPred, taus) % Y: [numResponses, numTimeSteps, numObservations] % YPred: [numQuantiles, numTimeSteps, numObservations] % taus: 行向量例如 [0.1 0.5 0.9] numQ numel(taus); loss 0; for i 1:numQ tau taus(i); err Y - YPred(i, :, :); loss loss mean(mean(max(tau * err, (tau - 1) * err), 2), 1); end loss loss / numQ; end这个循环对分位点遍历每个 $\tau$ 计算独立的 pinball loss 后取平均。维度上 YPred 第一维是分位点MATLAB 的自动广播机制会把 Y 广播到每个分位点对应的通道上。这里的mean(mean(..., 2), 1)与上一节一致。这种做法比独立训练三个模型更稳定训练时间也短得多代价是最后一个全连接层要让每个分位点有独立的权重通道。3. MATLAB 里复现 QRLSTM 的最小完整流程3.1 数据构造先验证可行性再换自己的数据完整复现 QRLSTM 不需要一上来就处理复杂的真实业务数据。用一个带趋势和周期性的模拟时间序列验证代码链路比直接拿真实数据调试更高效。下面是构造示例数据的代码rng(42); T 1000; t (1:T); X sin(2*pi*t/50) 0.05*t 0.5*randn(T, 1);这段代码生成 1000 个时间步的合成数据包含周期项、线性趋势和高斯噪声。用它可以快速检验整个训练链路是否跑得通。数据集划分时用前 800 个时间步训练后 200 个测试。构造 LSTM 训练样本需要把时间序列转换成「输入窗口到输出」的监督学习格式。常用做法是用前numTimeSteps个连续值预测下一位numTimeSteps 20; XTrain zeros(numTimeSteps, 780, 1); YTrain zeros(1, 780, 1); for i 1:780 XTrain(:, i, 1) X(i:inumTimeSteps-1, 1); YTrain(:, i, 1) X(inumTimeSteps, 1); end这段代码里XTrain的维度是[时间步, 样本数, 特征数]对应 MATLAB 深度学习工具箱要求的序列格式。循环中每个样本取 20 个连续时间点作为特征第 21 个点作为目标值。得到 780 个样本是因为前 800 个时间步中最后一个完整窗口需要 20 个点做输入因此最多到第 780 个索引处取目标。测试集同理XTest zeros(numTimeSteps, 180, 1); YTest zeros(1, 180, 1); for i 1:180 XTest(:, i, 1) X(800i-1 : 800inumTimeSteps-2, 1); YTest(:, i, 1) X(800inumTimeSteps-1, 1); end这里测试起点选在 800样本从800i-1开始取与训练集的结构保持一致。注意YTest中第一个样本的真实值对应原始序列的第 821 个点后面的预测评估要以这个对齐关系为准忽略这个问题很容易在画图时把预测曲线平移一个窗口。3.2 网络定义LSTM 层后面接什么输出层QRLSTM 的网络主体与标准 LSTM 回归网络一致唯一区别在输出维度。以三分位点为例输出层需要 3 个神经元。用 MATLAB 的层数组定义numFeatures 1; numHiddenUnits 64; numQuantiles numel(taus); % taus 已在上面定义 layers [ sequenceInputLayer(numFeatures, Name, in) lstmLayer(numHiddenUnits, Name, lstm1) dropoutLayer(0.2, Name, drop) fullyConnectedLayer(numQuantiles, Name, fc_out) regressionLayer(Name, rout) ];sequenceInputLayer接受[特征数]参数对应XTrain第三维为 1。lstmLayer的numHiddenUnits是核心超参数64 个隐藏单元在中小规模数据上足够建模复杂依赖。dropoutLayer(0.2)放在 LSTM 之后全连接之前防止过拟合0.2表示随机丢弃 20% 的神经元输出。最终fullyConnectedLayer(numQuantiles)输出 3 个值。问题在于regressionLayer在 MATLAB 中内置的损失是均方误差不能直接用于分位数回归。因此这里不能用标准的trainNetwork一步训练而要用dlnetwork配合自定义训练循环dlnet dlnetwork(layers(1:end-1)); % 去掉 regressionLayerdlnetwork接受Layer数组但要求最后一层不能是输出层。layers(1:end-1)截断掉最后的regressionLayer保留从输入到全连接的完整计算图。后面的训练循环用modelGradients函数计算损失和梯度再调用adamupdate更新参数。3.3 训练循环手动迭代才能自定义损失训练函数的核心是数据转换成dlarray在前向传播后计算 pinball loss再用自动微分求梯度。完整的训练代码片段如下X dlarray(XTrain, CTB); Y dlarray(YTrain, CTB); numEpochs 200; miniBatchSize 64; learningRate 0.01; trailingAvg []; trailingAvgSq []; for epoch 1:numEpochs % 打乱数据顺序 idx randperm(size(X, 2)); X X(:, idx, :); Y Y(:, idx, :); for i 1:floor(size(X, 2)/miniBatchSize) batchIdx (i-1)*miniBatchSize1 : i*miniBatchSize; XBatch X(:, batchIdx, :); YBatch Y(:, batchIdx, :); % 计算梯度 [loss, gradients] dlfeval(modelGradients, dlnet, XBatch, YBatch, taus); % Adam 更新 [dlnet, trailingAvg, trailingAvgSq] adamupdate(dlnet, gradients, ... trailingAvg, trailingAvgSq, epoch, learningRate); end enddlarray(XTrain, CTB)中的格式标签说明C为通道维、T为时间维、B为批量维。dlfeval是 MATLAB 自动微分入口它会把modelGradients函数内部所有操作记录在计算图中。adamupdate实现了完整的 Adam 优化器trailingAvg和trailingAvgSq分别保存一阶和二阶动量初始为空数组时函数会自动初始化。modelGradients函数定义如下function [loss, gradients] modelGradients(dlnet, X, Y, taus) YPred forward(dlnet, X); loss multiQuantileLoss(Y, YPred, taus); gradients dlgradient(loss, dlnet.Learnables); endforward在前向传播时不累积梯度比predict更高效。multiQuantileLoss函数沿用 2.2 节的定义。关键在dlgradient(loss, dlnet.Learnables)它返回损失对所有可学习参数的梯度张量结构由 MATLAB 自动匹配。注意multiQuantileLoss内部的所有操作必须保持可微max和mean都是可微操作不需要额外处理。3.4 预测与区间组装三层境界训练完成后预测代码与普通 LSTM 没有区别dlYPred predict(dlnet, dlarray(XTest, CTB)); YPred extractdata(dlYPred); % YPred 尺寸: [3, 180, 1] lower squeeze(YPred(1, :, :)); median squeeze(YPred(2, :, :)); upper squeeze(YPred(3, :, :));extractdata把dlarray转回普通数值数组。矩阵重构后lower、median、upper是三个长度为 180 的列向量。画区间图时用fill函数figure; t 1:180; fill([t fliplr(t)], [lower fliplr(upper)], [0.85 0.9 1], EdgeColor, none); hold on; plot(t, YTest, k-, LineWidth, 1.5); plot(t, median, r--, LineWidth, 1);fill的第一个参数是 x 坐标序列[t fliplr(t)]把正向和反向的时间轴拼接构成多边形的回程路径。第二个参数对应 y 坐标前半段用下界、后半段用翻转的上界闭合出一个填充区域。颜色[0.85 0.9 1]是浅蓝色虚线红表示中位数预测黑色实线是真实值。这样一张图就能直观看出区间覆盖率与宽度。4. 训练 QRLSTM 的三个必调参数与常见坑4.1 学习率分位数损失对学习率更敏感配合学习率调度分位数损失的梯度没有 MSE 那样平缓的二次形式。在分位点接近 0 或 1 时pinball loss 在误差为负区间的斜率是 $1-\tau$接近 0梯度极小在误差为正区间斜率是 $\tau$很大。这种不对称梯度分布使得固定学习率容易在两个区域间震荡。实际调试中初始学习率建议设置在 0.001 到 0.01 之间配合piecewiseLearningRate调度learnRate 0.01; dropFactor 0.5; dropPeriod 50; if mod(epoch, dropPeriod) 0 learnRate learnRate * dropFactor; end每 50 轮学习率减半分位数损失在这种调度下比固定学习率收敛更快。调试时观察训练损失曲线若损失震荡剧烈且不下降优先把初始学习率调低十倍再试。4.2 分位点个数五个分位点通常够用九分位留作验证分位点选得越多区间边界更平滑能更精确地描述分布形状。但每个分位点对应输出层的一个神经元分位点太多会让全连接层参数量线性增长小数据集容易过拟合。常见做法是取奇数个对称分位点比如 [0.1 0.5 0.9] 或 [0.05 0.25 0.5 0.75 0.95]。选三个点时训练最快选五个点时区间边界更细可以画出不同置信水平的嵌套区间。增加分位点数量不需要修改网络结构之外的代码只需改taus向量。训练时间大致线性增长因为前向和反向传播的矩阵乘法维度变大了但 LSTM 本身的参数量不变所以增幅有限。经验上五个点对比三个点的训练时间增加约 20% 到 30%。4.3 时间窗口长度20 到 50 步是安全区间长了不一定好LSTM 理论上有长时间记忆能力但 QRLSTM 的分位数回归目标让网络把更多容量分配给分布尾部而不是序列模式记忆。窗口设得越长输入维度越高训练数据量需求越大尾部估计反而越不稳定。对日频或小时频数据窗口 20 到 50 步一般足够捕捉短期自相关。判断窗口是否合适看训练损失下降曲线如果验证集区间覆盖率在训练后期抖动剧烈可能是窗口过长导致过拟合。也可以做一个快速实验同一数据集分别用 10、20、40 步训练三个模型比较测试集平均区间宽度。宽度最小且覆盖率达标的窗口即为当前数据的最佳窗口。4.4 常见坑分位数穿越、NaN 梯度、标签错位分位数穿越问题在多头输出模型中依然可能发生。排查方法很简单预测完成后统计lower upper的比例如果超过 1%说明模型对尾部学习不充分。解决手段除了共享隐藏层外还可以在损失函数中加一个惩罚项强制分位数按序排列function loss orderedQuantileLoss(Y, YPred, taus) loss multiQuantileLoss(Y, YPred, taus); penalty 0; for i 1:numel(taus)-1 penalty penalty mean(mean(max(0, YPred(i,:,:) - YPred(i1,:,:)))); end loss loss 0.1 * penalty; endpenalty对相邻分位点的逆序差值取正部即只有逆序时才产生惩罚排序正确时惩罚为 0。系数0.1是推荐起点太大则会让网络牺牲损失精度来强制排序。NaN 梯度是训练循环常见问题通常由数据中包含NaN或Inf导致。在dlfeval抛出梯度计算错误时第一反应是检查原始数据是否干净而不是怀疑网络结构。用any(isnan(XTrain(:)))一行代码排除。5. 区间校准与一致性QRLSTM 上线前的最后一公里区间预测做出来只是第一步上线前要验证一件事预测区间在统计上是否可信。如果 90% 区间实际覆盖率只有 70%业务方据此做决策会高估风险。这引出一个叫区间校准的概念测试集上的经验覆盖率要接近名义覆盖率。经验覆盖率计算方式是真实值落在预测区间内的样本比例。结合 3.4 节的预测结果用如下代码计算coverage mean(YTest lower YTest upper, all); fprintf(名义覆盖率 80%%实际覆盖率 %.2f%%\n, coverage * 100);理想情况下输出实际覆盖率接近 80。如果偏差超过 5 个百分点优先检查训练轮数是否足够。分位数回归的收敛速度通常比点预测慢因为尾部区域的梯度信号稀疏。200 轮不够就加到 500 轮观察验证集覆盖率是否趋于稳定。覆盖率达标后的另一个诊断指标是 Pinball Score它同时惩罚过窄和过宽的区间是区间预测任务的标准评估指标。为了实现上的一致性验证对比两个 QRLSTM 变体或对比 QRLSTM 与 Bootstrap LSTM 区间方法Pinball Score 提供了可比较的量化基准。计算代码pinballScores zeros(numel(taus), size(YPred, 2)); for i 1:numel(taus) err YTest - YPred(i, :); pinballScores(i, :) max(taus(i) * err, (taus(i) - 1) * err); end meanPinball mean(pinballScores, all);这个分数结合了区间覆盖率和宽度两个维度的信息数值越低越好。切换数据集或调整参数时对比meanPinball比单独看覆盖率更稳健。例如区间特别宽时覆盖率肯定高但 Pinball Score 会恶化这能防止过度保守的预测。最后检查分位点之间的单调性用all(YPred(1,:) YPred(2,:) YPred(2,:) YPred(3,:))一键验证。若穿越比例不低4.4 节的排序惩罚项要在训练时加入而不是事后补救。经过覆盖率与 Pinball Score 双重验证的 QRLSTM才具备进入生产环境的资格。本文还有配套的精品资源点击获取