尧图网站建设 的Matlab实现与交叉验证建模详解)
简介这组MATLAB源码用于实现偏最小二乘法PLS面向需要开展高维数据降维、回归建模与统计分析的研究人员和工程师尤其适合处理变量间存在多重共线性的数据场景。压缩包体积仅1024B包含1个m文件核心程序封装了从数据标准化、协方差矩阵计算、主成分提取、交叉验证主成分数到模型训练与预测的完整PLS流程代码量精炼且分段注释便于理解。已有219人学习参考可直接调用或改造后嵌入实际科研与工程建模任务。通过对该源码的研读能够直观掌握PLS在降维与回归中的核心思路熟悉迭代更新载荷向量、最大化潜在变量相关性的实现细节并复用于化学计量学、经济预测、工业过程监控等多元数据分析场景。1. 拿到“pls.rar”不代表会用偏最小二乘法从各路论坛上下载的“pls.rar_matlab源码”解压后通常有一堆pls.crossval、pls.nipals、pls.plot的函数很多人把数据一读套上demo就出曲线兴奋五分钟换一组数据立刻崩盘。这不是源码有bug而是偏最小二乘法PLS这把刀对刃口要求极高潜变量个数、数据预处理、交叉验证策略任何一个环节错模型都会变成对噪声的复读机。接下来不粉某个“最佳”源码包而是把PLS建模分析的原理和Matlab代码串起来让你无论用自带plsregress还是手写的pls.rar都能调出可信的结果。2. PLS的计算逻辑与预处理它凭什么能压住共线性数据2.1 从主成分回归到偏最小二乘建模思路转变在近红外光谱、化工过程监测这类数据里X的列数常常大于样本数列与列之间还有强共线性。直接做多元线性回归XX矩阵不可逆普通最小二乘的估计就不稳定。主成分回归先用PCA把X压缩成少数主成分再用这些主成分对y回归解决了秩亏问题但PCA分解X的时候并不看y前几个主成分可能只描述X自身的散射变化对y的预测几乎无贡献。偏最小二乘PLS在分解X的同时让得分向量的方向参考y让每个潜变量都尽量解释X和y的协方差所以同样的潜变量个数下PLS往往比PCR更能用在小样本建模中。如果要说得更技术一点PLS1单响应问题的目标可以写成找到单位方向向量w使t Xw与y的协方差最大。用样本数据表示就是最大化 Cov(t, y)同时在约束||w||1下。多响应时类似但需要同时衡量多个y。这个最优化问题没有闭合解所以常用NIPALS或SIMPLS迭代求解。SIMPLS用矩阵奇异值分解直接计算权重Matlab的plsregress内部实现更接近SIMPLS而很多开源的pls.rar包仍保留NIPALS写法。两者结果在标准化数据下基本一致差别体现在算法的数值稳定性上。方法是否使用y信息典型问题MLR否X严重共线时系数不稳定PCR否PCA阶段主成分可能和y不相关PLS是适合高维共线、小样本但需要调成分数2.2 NIPALS迭代潜变量提取过程的每一步NIPALS的核心思路是逐个提取潜变量每提取一个就消去X和Y中的这部分贡献。以多响应PLS为例第一个潜变量的计算流程大致如下。选Y的某一列作为u的初值通常取Y的第一列也可用随机列。计算X的权重 w Xu / ||Xu||这一步让w与当前u对齐。计算X的得分 t Xw。计算Y的载荷 q Yt / (tt)再更新u Yq。比较新旧t如果变化不超过阈值就继续否则回到第2步。计算X的载荷 p Xt / (tt)。更新残差X ← X - t pY ← Y - t q。对下一个潜变量重复整个过程。每一步都需要说清楚为什么第2步把X投影到u的方向上单位化是为了避免w无限扩大第4步的q是从Y中回回归出与t最匹配的载荷这样t携带的X信息不会与Y脱节。第7步的残差更新保证下一个潜变量专门解释之前没被解释的部分这也是PLS得分向量之间虽然正交但载荷不是正交的原因。单响应情况下y只有一列上述迭代在第1步就已经让u等于y所以等价于直接用w Xy不需要内部循环。2.3 均值中心化还是标准化预处理怎么选PLS实现里几乎都会有中心化这一步因为中心化之后模型可以写成 y Xb b0 的形式。对同量纲的连续变量均值中心化就够了例如近红外光谱里各波长的吸光度都在0到1附近直接中心化不会让某个波长因量纲过大而抢走权重。如果X里同时有温度、压力、流量这些量纲完全不同的变量就必须在中心化之后按每个变量的标准差缩放也就是标准化到单位方差否则计算w时量纲大的变量会天然占主导得到的潜变量会被“大数”绑架。这里出现一个常见错误整个数据集一起做标准化。正确的做法是先用训练集计算均值和标准差再用同一套统计量缩放验证集。更严格的说法是所有预处理的参数都必须在训练集内估计否则测试集的信息就污染了训练过程后面交叉验证的结果会虚高。光谱数据还经常用到SNV标准正态变量变换或一阶导数做平滑这些变换同样要先在训练集上拟合平滑窗口参数。场景预处理原因同量纲连续变量均值中心化保留变量尺度差很小不同量纲变量中心化标准化避免大尺度变量支配w光谱反射/吸光度再加上SNV或导数平滑消除散射和基线漂移分布偏斜很大log或Box-Cox后标准化让分布更接近对称3. 用Matlab手写NIPALS再读懂plsregress的返回参数3.1 一个可运行的单响应NIPALS实现很多包里的pls.m对多响应Y做了通用实现但调试时并不直观。我常用的办法是先写一版单响应PLS1跑通结果后再换成plsregress确认。PLS1的NIPALS比多响应简单得多因为y只有一列没有内部收敛循环代码不到40行。function [B, W, P, T] my_pls1(X, y, ncomp) % 单响应偏最小二乘回归的NIPALS实现 % 输入: X为n行p列预测矩阵, y为n行1列响应, ncomp为潜变量个数 % 输出: B为p行1列回归系数, W为p行ncomp权重, P为p行ncomp载荷, T为n行ncomp得分 % 注意: 调用前需要对X和y做中心化或标准化 [n, p] size(X); T zeros(n, ncomp); W zeros(p, ncomp); P zeros(p, ncomp); Q zeros(1, ncomp); Xk X; yk y; for i 1:ncomp w Xk * yk; % 协方差方向 w w / norm(w); % 单位化权重向量 t Xk * w; % 得分 q (t * yk) / (t * t); % y方向的载荷标量 p (Xk * t) / (t * t); % X方向的载荷向量 T(:, i) t; W(:, i) w; P(:, i) p; Q(1, i) q; Xk Xk - t * p; % 从X中消去该潜变量解释的部分 yk yk - t * q; % 从y中消去该潜变量解释的部分 end B W * ((P * W) \ Q); % 回归系数 end第5行的w计算本质上是当前X与y的协方差因为yk是残差所以这里的协方差排除了此前潜变量的影响。第8行q是把y对t做一元回归得到斜率这一步决定了该得分在y方向上的贡献。最后的B可以不在循环内累加而是等所有潜变量提取完后统一计算因为PLS的回归系数在X已经被收缩后无法逐分量简单地累加用P和W的外积公式更准确。实际使用中预测值为y_hat (new_x - mu_x) * B mu_yB只对应标准化后的空间。如果你在函数外做过标准化记得把截距和系数变换回原空间否则新手很容易在这里出错。3.2 对比plsregress的输出XL、YL、XS、YS、PCTVAR、MSEMatlab自带的plsregress是“官方做法”适合快速验证。它默认对X和Y做中心化但不会自动标准化。返回值比较多下面用spectra数据演示load spectra; X NIR; y octane; [XL, YL, XS, YS, BETA, PCTVAR, MSE] plsregress(X, y, 5); bar(cumsum(PCTVAR(1,:))); % 前5个潜变量解释X的比例返回矩阵里XL是X的载荷矩阵尺寸是p×ncompYL是Y的载荷尺寸是m×ncompXS是n×ncomp的X得分YS是n×m的Y得分。BETA是(p1)×m的回归系数矩阵注意第一行是截距。因此预测要用[ones(n,1), X] * BETA。PCTVAR是2×ncomp矩阵第一行X方差解释率第二行Y方差解释率第二行更值得关注因为Y的解释率才是建模目标。MSE包含训练集残差方差和交叉验证残差方差plsregress内部已经用10折交叉验证计算了MSE但折数固定如果想改成5折需要自己写。输出尺寸用途XLp×ncompX载荷观察变量重要性YLm×ncompY载荷响应与潜变量相关XSn×ncompX得分可用于样本聚类YSn×mY得分诊断离群点BETA(p1)×m回归系数含截距PCTVAR2×ncompX/Y方差解释百分比MSE2×(ncomp1)训练和交叉验证的MSE3.3 解压源码包后的路径检查与调用顺序很多网上下载的pls.rar解压后里面至少会有pls.m、cross_val.m、preprocess.m和demo.m。在跑demo之前先确认两件事。第一用which pls.m查看当前pls.m来自哪个路径避免Matlab工具箱里的同名函数被覆盖。第二把压缩包解压到不含中文和空格的纯英文路径否则addpath之后部分老代码会因为字符编码问题读不到数据文件。第三看看preprocess.m默认做了什么变换不少源码包会默认对X做标准化而y只做中心化。如果你的数据和作者当时测试的对象不同这个默认值可能就是“跑demo很顺、换数据就踩坑”的根源。4. 用交叉验证和RMSE把最佳潜变量数定在靠谱区间4.1 训练集R²高不代表模型好PCTVAR的陷阱当ncomp接近min(n,p)时训练集的R²会超过0.99PCTVAR第二行也接近1但这只是模型记住了X和y的对应关系。PLS不是没有过拟合只是它的潜变量比原始变量少过拟合来得晚一点。判断模型好坏要用没进过训练集的样本。最常用的指标是交叉验证的RMSE或Q²。Q²的定义是1减去预测残差平方和除以y的离差平方和实际计算时每个样本的预测来自不同折的模型和训练R²不是一回事。4.2 用K折交叉验证选择ncomp的Matlab模板我一般用10折交叉验证样本量小于100时会把折数改成5。下面这段直接跑会输出每个ncomp的CV RMSE。rng(2026); K 10; c cvpartition(size(X,1), KFold, K); ncomp_list 1:15; cv_rmse zeros(length(ncomp_list), 1); for j 1:length(ncomp_list) ncomp ncomp_list(j); rmse_sum 0; for k 1:K trIdx training(c, k); teIdx test(c, k); % 标准化统计量只从训练集计算 mu_x mean(X(trIdx,:), 1); s_x std(X(trIdx,:), 0, 1); mu_y mean(y(trIdx,:), 1); s_y std(y(trIdx,:), 0, 1); Xtr (X(trIdx,:) - mu_x) ./ s_x; ytr (y(trIdx,:) - mu_y) ./ s_y; Xte (X(teIdx,:) - mu_x) ./ s_x; yte (y(teIdx,:) - mu_y) ./ s_y; [XL, YL, XS, YS, BETA] plsregress(Xtr, ytr, ncomp); yhat_std [ones(size(Xte,1),1), Xte] * BETA; yhat yhat_std * s_y mu_y; % 还原到原始尺度 rmse_sum rmse_sum sqrt(mean((y(teIdx,:) - yhat).^2)); end cv_rmse(j) rmse_sum / K; end [best_rmse, idx] min(cv_rmse); fprintf(best ncomp %d, CV RMSE %.4f\n, ncomp_list(idx), best_rmse); plot(ncomp_list, cv_rmse, -o); xlabel(ncomp); ylabel(CV RMSE);cvpartition产生K折索引training/test返回逻辑索引。中心化y后预测再乘回s_y加回mu_y是为了让RMSE保持在原始量纲下。如果数据本身就是经标准化后的就不需要这一步。另外这段代码里用了std(...,0,1)0表示除以n-11表示按列写完整总比依赖默认值要稳妥。如果你还在老版本Matlab上跑减法和除法可以换成bsxfun结果一致。这段逻辑同样适用于光谱数据但前提是样本划分要合理见4.4。4.3 模型评价指标RMSE、R²、Q²怎么用不踩坑评估指标需要放在一起看。RMSE与y量纲一致只能作为绝对误差参考R²通常指测试集的R²计算公式是1 - 测试集残差平方和 / 测试集y的中心化平方和取值范围可以小于0说明模型比直接拿均值差。Q²特指交叉验证R²是模型稳定性的常用门槛。表格如下指标公式含义可接受范围注意RMSE预测残差平方根均值越小越好和y的量纲一致跨数据集比较无意义R²1 - SS_res/SS_tot越接近1越好如果为负说明预测不如直接平均Q²交叉验证版本R²大于0.5较好回归和分类标准不同不能死套0.5在化学计量学里R²0.9且Q²0.8常被认为是可用模型但这是针对近红外光谱这类信号较好的情况。如果数据本身噪声大R²在0.6附近也未必没用。所以要报告模型时至少要同时给出R²、RMSE、样本数和ncomp缺一个别人都无法复现评价过程。4.4 数据划分考虑光谱数据别用随机划分随机划分只适合样本独立同分布的数据。如果样本来自不同的批次、不同的时间随机划分会把时间漂移信息藏进训练集交叉验证结果会偏高。光谱数据常见的做法是按X的欧氏距离选择覆盖空间范围的样本比如KS算法先用样本向量间的距离选出距离最远的两个样本再逐步选择到已选样本最近距离最大的点。很多老源码包里会带着ks.m没有的话手写一个也不难。具体操作顺序是先做重现性分组同一批样品优先放进同一折再做KS划分。如果只有一个批次的样本随机划分也可以但要固定随机种子。5. 用置换检验给PLS模型的交叉验证结果补一个显著性证据5.1 置换检验的Matlab实现当样本量很小或Q²处于0.3这种“灰色地带”时单靠交叉验证很难判断模型是否真的存在可泛化的关系。置换检验的做法是保持X不变把y的标签随机打乱重复建模看看随机数据下能否达到原来的Q²。如果原始Q²出现在随机分布的高位说明模型与随机模式区分明显。下面是一个精简模板假定cv_q2函数就是上一节交叉验证代码的封装。rng(777); n_perm 200; perm_q2 zeros(n_perm, 1); for i 1:n_perm y_perm y(randperm(length(y)), :); perm_q2(i) cv_q2(X, y_perm, best_ncomp); % 返回10折CV的Q² end real_q2 cv_q2(X, y, best_ncomp); p_val (sum(perm_q2 real_q2) 1) / (n_perm 1); fprintf(原始Q²%.3f, 置换分布95%%分位%.3f, p%.3f\n, ... real_q2, quantile(perm_q2, 0.95), p_val);分子加1是习惯处理避免p值刚好为0也等价于把真实模型纳入排列分布累计。200次置换只够判断p是否小于0.05如果要报告p0.01至少跑999次。置换检验本质上是非参数检验对PLS模型分布不作正态假设因此它比单纯的F检验更适合高维共线模型。在论文或实验报告中写清楚“ncomp固定为交叉验证最优值预处理参数在每次置换中都重新从训练折内估计”代码才能被别人复现。这个细节也是区分“会调plsregress”和“懂偏最小二乘法”的常见分界线。本文还有配套的精品资源点击获取 偏最小二乘法(PLS)的Matlab实现与交叉验证建模详解
简介这组MATLAB源码用于实现偏最小二乘法PLS面向需要开展高维数据降维、回归建模与统计分析的研究人员和工程师尤其适合处理变量间存在多重共线性的数据场景。压缩包体积仅1024B包含1个m文件核心程序封装了从数据标准化、协方差矩阵计算、主成分提取、交叉验证主成分数到模型训练与预测的完整PLS流程代码量精炼且分段注释便于理解。已有219人学习参考可直接调用或改造后嵌入实际科研与工程建模任务。通过对该源码的研读能够直观掌握PLS在降维与回归中的核心思路熟悉迭代更新载荷向量、最大化潜在变量相关性的实现细节并复用于化学计量学、经济预测、工业过程监控等多元数据分析场景。1. 拿到“pls.rar”不代表会用偏最小二乘法从各路论坛上下载的“pls.rar_matlab源码”解压后通常有一堆pls.crossval、pls.nipals、pls.plot的函数很多人把数据一读套上demo就出曲线兴奋五分钟换一组数据立刻崩盘。这不是源码有bug而是偏最小二乘法PLS这把刀对刃口要求极高潜变量个数、数据预处理、交叉验证策略任何一个环节错模型都会变成对噪声的复读机。接下来不粉某个“最佳”源码包而是把PLS建模分析的原理和Matlab代码串起来让你无论用自带plsregress还是手写的pls.rar都能调出可信的结果。2. PLS的计算逻辑与预处理它凭什么能压住共线性数据2.1 从主成分回归到偏最小二乘建模思路转变在近红外光谱、化工过程监测这类数据里X的列数常常大于样本数列与列之间还有强共线性。直接做多元线性回归XX矩阵不可逆普通最小二乘的估计就不稳定。主成分回归先用PCA把X压缩成少数主成分再用这些主成分对y回归解决了秩亏问题但PCA分解X的时候并不看y前几个主成分可能只描述X自身的散射变化对y的预测几乎无贡献。偏最小二乘PLS在分解X的同时让得分向量的方向参考y让每个潜变量都尽量解释X和y的协方差所以同样的潜变量个数下PLS往往比PCR更能用在小样本建模中。如果要说得更技术一点PLS1单响应问题的目标可以写成找到单位方向向量w使t Xw与y的协方差最大。用样本数据表示就是最大化 Cov(t, y)同时在约束||w||1下。多响应时类似但需要同时衡量多个y。这个最优化问题没有闭合解所以常用NIPALS或SIMPLS迭代求解。SIMPLS用矩阵奇异值分解直接计算权重Matlab的plsregress内部实现更接近SIMPLS而很多开源的pls.rar包仍保留NIPALS写法。两者结果在标准化数据下基本一致差别体现在算法的数值稳定性上。方法是否使用y信息典型问题MLR否X严重共线时系数不稳定PCR否PCA阶段主成分可能和y不相关PLS是适合高维共线、小样本但需要调成分数2.2 NIPALS迭代潜变量提取过程的每一步NIPALS的核心思路是逐个提取潜变量每提取一个就消去X和Y中的这部分贡献。以多响应PLS为例第一个潜变量的计算流程大致如下。选Y的某一列作为u的初值通常取Y的第一列也可用随机列。计算X的权重 w Xu / ||Xu||这一步让w与当前u对齐。计算X的得分 t Xw。计算Y的载荷 q Yt / (tt)再更新u Yq。比较新旧t如果变化不超过阈值就继续否则回到第2步。计算X的载荷 p Xt / (tt)。更新残差X ← X - t pY ← Y - t q。对下一个潜变量重复整个过程。每一步都需要说清楚为什么第2步把X投影到u的方向上单位化是为了避免w无限扩大第4步的q是从Y中回回归出与t最匹配的载荷这样t携带的X信息不会与Y脱节。第7步的残差更新保证下一个潜变量专门解释之前没被解释的部分这也是PLS得分向量之间虽然正交但载荷不是正交的原因。单响应情况下y只有一列上述迭代在第1步就已经让u等于y所以等价于直接用w Xy不需要内部循环。2.3 均值中心化还是标准化预处理怎么选PLS实现里几乎都会有中心化这一步因为中心化之后模型可以写成 y Xb b0 的形式。对同量纲的连续变量均值中心化就够了例如近红外光谱里各波长的吸光度都在0到1附近直接中心化不会让某个波长因量纲过大而抢走权重。如果X里同时有温度、压力、流量这些量纲完全不同的变量就必须在中心化之后按每个变量的标准差缩放也就是标准化到单位方差否则计算w时量纲大的变量会天然占主导得到的潜变量会被“大数”绑架。这里出现一个常见错误整个数据集一起做标准化。正确的做法是先用训练集计算均值和标准差再用同一套统计量缩放验证集。更严格的说法是所有预处理的参数都必须在训练集内估计否则测试集的信息就污染了训练过程后面交叉验证的结果会虚高。光谱数据还经常用到SNV标准正态变量变换或一阶导数做平滑这些变换同样要先在训练集上拟合平滑窗口参数。场景预处理原因同量纲连续变量均值中心化保留变量尺度差很小不同量纲变量中心化标准化避免大尺度变量支配w光谱反射/吸光度再加上SNV或导数平滑消除散射和基线漂移分布偏斜很大log或Box-Cox后标准化让分布更接近对称3. 用Matlab手写NIPALS再读懂plsregress的返回参数3.1 一个可运行的单响应NIPALS实现很多包里的pls.m对多响应Y做了通用实现但调试时并不直观。我常用的办法是先写一版单响应PLS1跑通结果后再换成plsregress确认。PLS1的NIPALS比多响应简单得多因为y只有一列没有内部收敛循环代码不到40行。function [B, W, P, T] my_pls1(X, y, ncomp) % 单响应偏最小二乘回归的NIPALS实现 % 输入: X为n行p列预测矩阵, y为n行1列响应, ncomp为潜变量个数 % 输出: B为p行1列回归系数, W为p行ncomp权重, P为p行ncomp载荷, T为n行ncomp得分 % 注意: 调用前需要对X和y做中心化或标准化 [n, p] size(X); T zeros(n, ncomp); W zeros(p, ncomp); P zeros(p, ncomp); Q zeros(1, ncomp); Xk X; yk y; for i 1:ncomp w Xk * yk; % 协方差方向 w w / norm(w); % 单位化权重向量 t Xk * w; % 得分 q (t * yk) / (t * t); % y方向的载荷标量 p (Xk * t) / (t * t); % X方向的载荷向量 T(:, i) t; W(:, i) w; P(:, i) p; Q(1, i) q; Xk Xk - t * p; % 从X中消去该潜变量解释的部分 yk yk - t * q; % 从y中消去该潜变量解释的部分 end B W * ((P * W) \ Q); % 回归系数 end第5行的w计算本质上是当前X与y的协方差因为yk是残差所以这里的协方差排除了此前潜变量的影响。第8行q是把y对t做一元回归得到斜率这一步决定了该得分在y方向上的贡献。最后的B可以不在循环内累加而是等所有潜变量提取完后统一计算因为PLS的回归系数在X已经被收缩后无法逐分量简单地累加用P和W的外积公式更准确。实际使用中预测值为y_hat (new_x - mu_x) * B mu_yB只对应标准化后的空间。如果你在函数外做过标准化记得把截距和系数变换回原空间否则新手很容易在这里出错。3.2 对比plsregress的输出XL、YL、XS、YS、PCTVAR、MSEMatlab自带的plsregress是“官方做法”适合快速验证。它默认对X和Y做中心化但不会自动标准化。返回值比较多下面用spectra数据演示load spectra; X NIR; y octane; [XL, YL, XS, YS, BETA, PCTVAR, MSE] plsregress(X, y, 5); bar(cumsum(PCTVAR(1,:))); % 前5个潜变量解释X的比例返回矩阵里XL是X的载荷矩阵尺寸是p×ncompYL是Y的载荷尺寸是m×ncompXS是n×ncomp的X得分YS是n×m的Y得分。BETA是(p1)×m的回归系数矩阵注意第一行是截距。因此预测要用[ones(n,1), X] * BETA。PCTVAR是2×ncomp矩阵第一行X方差解释率第二行Y方差解释率第二行更值得关注因为Y的解释率才是建模目标。MSE包含训练集残差方差和交叉验证残差方差plsregress内部已经用10折交叉验证计算了MSE但折数固定如果想改成5折需要自己写。输出尺寸用途XLp×ncompX载荷观察变量重要性YLm×ncompY载荷响应与潜变量相关XSn×ncompX得分可用于样本聚类YSn×mY得分诊断离群点BETA(p1)×m回归系数含截距PCTVAR2×ncompX/Y方差解释百分比MSE2×(ncomp1)训练和交叉验证的MSE3.3 解压源码包后的路径检查与调用顺序很多网上下载的pls.rar解压后里面至少会有pls.m、cross_val.m、preprocess.m和demo.m。在跑demo之前先确认两件事。第一用which pls.m查看当前pls.m来自哪个路径避免Matlab工具箱里的同名函数被覆盖。第二把压缩包解压到不含中文和空格的纯英文路径否则addpath之后部分老代码会因为字符编码问题读不到数据文件。第三看看preprocess.m默认做了什么变换不少源码包会默认对X做标准化而y只做中心化。如果你的数据和作者当时测试的对象不同这个默认值可能就是“跑demo很顺、换数据就踩坑”的根源。4. 用交叉验证和RMSE把最佳潜变量数定在靠谱区间4.1 训练集R²高不代表模型好PCTVAR的陷阱当ncomp接近min(n,p)时训练集的R²会超过0.99PCTVAR第二行也接近1但这只是模型记住了X和y的对应关系。PLS不是没有过拟合只是它的潜变量比原始变量少过拟合来得晚一点。判断模型好坏要用没进过训练集的样本。最常用的指标是交叉验证的RMSE或Q²。Q²的定义是1减去预测残差平方和除以y的离差平方和实际计算时每个样本的预测来自不同折的模型和训练R²不是一回事。4.2 用K折交叉验证选择ncomp的Matlab模板我一般用10折交叉验证样本量小于100时会把折数改成5。下面这段直接跑会输出每个ncomp的CV RMSE。rng(2026); K 10; c cvpartition(size(X,1), KFold, K); ncomp_list 1:15; cv_rmse zeros(length(ncomp_list), 1); for j 1:length(ncomp_list) ncomp ncomp_list(j); rmse_sum 0; for k 1:K trIdx training(c, k); teIdx test(c, k); % 标准化统计量只从训练集计算 mu_x mean(X(trIdx,:), 1); s_x std(X(trIdx,:), 0, 1); mu_y mean(y(trIdx,:), 1); s_y std(y(trIdx,:), 0, 1); Xtr (X(trIdx,:) - mu_x) ./ s_x; ytr (y(trIdx,:) - mu_y) ./ s_y; Xte (X(teIdx,:) - mu_x) ./ s_x; yte (y(teIdx,:) - mu_y) ./ s_y; [XL, YL, XS, YS, BETA] plsregress(Xtr, ytr, ncomp); yhat_std [ones(size(Xte,1),1), Xte] * BETA; yhat yhat_std * s_y mu_y; % 还原到原始尺度 rmse_sum rmse_sum sqrt(mean((y(teIdx,:) - yhat).^2)); end cv_rmse(j) rmse_sum / K; end [best_rmse, idx] min(cv_rmse); fprintf(best ncomp %d, CV RMSE %.4f\n, ncomp_list(idx), best_rmse); plot(ncomp_list, cv_rmse, -o); xlabel(ncomp); ylabel(CV RMSE);cvpartition产生K折索引training/test返回逻辑索引。中心化y后预测再乘回s_y加回mu_y是为了让RMSE保持在原始量纲下。如果数据本身就是经标准化后的就不需要这一步。另外这段代码里用了std(...,0,1)0表示除以n-11表示按列写完整总比依赖默认值要稳妥。如果你还在老版本Matlab上跑减法和除法可以换成bsxfun结果一致。这段逻辑同样适用于光谱数据但前提是样本划分要合理见4.4。4.3 模型评价指标RMSE、R²、Q²怎么用不踩坑评估指标需要放在一起看。RMSE与y量纲一致只能作为绝对误差参考R²通常指测试集的R²计算公式是1 - 测试集残差平方和 / 测试集y的中心化平方和取值范围可以小于0说明模型比直接拿均值差。Q²特指交叉验证R²是模型稳定性的常用门槛。表格如下指标公式含义可接受范围注意RMSE预测残差平方根均值越小越好和y的量纲一致跨数据集比较无意义R²1 - SS_res/SS_tot越接近1越好如果为负说明预测不如直接平均Q²交叉验证版本R²大于0.5较好回归和分类标准不同不能死套0.5在化学计量学里R²0.9且Q²0.8常被认为是可用模型但这是针对近红外光谱这类信号较好的情况。如果数据本身噪声大R²在0.6附近也未必没用。所以要报告模型时至少要同时给出R²、RMSE、样本数和ncomp缺一个别人都无法复现评价过程。4.4 数据划分考虑光谱数据别用随机划分随机划分只适合样本独立同分布的数据。如果样本来自不同的批次、不同的时间随机划分会把时间漂移信息藏进训练集交叉验证结果会偏高。光谱数据常见的做法是按X的欧氏距离选择覆盖空间范围的样本比如KS算法先用样本向量间的距离选出距离最远的两个样本再逐步选择到已选样本最近距离最大的点。很多老源码包里会带着ks.m没有的话手写一个也不难。具体操作顺序是先做重现性分组同一批样品优先放进同一折再做KS划分。如果只有一个批次的样本随机划分也可以但要固定随机种子。5. 用置换检验给PLS模型的交叉验证结果补一个显著性证据5.1 置换检验的Matlab实现当样本量很小或Q²处于0.3这种“灰色地带”时单靠交叉验证很难判断模型是否真的存在可泛化的关系。置换检验的做法是保持X不变把y的标签随机打乱重复建模看看随机数据下能否达到原来的Q²。如果原始Q²出现在随机分布的高位说明模型与随机模式区分明显。下面是一个精简模板假定cv_q2函数就是上一节交叉验证代码的封装。rng(777); n_perm 200; perm_q2 zeros(n_perm, 1); for i 1:n_perm y_perm y(randperm(length(y)), :); perm_q2(i) cv_q2(X, y_perm, best_ncomp); % 返回10折CV的Q² end real_q2 cv_q2(X, y, best_ncomp); p_val (sum(perm_q2 real_q2) 1) / (n_perm 1); fprintf(原始Q²%.3f, 置换分布95%%分位%.3f, p%.3f\n, ... real_q2, quantile(perm_q2, 0.95), p_val);分子加1是习惯处理避免p值刚好为0也等价于把真实模型纳入排列分布累计。200次置换只够判断p是否小于0.05如果要报告p0.01至少跑999次。置换检验本质上是非参数检验对PLS模型分布不作正态假设因此它比单纯的F检验更适合高维共线模型。在论文或实验报告中写清楚“ncomp固定为交叉验证最优值预处理参数在每次置换中都重新从训练折内估计”代码才能被别人复现。这个细节也是区分“会调plsregress”和“懂偏最小二乘法”的常见分界线。本文还有配套的精品资源点击获取