简介多元线性回归是一种经典的多输入单输出回归方法该资源基于MATLAB实现完整预测流程面向计算机、电子信息工程、数学等专业学生适用于课程设计、期末大作业与毕业设计等学习场景。资源包共包含9个文件主要类型有3个MATLAB源码文件、3个结果图像文件、1个误差指标说明文本、1个CSV数据文件与1个MAT数据文件文件组织清晰压缩包仅247KB便于下载和快速浏览。代码采用参数化编程参数可灵活调整注释详细运行环境为MATLAB 2023b及以上版本运行后直接输出平均绝对误差、平均绝对百分比误差、均方误差、均方根误差和决定系数等多项评价指标可帮助读者系统走通数据预处理、模型训练、预测与误差分析完整流程。目前已有105人学习对刚接触回归预测的学生来说兼顾完整性与实用性是一份适合模仿与二次开发的入门实践资料。1. MLR多元线性回归到底是什么——一个多输入单输出的线性建模视角在工业预测和业务分析里多输入单输出回归几乎无处不在根据温度、压力、转速去预测设备剩余寿命根据访问量、加购数、促销标记去预测当日成交额。这类任务输入往往不止一个输出只有一个连续值MLR多元线性回归是解决这个问题最经典的起点。它把所有输入按不同权重相加再用一个截距项修正基准值得到一个可以直接写出公式的预测模型。和BP神经网络、随机森林相比MLR的精度上限可能不是最高但它参数少、训练快、结果完全可解释。当样本量只有几百甚至几十时MLR很少像复杂模型那样彻底失控这让它特别适合做基线模型和业务风控模型。接下来我按一份Matlab完整源码和数据的思路拆解从矩阵求解讲到数据预处理、共线性诊断最后给出一个交叉验证技巧帮你判断模型是真的可用还是只对训练集有效。2. MLR模型原理与最小二乘求解——为什么多输入单输出要这样建模2.1 从一元到多元MLR的数学表达与设计矩阵一元回归是 y b0 b1*x。到了多输入单输出场景假设有p个输入对第i个样本就有yi beta0 beta1xi1 beta2xi2 ... betap*xip eps_i把所有样本堆起来可以写成矩阵形式 Y X*beta epsilon。这里的X不是原始特征矩阵而是设计矩阵如果原始特征是X_raw那么设计矩阵要在第一列加上一列全1才能把截距项beta0放进系数向量。很多新手直接写下 X_raw \ y得到的是过原点的拟合当特征均值不为0时预测会出现系统性偏移。设计矩阵的维度是n行乘(p1)列每一行对应一个样本第j1列对应第j个特征的n个观测值。这个矩阵把回归问题转成了线性方程组求解。实际数据几乎不可能严格满足等式所以我们通过最小化残差平方和来求解这就是最小二乘的思想。2.2 最小二乘估计的闭式解与rank条件目标函数是 J(beta) ||Y - X*beta||_2^2。对beta求导并令导数为0可以得到正规方程(XX) * beta XY如果XX满秩解就是 beta_hat (XX)^(-1) * XY。在Matlab里不推荐直接写inv(X*X)*X*y更常见的做法是使用反斜杠 X\y。反斜杠会根据矩阵结构自动选择LU分解、QR分解或最小二乘算法数值上比显式求逆可靠得多。% 构造设计矩阵第一列全1对应截距项 % X_raw: n*py: n*1 X_design [ones(size(X_raw, 1), 1), X_raw]; % 方法1正规方程形式教学直观但会放大条件数 beta_normal (X_design * X_design) \ (X_design * y); % 方法2直接反斜杠推荐日常使用 beta X_design \ y;第一行在原始特征左侧拼接一列1保证截距项被纳入求解。方法1虽然写起来对称但XX会放大矩阵条件数遇到近似共线时数值结果不稳定。方法2用反斜杠直接对设计矩阵做QR分解稳定性更好后面所有代码都建议用这一种。当XX接近奇异时Matlab会输出警告这时需要回看特征之间的相关性而不是继续强行求解。2.3 与BP神经网络等非线性模型的边界MLR最核心的假设是输出与输入之间近似线性噪声独立同分布。一旦数据存在曲线关系、交互效应或者异方差MLR残差就会出现结构。很多人在这个阶段会考虑用pytorch实现bp神经网络回归预测与shap分析用深度模型拟合非线性再借助SHAP解释特征贡献。这没有问题但BP网络需要更大样本量、更仔细的调参而且多次训练结果会有随机波动。我一般会在动手前先画一组散点图把每个输入和输出放一起看趋势。如果关系近似直线MLR已经足够如果看到明显的弯折、分段或周期性再考虑加入二次项、交互项或者直接换树模型。下表列了几个维度方便理解MLR和其他模型的边界维度MLR多元线性回归BP神经网络树模型如随机森林可解释性系数直接读出影响方向和大小黑箱需额外做SHAP特征重要性可读样本量要求约为特征数的10倍即可起步通常需要几千条以上几百条到几千条训练时间矩阵求解秒级迭代训练可能分钟级构建多棵树中等非线性能力弱需手动构造特征变换强强主要风险共线性、线性假设不成立过拟合、调参成本高过拟合、解释不如系数直观选型时如果数据量不大业务上又需要向非技术方解释每个输入的影响MLR应该排在第一位。表格里的这些差异也决定了后面源码实现的侧重点先保证矩阵求解正确再处理好数据质量和评估。3. Matlab完整源码实现从数据到预测的核心代码3.1 数据集的组织readmatrix读取Excel或MAT格式一份可复现的Matlab源码第一步是搞清楚数据文件怎么组织。常见的Excel格式是每一行一个样本前p列是输入特征最后一列是输出目标。假设“数据.xlsx”里有100行、6列前5列是输入第6列是输出读取代码如下% 如果第一行是表头用NumHeaderLines跳过 data readmatrix(数据.xlsx, NumHeaderLines, 1); % 没有表头时可以直接 % data readmatrix(数据.xlsx); X_raw data(:, 1:end-1); % 所有行去掉最后一列 y data(:, end); % 输出列向量 [n, p] size(X_raw);readmatrix返回的是double矩阵Excel里的空单元格会被补成NaN。带着NaN进反斜杠结果会整片变成NaN所以读取之后要先单独检查。如果数据是.mat格式直接load(数据.mat)然后按变量名取出X_raw和y。这里要注意readmatrix的自动类型识别如果某一列里混入了文本整个矩阵会被读成cell或字符串此时需要先用double转换或把那一列删掉。3.2 最小二乘求解封装写一个mlr_fit函数为了后面多次调用我会把核心回归计算封装成一个函数mlr_fit放进mlr_fit.m文件里。函数同时返回系数、拟合值和评估指标训练集和测试集都能复用。function [beta, y_hat, stats] mlr_fit(X, y) % X: n*p 特征矩阵 % y: n*1 输出列向量 % beta: (p1)*1第一个元素是截距项 X_design [ones(size(X,1),1), X]; % 反斜杠求解最小二乘 beta X_design \ y; % 拟合值与残差 y_hat X_design * beta; resid y - y_hat; n size(X,1); p size(X,2); % 残差方差自由度为 n-p-1 sigma2 resid * resid / (n - p - 1); % 系数协方差矩阵共线性严重时会被放大 cov_beta sigma2 * inv(X_design * X_design); % 基本评估指标 SS_tot sum((y - mean(y)).^2); SS_res sum(resid.^2); R2 1 - SS_res / SS_tot; RMSE sqrt(SS_res / n); MAE mean(abs(resid)); stats struct(R2, R2, RMSE, RMSE, MAE, MAE, ... sigma2, sigma2, cov_beta, cov_beta, ... resid, resid, y_hat, y_hat); end函数首先构造设计矩阵再用反斜杠求系数。注意sigma2除以的是n-p-1目的是做自由度校正这样在样本量不大时残差方差的估计不会过于乐观。cov_beta这行不是必须的但如果你准备输出系数的t统计量或置信区间它就有用了。stats结构体里的R2是训练集上的决定系数只能说明拟合程度不能当作泛化精度。3.3 训练/测试集划分与预测输出模型评估必须在没见过的数据上进行所以要先把数据划成训练集和测试集。最常见的做法是随机打乱索引然后用前80%训练、后20%测试。rng(42); % 固定随机种子保证结果可复现 idx randperm(n); train_ratio 0.8; n_train floor(n * train_ratio); train_idx idx(1:n_train); test_idx idx(n_train1:end); % 训练 [beta_train, ~, ~] mlr_fit(X_raw(train_idx, :), y(train_idx)); % 测试集构造设计矩阵注意同样要加截距列 X_test [ones(length(test_idx),1), X_raw(test_idx, :)]; y_pred X_test * beta_train; % 测试误差 err y(test_idx) - y_pred; rmse_test sqrt(mean(err.^2)); fprintf(Test RMSE: %.4f\n, rmse_test);rng(42)让每次运行都得到同样的随机划分这对调试和复现至关重要。train_idx和test_idx是行索引数组直接作用于X_raw的行。测试集做预测时设计矩阵的构造必须和训练时完全一致也就是在最左边加一列1否则截距项会被丢掉。train_ratio设为0.8意味着80%的样本用于训练如果样本量只有50个测试集只剩10个评估波动会很大这时可以把比例降到0.7甚至0.6。注意如果数据是时间序列不能随机划分。过去的数据和未来的数据一旦混进同一个训练集预测相当于“偷看未来”结果会虚高。这时要按时间顺序前80%训练、后20%测试。命令/函数作用注意点readmatrix读取Excel/CSV数值数据有表头时设置NumHeaderLinesones构造截距项列必须在特征矩阵左侧拼接\最小二乘求解避免显式inv(XX)randperm随机划分样本先rng固定种子rmmissing删除含NaN行在划分前执行避免索引错位4. 数据预处理与模型评估决定MLR预测精度的几个关键点4.1 标准化与中心化zscore还是mapminmax在纯最小二乘MLR里系数会自适应特征尺度标准化不改变预测值。但一旦要做系数比较、VIF诊断或后续加岭回归量纲差异就会带来麻烦。常见做法是使用zscore标准化而不是mapminmax。mapminmax在神经网络代码里更常见它把数据映射到[-1,1]但如果测试集出现超出训练集范围的极值映射关系就会外推混乱。zscore用均值减除和标准差缩放对极值的敏感度相对更低。X_train X_raw(train_idx, :); X_test X_raw(test_idx, :); y_train y(train_idx); y_test y(test_idx); % 用训练集的均值/标准差做标准化 mu_X mean(X_train, 1); sigma_X std(X_train, 0, 1); X_train_std (X_train - mu_X) ./ sigma_X; X_test_std (X_test - mu_X) ./ sigma_X; % 输出变量也可以标准化 mu_y mean(y_train); sigma_y std(y_train); y_train_std (y_train - mu_y) / sigma_y;第5行的代码是最关键的一步测试集必须沿用训练集的mu_X和sigma_X不能在测试集上重新计算均值和标准差。如果重新计算测试集的信息就提前参与到模型流程中这叫数据泄漏会让评估结果偏乐观。std(X_train,0,1)中的0表示除以n-11表示沿着列的方向计算和corrcoef的默认行为一致。标准化之后训练出的beta就不再是原始量纲了预测时要反过来恢复到原始尺度。恢复代码如下% 标准化后的模型预测 X_test_std_design [ones(size(X_test_std,1),1), X_test_std]; y_pred_std X_test_std_design * beta_std; % 恢复到原始尺度 y_pred y_pred_std * sigma_y mu_y;这里有个很容易忽略的点如果对输出做了标准化训练时y_train_std使用mu_y和sigma_y那么测试集预测值要先乘sigma_y再加mu_y。只对输入标准化而对输出不做标准化系数会变得很难看而且截距项会包含两个不同量纲的信息。4.2 多重共线性诊断VIF阈值与处理方式MLR最怕的不是特征少而是特征之间高度相关。假设X中有两列几乎成比例XX就接近奇异系数估计的方差迅速变大可能得到正负号颠倒、数值巨大的权重。诊断多重共线性的常用指标是方差膨胀因子VIF。% 在标准化特征上计算相关系数矩阵 R corrcoef(X_train_std); % VIF是相关系数矩阵逆矩阵的对角线 VIF diag(inv(R)); for i 1:p fprintf(特征 %d 的VIF: %.2f\n, i, VIF(i)); endcorrcoef计算特征间的Pearson相关系数矩阵inv对R求逆diag取出对角线。VIF的含义是该特征被其他特征线性表示后剩下的独立信息还剩多少。VIF越大说明该特征的信息冗余越严重。VIF范围判断常见处理1无共线性无需处理1~5可接受正常建模5~10中等检查特征含义考虑删除或合并10严重优先删除相关性最高的特征之一或改用岭回归/主成分回归如果VIF显示严重共线性最直接的处理是删掉相关性最高的两个特征之一。删除后必须重新训练模型因为剩余特征的系数会重新分配。如果业务上每个特征都必须保留可以改用岭回归它通过给对角线加一个小常数来稳定求解。4.3 回归评估指标与残差模式判断测试集上的预测误差需要统一衡量标准。R2反映模型解释的方差比例RMSE是原量纲的根均方误差MAE对离群值相对不敏感MAPE则是百分比误差但在输出值接近0时会爆炸。% 测试集评估指标 R2_test 1 - sum((y_test - y_pred).^2) / sum((y_test - mean(y_test)).^2); RMSE_test sqrt(mean((y_test - y_pred).^2)); MAE_test mean(abs(y_test - y_pred)); MAPE_test mean(abs((y_test - y_pred)./y_test)) * 100; % 残差图 resid_test y_test - y_pred; figure; scatter(y_pred, resid_test, 15, filled); xlabel(预测值); ylabel(残差); grid on;残差图比指标数字更能说明问题。如果残差在零线附近随机散布MLR的线性假设基本成立如果残差随预测值增大呈喇叭形发散说明存在异方差如果残差呈现明显曲线说明漏掉了非线性项或交互项。遇到曲线时可以增加x1.*x2这种交互项但每增加一项就会消耗一个自由度样本量不大的时候反而加剧过拟合。MAPE在输出有负值或零值时没有意义这种情况建议直接用RMSE和MAE。R2高并不代表模型一定可靠尤其时间序列里两个共同趋势的变量可能产生虚假回归。所以最后还要看残差是否服从零均值正态分布这直接关系到预测区间是否可信。5. 验证模型有效性的一个实用技巧交叉验证与残差正态性检验5.1 为什么单次划分不够单次随机划分的测试RMSE依赖划分运气。有时候训练集恰好包含了所有极端值测试集又都落在均值附近RMSE会异常低反过来极端值进了测试集RMSE又会异常高。一次划分的结果不具备统计学稳定性所以完整源码里至少要做一个k折交叉验证。5.2 用cvpartition做五折交叉验证Matlab的cvpartition可以把样本划分为k份每轮轮流保留一份做测试。下面的代码输出每一折的RMSE以及均值与标准差rng(42); cv cvpartition(n, KFold, 5); rmse_fold zeros(cv.NumTestSets, 1); for k 1:cv.NumTestSets trIdx cv.training(k); teIdx cv.test(k); % 训练集设计矩阵 X_tr [ones(sum(trIdx), 1), X_raw(trIdx, :)]; beta_k X_tr \ y(trIdx); % 测试集预测 X_te [ones(sum(teIdx), 1), X_raw(teIdx, :)]; y_pred_k X_te * beta_k; rmse_fold(k) sqrt(mean((y(teIdx) - y_pred_k).^2)); end fprintf(5折CV RMSE: 均值 %.4f, 标准差 %.4f\n, ... mean(rmse_fold), std(rmse_fold));cv.training(k)和cv.test(k)返回逻辑索引sum(teIdx)统计测试样本数逻辑索引直接作用在X_raw行上。每一折训练和测试的数据都不重叠最终得到的均值就是模型泛化误差的更可靠估计。如果标准差接近均值数量级说明模型对数据划分很敏感要排查是否有异常样本或特征数量过少。如果后续要计算预测区间还需要检验残差正态性。可以用lillietest直接检查H lillietest(resid_test); % H0 表示不能拒绝残差服从正态分布的假设把交叉验证的RMSE均值和单次划分的测试RMSE放在一起对比二者接近说明模型性能可信单次划分远低于交叉验证均值说明这次划分偏乐观真正部署时要按均值预留误差余量。这一步是MLR回归项目交付前最值得做的验收动作。本文还有配套的精品资源点击获取 Matlab实现MLR多元线性回归：从最小二乘到交叉验证
简介多元线性回归是一种经典的多输入单输出回归方法该资源基于MATLAB实现完整预测流程面向计算机、电子信息工程、数学等专业学生适用于课程设计、期末大作业与毕业设计等学习场景。资源包共包含9个文件主要类型有3个MATLAB源码文件、3个结果图像文件、1个误差指标说明文本、1个CSV数据文件与1个MAT数据文件文件组织清晰压缩包仅247KB便于下载和快速浏览。代码采用参数化编程参数可灵活调整注释详细运行环境为MATLAB 2023b及以上版本运行后直接输出平均绝对误差、平均绝对百分比误差、均方误差、均方根误差和决定系数等多项评价指标可帮助读者系统走通数据预处理、模型训练、预测与误差分析完整流程。目前已有105人学习对刚接触回归预测的学生来说兼顾完整性与实用性是一份适合模仿与二次开发的入门实践资料。1. MLR多元线性回归到底是什么——一个多输入单输出的线性建模视角在工业预测和业务分析里多输入单输出回归几乎无处不在根据温度、压力、转速去预测设备剩余寿命根据访问量、加购数、促销标记去预测当日成交额。这类任务输入往往不止一个输出只有一个连续值MLR多元线性回归是解决这个问题最经典的起点。它把所有输入按不同权重相加再用一个截距项修正基准值得到一个可以直接写出公式的预测模型。和BP神经网络、随机森林相比MLR的精度上限可能不是最高但它参数少、训练快、结果完全可解释。当样本量只有几百甚至几十时MLR很少像复杂模型那样彻底失控这让它特别适合做基线模型和业务风控模型。接下来我按一份Matlab完整源码和数据的思路拆解从矩阵求解讲到数据预处理、共线性诊断最后给出一个交叉验证技巧帮你判断模型是真的可用还是只对训练集有效。2. MLR模型原理与最小二乘求解——为什么多输入单输出要这样建模2.1 从一元到多元MLR的数学表达与设计矩阵一元回归是 y b0 b1*x。到了多输入单输出场景假设有p个输入对第i个样本就有yi beta0 beta1xi1 beta2xi2 ... betap*xip eps_i把所有样本堆起来可以写成矩阵形式 Y X*beta epsilon。这里的X不是原始特征矩阵而是设计矩阵如果原始特征是X_raw那么设计矩阵要在第一列加上一列全1才能把截距项beta0放进系数向量。很多新手直接写下 X_raw \ y得到的是过原点的拟合当特征均值不为0时预测会出现系统性偏移。设计矩阵的维度是n行乘(p1)列每一行对应一个样本第j1列对应第j个特征的n个观测值。这个矩阵把回归问题转成了线性方程组求解。实际数据几乎不可能严格满足等式所以我们通过最小化残差平方和来求解这就是最小二乘的思想。2.2 最小二乘估计的闭式解与rank条件目标函数是 J(beta) ||Y - X*beta||_2^2。对beta求导并令导数为0可以得到正规方程(XX) * beta XY如果XX满秩解就是 beta_hat (XX)^(-1) * XY。在Matlab里不推荐直接写inv(X*X)*X*y更常见的做法是使用反斜杠 X\y。反斜杠会根据矩阵结构自动选择LU分解、QR分解或最小二乘算法数值上比显式求逆可靠得多。% 构造设计矩阵第一列全1对应截距项 % X_raw: n*py: n*1 X_design [ones(size(X_raw, 1), 1), X_raw]; % 方法1正规方程形式教学直观但会放大条件数 beta_normal (X_design * X_design) \ (X_design * y); % 方法2直接反斜杠推荐日常使用 beta X_design \ y;第一行在原始特征左侧拼接一列1保证截距项被纳入求解。方法1虽然写起来对称但XX会放大矩阵条件数遇到近似共线时数值结果不稳定。方法2用反斜杠直接对设计矩阵做QR分解稳定性更好后面所有代码都建议用这一种。当XX接近奇异时Matlab会输出警告这时需要回看特征之间的相关性而不是继续强行求解。2.3 与BP神经网络等非线性模型的边界MLR最核心的假设是输出与输入之间近似线性噪声独立同分布。一旦数据存在曲线关系、交互效应或者异方差MLR残差就会出现结构。很多人在这个阶段会考虑用pytorch实现bp神经网络回归预测与shap分析用深度模型拟合非线性再借助SHAP解释特征贡献。这没有问题但BP网络需要更大样本量、更仔细的调参而且多次训练结果会有随机波动。我一般会在动手前先画一组散点图把每个输入和输出放一起看趋势。如果关系近似直线MLR已经足够如果看到明显的弯折、分段或周期性再考虑加入二次项、交互项或者直接换树模型。下表列了几个维度方便理解MLR和其他模型的边界维度MLR多元线性回归BP神经网络树模型如随机森林可解释性系数直接读出影响方向和大小黑箱需额外做SHAP特征重要性可读样本量要求约为特征数的10倍即可起步通常需要几千条以上几百条到几千条训练时间矩阵求解秒级迭代训练可能分钟级构建多棵树中等非线性能力弱需手动构造特征变换强强主要风险共线性、线性假设不成立过拟合、调参成本高过拟合、解释不如系数直观选型时如果数据量不大业务上又需要向非技术方解释每个输入的影响MLR应该排在第一位。表格里的这些差异也决定了后面源码实现的侧重点先保证矩阵求解正确再处理好数据质量和评估。3. Matlab完整源码实现从数据到预测的核心代码3.1 数据集的组织readmatrix读取Excel或MAT格式一份可复现的Matlab源码第一步是搞清楚数据文件怎么组织。常见的Excel格式是每一行一个样本前p列是输入特征最后一列是输出目标。假设“数据.xlsx”里有100行、6列前5列是输入第6列是输出读取代码如下% 如果第一行是表头用NumHeaderLines跳过 data readmatrix(数据.xlsx, NumHeaderLines, 1); % 没有表头时可以直接 % data readmatrix(数据.xlsx); X_raw data(:, 1:end-1); % 所有行去掉最后一列 y data(:, end); % 输出列向量 [n, p] size(X_raw);readmatrix返回的是double矩阵Excel里的空单元格会被补成NaN。带着NaN进反斜杠结果会整片变成NaN所以读取之后要先单独检查。如果数据是.mat格式直接load(数据.mat)然后按变量名取出X_raw和y。这里要注意readmatrix的自动类型识别如果某一列里混入了文本整个矩阵会被读成cell或字符串此时需要先用double转换或把那一列删掉。3.2 最小二乘求解封装写一个mlr_fit函数为了后面多次调用我会把核心回归计算封装成一个函数mlr_fit放进mlr_fit.m文件里。函数同时返回系数、拟合值和评估指标训练集和测试集都能复用。function [beta, y_hat, stats] mlr_fit(X, y) % X: n*p 特征矩阵 % y: n*1 输出列向量 % beta: (p1)*1第一个元素是截距项 X_design [ones(size(X,1),1), X]; % 反斜杠求解最小二乘 beta X_design \ y; % 拟合值与残差 y_hat X_design * beta; resid y - y_hat; n size(X,1); p size(X,2); % 残差方差自由度为 n-p-1 sigma2 resid * resid / (n - p - 1); % 系数协方差矩阵共线性严重时会被放大 cov_beta sigma2 * inv(X_design * X_design); % 基本评估指标 SS_tot sum((y - mean(y)).^2); SS_res sum(resid.^2); R2 1 - SS_res / SS_tot; RMSE sqrt(SS_res / n); MAE mean(abs(resid)); stats struct(R2, R2, RMSE, RMSE, MAE, MAE, ... sigma2, sigma2, cov_beta, cov_beta, ... resid, resid, y_hat, y_hat); end函数首先构造设计矩阵再用反斜杠求系数。注意sigma2除以的是n-p-1目的是做自由度校正这样在样本量不大时残差方差的估计不会过于乐观。cov_beta这行不是必须的但如果你准备输出系数的t统计量或置信区间它就有用了。stats结构体里的R2是训练集上的决定系数只能说明拟合程度不能当作泛化精度。3.3 训练/测试集划分与预测输出模型评估必须在没见过的数据上进行所以要先把数据划成训练集和测试集。最常见的做法是随机打乱索引然后用前80%训练、后20%测试。rng(42); % 固定随机种子保证结果可复现 idx randperm(n); train_ratio 0.8; n_train floor(n * train_ratio); train_idx idx(1:n_train); test_idx idx(n_train1:end); % 训练 [beta_train, ~, ~] mlr_fit(X_raw(train_idx, :), y(train_idx)); % 测试集构造设计矩阵注意同样要加截距列 X_test [ones(length(test_idx),1), X_raw(test_idx, :)]; y_pred X_test * beta_train; % 测试误差 err y(test_idx) - y_pred; rmse_test sqrt(mean(err.^2)); fprintf(Test RMSE: %.4f\n, rmse_test);rng(42)让每次运行都得到同样的随机划分这对调试和复现至关重要。train_idx和test_idx是行索引数组直接作用于X_raw的行。测试集做预测时设计矩阵的构造必须和训练时完全一致也就是在最左边加一列1否则截距项会被丢掉。train_ratio设为0.8意味着80%的样本用于训练如果样本量只有50个测试集只剩10个评估波动会很大这时可以把比例降到0.7甚至0.6。注意如果数据是时间序列不能随机划分。过去的数据和未来的数据一旦混进同一个训练集预测相当于“偷看未来”结果会虚高。这时要按时间顺序前80%训练、后20%测试。命令/函数作用注意点readmatrix读取Excel/CSV数值数据有表头时设置NumHeaderLinesones构造截距项列必须在特征矩阵左侧拼接\最小二乘求解避免显式inv(XX)randperm随机划分样本先rng固定种子rmmissing删除含NaN行在划分前执行避免索引错位4. 数据预处理与模型评估决定MLR预测精度的几个关键点4.1 标准化与中心化zscore还是mapminmax在纯最小二乘MLR里系数会自适应特征尺度标准化不改变预测值。但一旦要做系数比较、VIF诊断或后续加岭回归量纲差异就会带来麻烦。常见做法是使用zscore标准化而不是mapminmax。mapminmax在神经网络代码里更常见它把数据映射到[-1,1]但如果测试集出现超出训练集范围的极值映射关系就会外推混乱。zscore用均值减除和标准差缩放对极值的敏感度相对更低。X_train X_raw(train_idx, :); X_test X_raw(test_idx, :); y_train y(train_idx); y_test y(test_idx); % 用训练集的均值/标准差做标准化 mu_X mean(X_train, 1); sigma_X std(X_train, 0, 1); X_train_std (X_train - mu_X) ./ sigma_X; X_test_std (X_test - mu_X) ./ sigma_X; % 输出变量也可以标准化 mu_y mean(y_train); sigma_y std(y_train); y_train_std (y_train - mu_y) / sigma_y;第5行的代码是最关键的一步测试集必须沿用训练集的mu_X和sigma_X不能在测试集上重新计算均值和标准差。如果重新计算测试集的信息就提前参与到模型流程中这叫数据泄漏会让评估结果偏乐观。std(X_train,0,1)中的0表示除以n-11表示沿着列的方向计算和corrcoef的默认行为一致。标准化之后训练出的beta就不再是原始量纲了预测时要反过来恢复到原始尺度。恢复代码如下% 标准化后的模型预测 X_test_std_design [ones(size(X_test_std,1),1), X_test_std]; y_pred_std X_test_std_design * beta_std; % 恢复到原始尺度 y_pred y_pred_std * sigma_y mu_y;这里有个很容易忽略的点如果对输出做了标准化训练时y_train_std使用mu_y和sigma_y那么测试集预测值要先乘sigma_y再加mu_y。只对输入标准化而对输出不做标准化系数会变得很难看而且截距项会包含两个不同量纲的信息。4.2 多重共线性诊断VIF阈值与处理方式MLR最怕的不是特征少而是特征之间高度相关。假设X中有两列几乎成比例XX就接近奇异系数估计的方差迅速变大可能得到正负号颠倒、数值巨大的权重。诊断多重共线性的常用指标是方差膨胀因子VIF。% 在标准化特征上计算相关系数矩阵 R corrcoef(X_train_std); % VIF是相关系数矩阵逆矩阵的对角线 VIF diag(inv(R)); for i 1:p fprintf(特征 %d 的VIF: %.2f\n, i, VIF(i)); endcorrcoef计算特征间的Pearson相关系数矩阵inv对R求逆diag取出对角线。VIF的含义是该特征被其他特征线性表示后剩下的独立信息还剩多少。VIF越大说明该特征的信息冗余越严重。VIF范围判断常见处理1无共线性无需处理1~5可接受正常建模5~10中等检查特征含义考虑删除或合并10严重优先删除相关性最高的特征之一或改用岭回归/主成分回归如果VIF显示严重共线性最直接的处理是删掉相关性最高的两个特征之一。删除后必须重新训练模型因为剩余特征的系数会重新分配。如果业务上每个特征都必须保留可以改用岭回归它通过给对角线加一个小常数来稳定求解。4.3 回归评估指标与残差模式判断测试集上的预测误差需要统一衡量标准。R2反映模型解释的方差比例RMSE是原量纲的根均方误差MAE对离群值相对不敏感MAPE则是百分比误差但在输出值接近0时会爆炸。% 测试集评估指标 R2_test 1 - sum((y_test - y_pred).^2) / sum((y_test - mean(y_test)).^2); RMSE_test sqrt(mean((y_test - y_pred).^2)); MAE_test mean(abs(y_test - y_pred)); MAPE_test mean(abs((y_test - y_pred)./y_test)) * 100; % 残差图 resid_test y_test - y_pred; figure; scatter(y_pred, resid_test, 15, filled); xlabel(预测值); ylabel(残差); grid on;残差图比指标数字更能说明问题。如果残差在零线附近随机散布MLR的线性假设基本成立如果残差随预测值增大呈喇叭形发散说明存在异方差如果残差呈现明显曲线说明漏掉了非线性项或交互项。遇到曲线时可以增加x1.*x2这种交互项但每增加一项就会消耗一个自由度样本量不大的时候反而加剧过拟合。MAPE在输出有负值或零值时没有意义这种情况建议直接用RMSE和MAE。R2高并不代表模型一定可靠尤其时间序列里两个共同趋势的变量可能产生虚假回归。所以最后还要看残差是否服从零均值正态分布这直接关系到预测区间是否可信。5. 验证模型有效性的一个实用技巧交叉验证与残差正态性检验5.1 为什么单次划分不够单次随机划分的测试RMSE依赖划分运气。有时候训练集恰好包含了所有极端值测试集又都落在均值附近RMSE会异常低反过来极端值进了测试集RMSE又会异常高。一次划分的结果不具备统计学稳定性所以完整源码里至少要做一个k折交叉验证。5.2 用cvpartition做五折交叉验证Matlab的cvpartition可以把样本划分为k份每轮轮流保留一份做测试。下面的代码输出每一折的RMSE以及均值与标准差rng(42); cv cvpartition(n, KFold, 5); rmse_fold zeros(cv.NumTestSets, 1); for k 1:cv.NumTestSets trIdx cv.training(k); teIdx cv.test(k); % 训练集设计矩阵 X_tr [ones(sum(trIdx), 1), X_raw(trIdx, :)]; beta_k X_tr \ y(trIdx); % 测试集预测 X_te [ones(sum(teIdx), 1), X_raw(teIdx, :)]; y_pred_k X_te * beta_k; rmse_fold(k) sqrt(mean((y(teIdx) - y_pred_k).^2)); end fprintf(5折CV RMSE: 均值 %.4f, 标准差 %.4f\n, ... mean(rmse_fold), std(rmse_fold));cv.training(k)和cv.test(k)返回逻辑索引sum(teIdx)统计测试样本数逻辑索引直接作用在X_raw行上。每一折训练和测试的数据都不重叠最终得到的均值就是模型泛化误差的更可靠估计。如果标准差接近均值数量级说明模型对数据划分很敏感要排查是否有异常样本或特征数量过少。如果后续要计算预测区间还需要检验残差正态性。可以用lillietest直接检查H lillietest(resid_test); % H0 表示不能拒绝残差服从正态分布的假设把交叉验证的RMSE均值和单次划分的测试RMSE放在一起对比二者接近说明模型性能可信单次划分远低于交叉验证均值说明这次划分偏乐观真正部署时要按均值预留误差余量。这一步是MLR回归项目交付前最值得做的验收动作。本文还有配套的精品资源点击获取