在 MATLAB 里做时间序列预测大家最先想到的往往是 ARIMA、BP、LSTM 这几类。但真正到工程现场尤其是环境里没装统计机器学习工具箱、plsregress根本调不出来的时候偏最小二乘算法PLS反而是非常实用的备选方案。最近我做一组库存时序数据预测就是在这种限制下被迫手写了一遍 PLS结果把那套“算法流程 回归系数换算 数据切分”想得比以往任何时候都清楚。这篇文章把我这次的完整过程整理出来PLS 做时间序列预测的核心逻辑、滞后特征怎么构造、不调用任何工具箱函数的 MATLAB 代码长什么样、以及我实际踩过哪几个坑。代码部分我会尽量逐行给注释方便直接抄走改数据就能用。1. PLS 为什么适合时间序列预测先搞懂它在算什么许多人对 PLS 的第一印象来自化学计量学标定光谱和浓度关系。但它的适用场景远不止这个。你把它用在时间序列预测上核心原因只有一个自变量之间高度相关时普通最小二乘会很脆弱而 PLS 通过提取潜变量在压缩维度的同时保留了与因变量最相关的信息。用滞后项做预测时这个优势尤其明显。假设你要用 t-1、t-2、……、t-10 这 10 个历史值去预测 t 时刻的值这些滞后项一定高度自相关尤其是相邻两阶之间的相关性往往极高。你用普通多元回归去拟合回归系数会被多重共线性搞得很不稳定换一段数据结果就跳得厉害。PLS 不是直接拿这 10 个变量做回归而是先构造若干“综合变量”每个综合变量都是原始滞后变量的线性组合并且保证这些综合变量与预测目标 y 的协方差最大。从算法层面看单因变量的 PLS也叫 PLS1核心迭代其实不复杂先初始化一个权重向量让 X 的列与 y 的相关性集中到这一方向上用这个权重算出得分向量 t它相当于从 X 中提取出的“第一根主脉”分别计算 X 和 y 在这个方向上的载荷从 X 和 y 中减去已经被提取出来的信息重复上述步骤直到提取出设定数量的潜变量。等到所有潜变量提取完再把结果换算回原始自变量空间得到一个可以直接用于预测的线性回归系数。这个过程不依赖任何高级函数只有矩阵转置、乘法、范数计算这些基础运算。我习惯用一个比方来理解它你可以把 X 看成一把乱掉的毛线直接扯其中一根线很容易拉断而且扯不出规律。PLS 先顺着最接近目标的方向“梳”出一缕把它抽走再梳下一缕最后得到的虽然是多缕毛线辫但它们每一缕都和目标更贴近整体预测更稳。对于时间序列来说这个特性让它天然适合成为一个基准模型。你可能最终还是要上 LSTM 或梯度提升树但在数据量不够大、需要快速出一个可解释预测结果时PLS 比神经网络更容易落地也比 OLS 稳得多。2. 时间序列预测用 PLS第一步是把时序变成回归样本PLS 是一个回归模型本身没有“时间记忆”。所以要想让它预测时间序列必须先做一个关键数据变换给序列构造滞后特征也叫时间窗口特征。举个例子如果我要用过去 7 天预测今天那么第 8 条样本的特征就是第 1 到第 7 天的值标签是第 8 天的实际值。第 9 条样本的特征则是第 2 到第 8 天的值标签是第 9 天的实际值。这样时间序列就被转换成了经典的有监督回归数据集。这个操作听着简单但有个很容易写错的细节特征行和标签行必须严格对齐。我下面给出的make_lagged函数按时间顺序构造不会出现错位问题。function [X, y] make_lagged(series, lag) % ------------------------------------------------------------------ % 输入: % series : n x 1 的时间序列数据 % lag : 滞后阶数例如 lag7 表示用前7期预测当期 % 输出: % X : (n-lag) x lag 的特征矩阵 % y : (n-lag) x 1 的标签向量 % ------------------------------------------------------------------ series series(:); % 保证是列向量 n length(series); X zeros(n-lag, lag); for k 1:lag X(:, k) series(lag1-k : n-k); end y series(lag1 : n); end这里最容易犯的错是索引偏移。比如第一行对应的是用第 1 到第 lag 期预测第 lag1 期那你填进去的值应该是series(1:lag)而不是其他区间。上面代码通过lag1-k : n-k这个索引切片保证了每一行的时间关系都是对齐的。你可以自己拿一个长度为 5 的序列试算一下把结果打印出来看比空想明白得多。构造完特征后下一步是切分数据集。这里必须强调时间序列不能随机切分。随机 K 折交叉验证会把未来信息混入训练集导致你的验证分数虚高得离谱。正确做法是严格按时间先后顺序切训练集用前面的 70%-80%验证集/测试集用最后一段。我给自己的切分方案是假设总样本 300 条滞后 8 阶那么特征矩阵有 292 行。取前 220 行为训练集后 72 行为测试集。注意测试集必须是不参与拟合、不参与组件选择的“新鲜数据”否则后续评估全是在自欺欺人。% 生成一段仿真时间序列好演示完整流程 rng(2024); n 300; t (1:n); series 1.2 0.03*t sin(0.06*t) randn(n,1)*0.4; lag 8; [X, y] make_lagged(series, lag); nTr 220; Xtr X(1:nTr, :); ytr y(1:nTr); Xte X(nTr1:end, :); yte y(nTr1:end);这组数据带了一个线性趋势和一个正弦周期噪声也加了一点比较接近现实中“有趋势、有波动、有噪声”的序列形态。接下来就可以正式写 PLS 拟合函数了。3. 不调用工具箱函数的 MATLAB 代码完整实现与逐行注释下面这套代码是核心。我没有用plsregress也没有用任何工具箱里的fit、cvpartition、predict之类的函数全部基于 MATLAB 基础矩阵运算。如果你之后要移植到 Octave 或者写成 Python numpy 版本基本可以一行行对应翻过去。3.1 PLS 拟合函数function model my_pls_fit(X, y, ncomp) % ------------------------------------------------------------------ % 手写偏最小二乘回归 PLS1单因变量版 % 不调用统计/机器学习工具箱函数 % % 输入: % X : N x p 特征矩阵 % y : N x 1 标签向量 % ncomp : 潜变量个数一般取 1~p 之间的整数 % % 输出: % model : 结构体保存回归系数、均值、偏置等信息 % ------------------------------------------------------------------ X X(:,:); y y(:); assert(size(X,1) length(y), X 与 y 的行数不一致); [N, p] size(X); ncomp max(1, min(ncomp, p)); % 防御性限制避免越界 % PLS 要求数据中心化。注意均值和方差必须在训练集上单独计算 meanX mean(X, 1); meanY mean(y, 1); Xc X - meanX; % 中心化后的 X yc y - meanY; % 中心化后的 y % 存储每次迭代的权重向量、载荷向量、y 的载荷标量 W zeros(p, ncomp); Pmat zeros(p, ncomp); C zeros(ncomp, 1); for k 1:ncomp % 第 1 步计算权重向量 w % Xc * yc 相当于 X 各列与 y 的协方差向量 % 归一化后作为潜变量方向。 w Xc * yc; w w / norm(w); % 第 2 步得分向量 t即 X 在方向 w 上的投影 t Xc * w; % 第 3 步X 的载荷向量 p % 它表示当前潜变量从 X 中“提取”出的特征模式 p Xc * t / (t * t); % 第 4 步y 的载荷系数 c % 表示当前潜变量对 y 的回归系数 c (t * yc) / (t * t); % 记录下来之后要换回原始空间 W(:, k) w; Pmat(:, k) p; C(k) c; % 第 5 步信息压缩 % 从 Xc 和 yc 中减去当前潜变量已经解释的部分 Xc Xc - t * p; yc yc - t * c; end % 潜变量空间系数换算到原始自变量空间 % 公式: beta W * (P * W)^(-1) * C % 不要直接对 Pmat 求逆用左除更稳定 beta W * ((Pmat * W) \ C); % 还原截距项 % 因为在中心化空间里做回归原始空间的截距需要补回来 bias meanY - meanX * beta; model.beta beta; model.bias bias; model.meanX meanX; model.meanY meanY; model.ncomp ncomp; end这里面的第 5 步是这个算法最核心的地方。它做的事情可以理解成“把已经解释过的信息从数据里抽走”剩下来的残差再去提取下一个潜变量。如果没有这一步后面提取到的方向会和前面高度重叠起不到降维和去相关的效果。3.2 预测函数拟合完得到的是回归系数预测就很简单了直接用矩阵乘function yhat my_pls_predict(model, Xnew) % ------------------------------------------------------------------ % 用训练好的手写 PLS 模型做预测 % % 输入: % model : my_pls_fit 返回的结构体 % Xnew : M x p 的新特征矩阵 % 输出: % yhat : M x 1 预测值 % ------------------------------------------------------------------ if size(Xnew, 2) ~ length(model.beta) error(新数据列数与训练时不一致请检查特征构造); end yhat Xnew * model.beta model.bias; end预测函数没有再做中心化因为偏置项bias已经在训练时把中心化造成的偏移吸收进去了。这一点很多人会弄混如果你在这里又减一次均值预测值就会整体偏移误差会突然变大。3.3 完整预测脚本把上面几个函数接起来在数据构造完成后可以直接跑这个流程% 训练 PLS 模型这里先手动指定潜变量个数为 4 model my_pls_fit(Xtr, ytr, 4); % 测试集预测 yhat my_pls_predict(model, Xte); % 计算评估指标 rmse sqrt(mean((yte - yhat).^2)); mae mean(abs(yte - yhat)); fprintf(潜变量数%d, RMSE%.4f, MAE%.4f\n, model.ncomp, rmse, mae); % 画图对比 figure; plot(yte, b-, LineWidth, 1.2); hold on; plot(yhat, r--, LineWidth, 1.2); legend(真实值, PLS预测值); xlabel(测试样本序号); ylabel(序列值); title(手写 PLS 时间序列预测结果); grid on;如果你的 MATLAB 版本支持脚本内局部函数R2016b 之后都支持可以直接把my_pls_fit、my_pls_predict、make_lagged三个函数放在同一个脚本末尾。如果是老版本就拆成三个.m文件保存文件名必须和函数名一致。4. 绕开工具箱之后怎么选潜变量个数和评估效果潜变量个数ncomp是手写 PLS 里最需要谨慎对待的超参数。它相当于你决定从原始特征里“梳出多少根辫子”。太少模型欠拟合趋势抓不住太多模型会把噪声也一起学进去测试集表现明显变差。我常用的方式是“预留一段尾部验证样本”。时间序列数据天然适合这么做训练集后面切一块连续的数据专门用来选参数。比如训练集 220 条再往后预留 30 条作为验证段选一个让验证误差最小的ncomp然后再回到全量训练集上以该参数重训模型。function bestL choose_ncomp(X, y, maxL) % ------------------------------------------------------------------ % 按时间序列方式选择 PLS 潜变量个数 % 不使用工具箱里的 cvpartition直接手写时序前向验证 % ------------------------------------------------------------------ N size(X, 1); valN max(10, round(N * 0.15)); % 预留末段15%作为验证 trN N - valN; XtrV X(1:trN, :); ytrV y(1:trN); XvaV X(trN1:end, :); yvaV y(trN1:end); bestE inf; bestL 1; for L 1:maxL m my_pls_fit(XtrV, ytrV, L); yv my_pls_predict(m, XvaV); e mean((yv - yvaV).^2); if e bestE bestE e; bestL L; end end end这里有一个很关键的原则用于选参的验证段和最后出指标的测试段必须是两段不同的数据。如果你用同一段数据既选参数、又报最终结果得到的是一个乐观偏差严重的结果看起来很好上线后立刻露馅。我的习惯是“训练-验证-测试”三段式训练拟合验证选参测试出最终指标。关于数据标准化我想额外说几句。如果你构造的特征全是同一个序列的不同滞后阶比如今天这个例子量纲一致那么做中心化就够了不必强行做标准化。但如果你把多个不同来源的特征拼进矩阵比如同时用了销售额、温度、库存量各自量纲差异很大那一定要在中心化的同时按列做标准差缩放否则 PLS 提取潜变量时会把量纲大的变量权重过度放大结果被它带偏。标准化要这样做在训练集上算均值mu和标准差sigma然后用同一组mu、sigma去处理测试集。绝对不能在全部数据上先算均值标准差再切分那会引入未来数据信息导致测试评估失真。5. 我在实操中踩过的三个坑排查过程与修正手写 PLS 的代码本身不难真正难的是数据构造和评估口径。下面这几个坑我几乎每一次给新场景写这套代码时都会碰上一两个。5.1 踩坑记录一中心化的“位置”错了几年前我第一次手写 PLS图省事直接把整个 X 和 y 在建模前一次性中心化然后才划分训练集和测试集。结果发现测试集预测曲线整体上移或下移RMSE 特别难看。后来排查了很久才意识到问题中心化时用了测试集自己的均值。这等于在训练之前就把测试集的“水平”透露给了模型听起来好像应该对结果有帮助才对但实际反而让训练集和测试集之间的均值偏移被错误地处理了。正确的做法是中心化参数只从训练集计算测试集预测时通过截距bias还原偏移。上面代码里meanX、meanY全部来自训练集预测函数只用Xnew * beta bias就是为了规避这个坑。5.2 踩坑记录二滞后矩阵索引错位但分数看起来不错有一个比中心化更隐蔽的问题滞后特征构造时索引写错。我有一版代码把第一行特征写成了series(1:lag)但标签却放到了series(lag1)看起来没问题实际第一行是用的“当场及之前”的数据预测未来并没有真正把预测点排除出去。这种泄漏会给出一个看似完美、实则虚假的预测曲线。后来我检查时写了一个小测试构造一个纯随机序列按道理任何模型都不该预测太好。如果 RMSE 显著低于序列标准差说明大概率有泄漏。这个排查技巧非常管用建议大家也保留。只要随机序列测试不过关就别急着调参先回头检查特征索引。5.3 踩坑记录三潜变量数过多把噪声学成了“规律”另一次我把ncomp直接设成和滞后阶数一样大想着“既然 PLS 能降维那就多取几个特征信息量更足”。结果训练集拟合漂亮测试误差反而比少取几个变量时还大。画出预测曲线后明显看到中后段跟着一些细小抖动在剧烈摆动这就是过拟合噪声的典型信号。当我改成用验证段试跑“潜变量数从 1 到 8”的误差变化后才发现这个数据集在 3 到 5 个潜变量附近误差最低再往上走验证误差不降反升。这个现象也说明潜变量个数不能拍脑袋必须放到数据上验证。我把这几个坑总结成下面这张表方便你对照自查问题现象可能原因正确处理预测曲线整体偏移中心化用了全样本均值只用训练集计算均值随机序列也能预测得很准滞后特征索引错位导致泄漏随机数据回归测试自检训练误差低、测试误差高ncomp 过大验证段选择潜变量个数测试指标虚高上线失效验证集与测试集混用三段式切分各司其职6. 把这段代码扩展成滚动预测增量更新的小技巧在很多实际业务里你需要的不是一次性预测一整段测试集而是每来一个新时刻就预测下一个时刻。这时可以用滚动预测模式保留一个长度固定的历史窗口每往前推一步就用最新真实值更新窗口重训一次模型或每隔几步重训一次。核心思路是先定义一个长度为window的训练窗口初始用第一批训练数据然后逐点预测window 180; bestL 4; Xhis Xtr(1:window, :); yhis ytr(1:window); nTest length(yte); yhatRoll zeros(nTest, 1); for i 1:nTest % 用当前窗口数据训练 m my_pls_fit(Xhis, yhis, bestL); % 预测下一个点 yhatRoll(i) my_pls_predict(m, Xte(i, :)); % 更新窗口加入真实观测踢掉最早一条 if i nTest Xhis [Xhis(2:end, :); Xte(i, :)]; yhis [yhis(2:end); yte(i)]; end end rmseRoll sqrt(mean((yte - yhatRoll).^2)); fprintf(滚动 PLS RMSE%.4f\n, rmseRoll);这个方案的好处是每次都让模型看到最近一段时间的动态趋势变化时响应更快。缺点是每步都要重新做一次 PLS 拟合数据多的时候计算量会上去。实际使用时我通常不会每一步都重训而是每 5 步或 10 步重训一次中间步骤直接用旧模型预测精度损失很小但时间能省下一大半。同样思路也适用于多步预测。如果一次要预测未来 3 步可以先预测第 1 步再把预测值作为特征拼到滞后窗口中滚动预测第 2 步、第 3 步。注意多步预测会把误差逐层累积步数越多结果越糙所以当预测步数较大时建议在每个预测时点都考虑用更保守的潜变量个数。这套代码我后来还移植到过 Python numpy 版本核心循环几乎原样保留只是把Xc * yc换成了Xc.T yc这种写法。手写一遍的价值就在这儿不是让你以后拒绝工具箱而是当工具箱不可用、或者你必须在另一种语言里复现算法时你脑子里已经有一张清晰的执行流程图不会慌。