QRBiLSTM分位数回归时序区间预测：MATLAB实战与避坑指南

简介：这份资源面向需要在MATLAB中实现时间序列区间预测的学习者与研究人员，聚焦风速、负荷、功率等场景，提供基于双向长短期记忆网络与分位数回归的QRBiLSTM完整实现方案。它解决的是传统点预测无法刻画不确定性的问题，通过分位数回归输出给定置信水平下的上下限，适合具备一定神经网络与MATLAB基础、希望深入掌握区间预测的中高级用户。压缩包共5个文件，包含4个m脚本与1个xlsx数据文件，整体约20KB，脚本分别承担主程序、分位数回归层、区间评价指标等核心功能，Excel数据便于直接替换与迁移实验。目前已有318人学习下载。读者可获得一套可直接运行的完整程序与配套数据，理解BiLSTM如何学习长期与短期依赖关系，并借助分位数回归层实现区间预测，同时通过PICP等指标评估区间质量，为股票、气象、交通等非线性时序研究提供可复用的代码框架与排错参考。

1. 从点预测到区间预测：QRBiLSTM 到底在解决什么问题

做过时间序列预测的人大多经历过这种尴尬：模型给出的单点预测值看着挺准，MAPE 也漂亮，但一放到实际决策场景里就露怯——电力负荷预测偏了 5%，调度计划就得推倒重来；风速预测少估了 2m/s，风机控制策略直接失效。点预测最大的问题是它不告诉你"这个数有多可信"。分位数回归（Quantile Regression, QR）就是冲着这个痛点来的：它不预测一个值，而是预测一组分位点，比如 0.05、0.5、0.95 三条线，中间那条是中位数预测，上下两条围成的区间就是预测区间。再配上双向长短期记忆网络（BiLSTM），让模型同时从前向和后向两个方向捕捉时序依赖，这就是 QRBiLSTM 的核心思路。

这套方案适合谁？做电力负荷、风速、光伏出力、交通流量、股价波动这类需要给出置信区间的时序预测任务的工程师和研究生。它不需要你从头推导分位数回归的数学证明，但需要你熟悉 MATLAB 的深度学习工具箱基本操作，知道怎么组织训练数据、怎么定义自定义损失函数。下面我从数据构造讲到模型训练再到区间评估，把整条链路拆开。

2. QRBiLSTM 的网络结构与分位数损失设计

2.1 为什么选 BiLSTM 而不是单向 LSTM

单向 LSTM 只利用当前时刻之前的信息，对时序预测来说这符合因果逻辑，但有个问题：它丢掉了未来上下文对当前状态的反向约束。BiLSTM 用两层 LSTM 分别跑正向和反向序列，然后把两个方向的隐状态拼接起来。注意，这里的"反向"不是让你用未来数据去预测未来，而是在训练阶段让模型学到更完整的序列表征。推理时你仍然只用历史窗口内的数据，不存在信息泄漏。

在 MATLAB 里搭建 BiLSTM 有两种方式。一种是直接用bilstmLayer，这是 R2021a 之后引入的封装层，用起来最省事。另一种是手动堆叠两个lstmLayer，一个设OutputMode='sequence'正向，另一个设Direction='reverse'，再拼接。我一般用第一种，代码简洁，不容易出错。

% 定义 BiLSTM 网络结构 numFeatures = 6; % 输入特征维度 numHiddenUnits = 128; % 隐藏单元数，常用 64~256 numQuantiles = 3; % 分位数个数，如 [0.05, 0.5, 0.95] layers = [ sequenceInputLayer(numFeatures, 'Name', 'input') bilstmLayer(numHiddenUnits, 'OutputMode', 'last', 'Name', 'bilstm') dropoutLayer(0.2, 'Name', 'dropout') % 防过拟合 fullyConnectedLayer(numQuantiles, 'Name', 'fc') ];

这段代码里几个参数值得说清楚。numHiddenUnits设 128 是个经验值，序列长度在 100~500 之间时够用；如果你的序列超过 1000 步，建议加到 256，但训练时间会明显变长。dropoutLayer的丢弃率设 0.2 是保守选择，数据量小于 5000 条时可以加到 0.3。fullyConnectedLayer的输出维度等于分位数个数，每个输出对应一个分位点的预测值。

2.2 分位数损失函数的 MATLAB 实现

分位数回归的关键在于损失函数。普通 MSE 只关心预测值和真实值的平方差，而分位数损失对高估和低估施加不对称惩罚。对于分位数 τ，损失函数定义为：

当 y_true ≥ y_pred 时，损失 = τ × (y_true - y_pred)；当 y_true < y_pred 时，损失 = (1-τ) × (y_pred - y_true)。

这个不对称性使得 τ=0.9 的模型倾向于高估，τ=0.1 的模型倾向于低估，从而自然形成区间。在 MATLAB 中，自定义损失需要返回损失值和梯度。下面是一个可直接用的实现：

function [loss, gradients] = quantileLoss(Y_pred, Y_true, quantiles) % Y_pred: numQuantiles x batchSize % Y_true: 1 x batchSize（广播到每行） % quantiles: 1 x numQuantiles，如 [0.05 0.5 0.95] numQ = size(Y_pred, 1); batchSize = size(Y_pred, 2); % 将真实值广播到每个分位数 Y_trueRep = repmat(Y_true, numQ, 1); % 计算残差 residual = Y_trueRep - Y_pred; % 分位数损失：不对称惩罚 tau = quantiles(:); % numQ x 1 tauRep = repmat(tau, 1, batchSize); lossMatrix = max(tauRep .* residual, (tauRep - 1) .* residual); loss = mean(lossMatrix, 'all'); % 梯度计算 gradients = zeros(size(Y_pred)); mask = residual >= 0; gradients(mask) = -tauRep(mask) / (numQ * batchSize); gradients(~mask) = (1 - tauRep(~mask)) / (numQ * batchSize); end

逻辑说明：max(tau*residual, (tau-1)*residual)这一行是整个损失的核心。当残差为正（真实值大于预测值，即低估）时，损失为 τ×残差；当残差为负（高估）时，损失为 (τ-1)×残差，由于 τ-1 是负数，乘以负残差得到正值。梯度部分按分位数分别计算，mask区分高估和低估两种情况。注意梯度要除以numQ * batchSize做归一化，否则学习率需要重新调。

参数方面，quantiles向量决定了你要预测哪些分位点。常见组合是[0.025, 0.5, 0.975]对应 95% 置信区间，或者[0.05, 0.25, 0.5, 0.75, 0.95]对应更细粒度的分布刻画。分位数越多，输出维度越大，训练难度也越高。我一般先用 3 个分位数跑通，确认损失收敛后再加。

2.3 训练选项与自定义训练循环

MATLAB 的trainNetwork默认用 MSE 损失，要换成自定义的分位数损失，需要用dlnetwork+ 自定义训练循环。这是很多人第一次做 QRBiLSTM 时卡住的地方——直接改trainingOptions里的损失函数是行不通的。

% 转换为 dlnetwork dlnet = dlnetwork(layers); % 训练超参数 numEpochs = 200; miniBatchSize = 64; learnRate = 0.001; decayRate = 0.95; decayPeriod = 20; % 优化器 averageGrad = []; averageSqGrad = []; gradDecay = 0.9; sqGradDecay = 0.999; % 训练循环 numIterations = floor(numTrainSamples / miniBatchSize); for epoch = 1:numEpochs shuffleIdx = randperm(numTrainSamples); for i = 1:numIterations idx = shuffleIdx((i-1)*miniBatchSize+1 : i*miniBatchSize); XBatch = XTrain(:,:,idx); YBatch = YTrain(idx); XBatch = dlarray(XBatch, 'CTB'); YBatch = dlarray(YBatch, 'CB'); [loss, grads] = dlfeval(@modelLoss, dlnet, XBatch, YBatch, quantiles); [dlnet, averageGrad, averageSqGrad] = adamupdate(dlnet, grads, ... averageGrad, averageSqGrad, epoch, learnRate, gradDecay, sqGradDecay); end % 学习率衰减 if mod(epoch, decayPeriod) == 0 learnRate = learnRate * decayRate; end end

modelLoss函数内部调用forward得到预测值，再调用前面写的quantileLoss。adamupdate是 MATLAB 内置的 Adam 优化器，gradDecay和sqGradDecay用默认值就行。学习率初始设 0.001，每 20 个 epoch 衰减 5%，这个策略在大多数时序任务上表现稳定。如果损失震荡厉害，把miniBatchSize加到 128 试试。

3. 数据构造与滑窗采样：把时序问题变成监督学习问题

3.1 滑窗构造的三种常见做法

时序预测的第一步是把原始序列切成 (输入窗口, 输出标签) 的样本对。假设原始序列长度为 N，输入窗口长度为 L，预测步长为 H，那么样本数为 N-L-H+1。MATLAB 里可以用循环构造，也可以用arrayDatastore配合自定义读取函数。我一般用循环，直观且方便调试。

function [X, Y] = createSlidingWindow(data, inputLen, predLen) % data: N x numFeatures 的归一化后数据 % inputLen: 输入窗口长度 % predLen: 预测步长（通常为1） N = size(data, 1); numFeatures = size(data, 2); numSamples = N - inputLen - predLen + 1; X = zeros(inputLen, numFeatures, numSamples); Y = zeros(1, numSamples); for i = 1:numSamples X(:,:,i) = data(i : i+inputLen-1, :); Y(i) = data(i+inputLen+predLen-1, 1); % 假设预测第1列 end end

inputLen的选择很关键。太短，模型看不到足够的上下文；太长，训练慢且容易过拟合。经验规则是至少覆盖两个完整周期。比如日周期数据，inputLen至少设 48（两天的半小时数据）；周周期数据，至少设 336。predLen通常设 1，做单步预测；如果要做多步，可以设 6 或 12，但输出维度要相应调整。

3.2 归一化与反归一化：别在这里翻车

归一化是时序预测里最容易出低级错误的地方。常见做法是把数据缩放到 [0,1] 或标准化到零均值单位方差。但有个坑：如果你用全量数据计算均值和标准差，然后划分训练集和测试集，就造成了信息泄漏。正确做法是只用训练集的统计量去归一化测试集。

% 正确做法：训练集统计量归一化测试集 [trainData, mu, sigma] = zscore(rawTrainData); testData = (rawTestData - mu) ./ sigma; % 反归一化预测结果 predOriginal = predNormalized .* sigma(1) + mu(1);

zscore返回的mu和sigma是训练集的均值和标准差，测试集直接用这两个值做变换。反归一化时注意只对预测目标列操作，别把其他特征也反了。这个细节看起来简单，但我在审阅别人代码时至少见过五次搞反的。

3.3 数据集划分与验证策略

时序数据不能随机打乱划分，必须按时间顺序切分。常见比例是 70% 训练、15% 验证、15% 测试。验证集用于早停和超参数选择，测试集只在最后评估时用一次。

numTrain = floor(0.7 * numSamples); numVal = floor(0.15 * numSamples); XTrain = X(:,:,1:numTrain); YTrain = Y(1:numTrain); XVal = X(:,:,numTrain+1:numTrain+numVal); YVal = Y(numTrain+1:numTrain+numVal); XTest = X(:,:,numTrain+numVal+1:end); YTest = Y(numTrain+numVal+1:end);

如果数据量本身就不大（少于 2000 条），可以考虑时间序列交叉验证：用前 60% 训练、接下来 10% 验证，然后窗口向前滑动。MATLAB 没有内置的时序交叉验证工具，需要自己写循环。我一般数据量够就不折腾这个，直接按比例切。

4. 训练过程中的避坑与排查清单

4.1 损失不下降或震荡剧烈

现象：训练 50 个 epoch 后损失仍在初始值附近波动，或者突然跳到 NaN。

原因通常有三个：学习率太大、梯度爆炸、数据里有异常值。先检查数据，用find(isnan(data))和find(isinf(data))扫一遍。如果有缺失值，用线性插值补上，别直接填零。学习率从 0.001 降到 0.0001 试试。梯度爆炸可以在modelLoss里加梯度裁剪：

gradients = dlupdate(@(g) max(min(g, 1), -1), gradients);

这行把梯度限制在 [-1, 1] 范围内，简单粗暴但有效。

4.2 分位数交叉：预测区间出现负宽度

现象：τ=0.05 的预测值大于 τ=0.95 的预测值，导致区间宽度为负。

这是分位数回归的经典问题。原因是各分位数独立训练，没有单调性约束。解决办法有两个：一是训练完后做后处理，对每个样本点强制排序；二是在损失函数里加惩罚项，当低分位数大于高分位数时施加额外惩罚。我一般用后处理，简单可靠：

% 后处理：强制分位数单调 predSorted = sort(predQuantiles, 1);

sort按行排序后，低分位数自然排在前面。这个操作不改变中位数预测，只调整上下界，对区间覆盖率影响很小。

4.3 区间覆盖率不达标

现象：设定 95% 置信区间，实际覆盖率只有 80% 左右。

原因可能是模型欠拟合，或者分位数损失权重需要调整。先检查训练集上的覆盖率，如果训练集也不达标，说明模型容量不够，加隐藏单元或加层。如果训练集达标但测试集不达标，说明过拟合，加 dropout 或减层。还有一种可能是数据分布偏移，测试集的统计特性和训练集差异大，这种情况需要重新审视数据划分策略。

4.4 训练速度慢到无法接受

现象：一个 epoch 要跑好几分钟，200 个 epoch 跑一整天。

先确认有没有用 GPU。MATLAB 里用gpuArray把数据和网络搬到 GPU 上：

if canUseGPU XBatch = gpuArray(XBatch); dlnet = dlupdate(@gpuArray, dlnet); end

如果没有 GPU，把miniBatchSize调大能提高 CPU 利用率，但内存占用也大。另外，bilstmLayer的numHiddenUnits从 128 降到 64 能省一半时间，精度损失通常在可接受范围内。

4.5 预测结果整体偏移

现象：预测曲线形状对，但整体偏高或偏低。

检查反归一化时用的mu和sigma是不是训练集的。如果测试集用了自己的统计量，就会产生系统性偏移。另一个可能是训练集和测试集的时间跨度差异大，比如训练集是淡季数据，测试集是旺季数据，模型没见过高负荷场景。这种情况要么扩充训练集覆盖更多工况，要么在损失函数里对高值区域加权。

5. 区间评估与进阶技巧：让预测区间真正可用

5.1 三个必须算的区间评估指标

区间预测不能只看点预测的 RMSE，必须评估区间的质量。我一般算三个指标：

PICP 的计算很简单：统计真实值落在上下界之间的比例。PINAW 是区间宽度的均值除以真实值的极差。CWC 把两者结合起来，覆盖率不达标时施加惩罚。

function [PICP, PINAW, CWC] = evaluateInterval(YTrue, YLower, YUpper, confidence) N = length(YTrue); inInterval = (YTrue >= YLower) & (YTrue <= YUpper); PICP = sum(inInterval) / N; rangeY = max(YTrue) - min(YTrue); PINAW = mean(YUpper - YLower) / rangeY; % CWC：覆盖率低于置信水平时惩罚 gamma = 10; % 惩罚系数 if PICP < confidence CWC = PINAW * (1 + gamma * exp(-gamma * (PICP - confidence))); else CWC = PINAW; end end

confidence设 0.95 对应 95% 区间。gamma控制惩罚力度，一般设 10~50。PICP 低于 0.9 时 CWC 会急剧增大，提醒你区间太窄了。

5.2 多分位数联合训练的一个实用技巧

如果你要预测 5 个以上的分位数，独立训练每个分位数的做法效率太低。更好的方式是在输出层一次性输出所有分位数，损失函数里对每个分位数分别计算再求和。这样共享底层 BiLSTM 特征，训练更快，分位数之间的一致性也更好。

% 多分位数输出层 numQuantiles = 5; quantiles = [0.05, 0.25, 0.5, 0.75, 0.95]; layers = [ sequenceInputLayer(numFeatures) bilstmLayer(128, 'OutputMode', 'last') dropoutLayer(0.2) fullyConnectedLayer(numQuantiles) ];

损失函数里对每个分位数分别算损失再取平均。注意fullyConnectedLayer的输出顺序要和quantiles向量的顺序对应，别搞混了。

5.3 用训练好的模型做滚动预测

实际部署时，你往往需要做多步滚动预测：用模型预测下一步，把预测值拼回输入窗口，再预测下一步。这个过程会累积误差，步数越多误差越大。缓解办法是每预测几步就用真实值校正一次，或者训练一个专门的多步预测模型。

% 滚动预测示例 numSteps = 24; % 预测未来24步 predSequence = zeros(1, numSteps); currentInput = XTest(:,:,end); % 最后一个测试样本的输入窗口 for step = 1:numSteps inputDL = dlarray(currentInput, 'CTB'); pred = predict(dlnet, inputDL); predValue = extractdata(pred(2, :)); % 取中位数分位数 predSequence(step) = predValue; % 更新输入窗口：去掉最早的时刻，拼入新预测 currentInput = [currentInput(2:end, :); ... [predValue, currentInput(end, 2:end)]]; end

这段代码里pred(2,:)取的是中位数分位数的输出，索引取决于你的分位数排列顺序。滚动预测的误差累积是固有缺陷，如果业务场景对长步预测精度要求高，建议直接训练一个输出维度等于预测步长的模型，而不是滚动。

5.4 一个我踩过的坑

最后说个血泪教训。我早期做 QRBiLSTM 时，把分位数损失里的梯度计算写错了符号，导致模型训练时损失确实在降，但预测出来的区间完全反了——高分位数预测值比低分位数还小。排查了两天才发现是gradients(~mask)那行少了个负号。所以自定义损失函数写完，一定要用数值梯度验证一遍：

% 数值梯度验证 epsilon = 1e-5; for i = 1:numel(Y_pred) Y_plus = Y_pred; Y_plus(i) = Y_plus(i) + epsilon; Y_minus = Y_pred; Y_minus(i) = Y_minus(i) - epsilon; lossPlus = quantileLoss(Y_plus, Y_true, quantiles); lossMinus = quantileLoss(Y_minus, Y_true, quantiles); numericalGrad = (lossPlus - lossMinus) / (2 * epsilon); % 和解析梯度对比，误差应在 1e-6 以内 end

这个验证步骤花不了十分钟，但能帮你省掉几天的调试时间。我现在养成了习惯：任何自定义损失函数，先跑数值梯度验证，通过了再上训练循环。希望这个习惯也能帮到你。

本文还有配套的精品资源，点击获取