从调包到通透：手把手教你用MATLAB自编程实现多元线性回归

1. 从“调包”到“通透”为什么你需要亲手实现多元线性回归如果你正在学习数学建模或者任何与数据分析、机器学习相关的课程那么“多元线性回归”这个词对你来说一定不陌生。在MATLAB里你很可能已经熟练地敲下fitlm或者regress这些函数几秒钟后模型系数、R²、p值等结果就整齐地呈现在眼前。这很方便对吧但问题也恰恰出在这里——这种“黑箱”式的操作让你错过了理解模型最核心、最精妙的部分。我见过太多同学在建模比赛中面对一个复杂的多元回归问题只会调用内置函数。当模型结果不理想、出现多重共线性、或者需要对残差进行深入分析时就完全束手无策了。他们知道“是什么”但完全不明白“为什么”以及“怎么办”。这就像你学会了开车但不知道引擎盖下发生了什么一旦抛锚就只能等待救援。所以今天我们不谈怎么用fitlm最快地跑出一个结果。我们要做的是从最根本的矩阵运算出发用MATLAB自编程一步步“徒手”实现多元线性回归的整个流程。这个过程会让你彻底明白那个神秘的回归系数β到底是怎么算出来的R²、调整R²这些评价指标背后的数学意义是什么如何从零开始进行模型的显著性检验F检验和系数的显著性检验t检验当遇到“设计矩阵X不是满秩”这种常见错误时其根源是什么又该如何诊断和处理通过亲手实现你收获的将不仅仅是一个可运行的代码而是一种对模型“庖丁解牛”般的掌控感。下次再遇到回归问题你将是那个能洞察数据本质、能灵活调整方案、能解决棘手问题的“司机”而不是只会按按钮的“乘客”。我们这就开始。2. 核心原理拆解多元线性回归的“骨架”与“灵魂”在动手写代码之前我们必须把多元线性回归的数学骨架搭清楚。这是后续一切操作的基础理解了它你就能看透大多数回归类模型的本质。2.1 模型表述从公式到矩阵多元线性回归模型的标准形式是y β₀ β₁x₁ β₂x₂ ... βₚxₚ ε其中y是因变量x₁, x₂, ..., xₚ是p个自变量β₀是截距项β₁到βₚ是各自变量的系数ε是随机误差项通常假设其服从均值为0的正态分布。这个公式很直观但不利于计算。为了能利用矩阵运算的强大力量这也是MATLAB的核心优势我们需要将其改写为矩阵形式。假设我们有n组观测数据那么因变量向量 Y一个n × 1的列向量Y [y₁, y₂, ..., yₙ]’。设计矩阵 X这是一个n × (p1)的矩阵。注意它多了一列这是为了容纳截距项β₀。通常我们会将这一列全部设为1称为“全1列”或“截距列”。所以X看起来是这样的X [1, x₁₁, x₁₂, ..., x₁ₚ; 1, x₂₁, x₂₂, ..., x₂ₚ; ... ... ... ... ...; 1, xₙ₁, xₙ₂, ..., xₙₚ]第一列全是1后面p列对应p个自变量的观测值。系数向量 β一个(p1) × 1的列向量β [β₀, β₁, β₂, ..., βₚ]’。误差向量 ε一个n × 1的列向量ε [ε₁, ε₂, ..., εₙ]’。于是整个模型可以优雅地写为Y Xβ ε。这就是多元线性回归的矩阵形式它把n个方程压缩成了一个简洁的矩阵方程。2.2 参数估计最小二乘法的几何与代数视角我们的目标是找到一组系数β使得模型预测值Ŷ Xβ与真实值Y之间的差距最小。这个差距用误差平方和Sum of Squared Errors, SSE来衡量SSE Σ(yᵢ - ŷᵢ)² (Y - Xβ)(Y - Xβ)。最小二乘法就是寻找使SSE达到最小的那个β。从几何上看Ŷ Xβ是Y在由X的列向量所张成的线性空间上的投影。最小二乘解就是找到Y在这个空间上的垂直投影点此时残差向量e Y - Ŷ与X的列空间垂直。从代数上求解需要对SSE关于β求导并令导数为零。经过推导这里不展开微积分过程我们可以得到著名的正规方程(XX) β XY其中X是X的转置。如果XX这个矩阵是可逆的即X是列满秩的那么我们就可以直接解出β的估计值β_hat (XX)⁻¹ XY这就是整个多元线性回归最核心的公式。我们后续的MATLAB代码将紧紧围绕这个公式展开。注意XX的可逆性是一个关键前提。如果X的列之间存在严格的线性关系例如一个变量是另外两个变量的和或者观测数n小于变量数p1那么XX就是奇异矩阵不可逆。这就是我们常说的“多重共线性”问题或“欠定”问题。在自编程时我们必须处理这种异常情况。2.3 模型评估不止是R²得到β_hat后我们可以计算预测值Ŷ X * β_hat和残差e Y - Ŷ。如何评价这个模型的好坏呢总平方和SSTSST Σ(yᵢ - ȳ)² (Y - ȳ)(Y - ȳ)其中ȳ是Y的均值。它代表了因变量的总波动。回归平方和SSRSSR Σ(ŷᵢ - ȳ)² (Ŷ - ȳ)(Ŷ - ȳ)。它代表了模型解释掉的那部分波动。误差平方和SSESSE Σ(yᵢ - ŷᵢ)² ee。它代表了模型未能解释的波动。决定系数 R²R² SSR / SST 1 - SSE/SST。它衡量了模型对数据波动的解释比例介于0到1之间越接近1越好。调整后 R²当自变量增加时R² 总会增加即使这个变量无关紧要。调整R² 引入了惩罚项Adj-R² 1 - [(SSE/(n-p-1)) / (SST/(n-1))]。它更适用于比较不同自变量数量的模型。残差标准误RSERSE sqrt(SSE / (n-p-1))。它可以理解为模型预测的平均误差大小其单位与Y相同非常直观。2.4 统计推断模型与系数是否真的有用算出了 R²模型就一定显著吗不一定。我们需要进行统计检验。模型的显著性检验F检验原假设 H₀所有自变量的系数都为0即β₁ β₂ ... βₚ 0模型没有意义。构造 F 统计量F (SSR/p) / (SSE/(n-p-1))。在原假设下F服从自由度为(p, n-p-1)的 F 分布。计算出的F值越大对应的 p-value 越小我们就越有理由拒绝原假设认为模型整体是显著的。系数的显著性检验t检验对于某一个系数βⱼ我们关心它是否显著不为0。首先需要估计系数估计值β_hat的方差-协方差矩阵Cov(β_hat) σ² (XX)⁻¹其中σ²是误差项方差的估计通常用MSE SSE/(n-p-1)来估计。系数βⱼ的标准误Standard Error就是Cov(β_hat)矩阵第j个对角线元素的平方根SE(βⱼ) sqrt( MSE * ((XX)⁻¹)[j,j] )。构造 t 统计量tⱼ β_hatⱼ / SE(βⱼ)。在原假设 H₀βⱼ 0下tⱼ服从自由度为n-p-1的 t 分布。据此可以计算每个系数的 p-value。看到这里你可能觉得头大。但请放心接下来的MATLAB代码会把这些抽象的公式全部变成具体的计算步骤。你会发现一旦理解了原理编程实现就是水到渠成的事情。3. MATLAB自编程实现一步步构建你的回归工具箱现在我们进入实战环节。我将带领你不依赖任何统计工具箱函数除了基础的矩阵运算完整实现多元线性回归。3.1 数据准备与设计矩阵构建首先我们需要一些数据。这里我们使用一个经典的例子波士顿房价数据集虽然现在不鼓励用了但用于教学非常清晰。我们会模拟一个类似的数据结构。% 清空环境 clear; clc; % 1. 模拟生成数据 n 100; % 100个样本 p 3; % 3个自变量 % 生成自变量X加入一些相关性以模拟真实情况 rng(2023); % 设定随机种子确保结果可复现 X_raw randn(n, p); % 让X2和X3有一定相关性 X_raw(:,3) 0.7 * X_raw(:,2) 0.3 * randn(n,1); % 定义真实系数 true_beta [2.5; -1.2; 0.8; 3.0]; % 第一个是截距项beta0 % 生成因变量Y: Y 2.5 -1.2*X1 0.8*X2 3.0*X3 noise noise 0.5 * randn(n,1); Y true_beta(1) X_raw * true_beta(2:end) noise; % 2. 构建设计矩阵X_design % 添加全1列以对应截距项beta0 X_design [ones(n,1), X_raw]; disp(设计矩阵X_design的前5行); disp(X_design(1:5, :)); disp(因变量Y的前5个值); disp(Y(1:5));这段代码的关键在于X_design [ones(n,1), X_raw]。这行代码构建了理论部分提到的n × (p1)的设计矩阵第一列全是1用于估计截距项β₀。这是实现中非常容易忘记但至关重要的一步。3.2 核心计算求解正规方程与系数估计接下来我们根据公式β_hat (XX)⁻¹ XY来计算系数。% 3. 核心利用最小二乘法求解回归系数 % 计算 XX 和 XY XtX X_design * X_design; XtY X_design * Y; % 检查 XX 是否可逆满秩 if rank(XtX) size(XtX, 1) warning(设计矩阵X不满秩存在多重共线性直接求逆可能不稳定。建议使用岭回归或剔除变量。); % 一种稳健的解法是使用伪逆 pinv beta_hat pinv(XtX) * XtY; else % 如果满秩直接求逆 beta_hat inv(XtX) * XtY; end fprintf(\n--- 回归系数估计结果 ---\n); fprintf(截距项 (beta0): %.4f\n, beta_hat(1)); for i 1:p fprintf(系数 beta%d (对应X%d): %.4f\n, i, i, beta_hat(i1)); end fprintf(真实系数为: [%.4f, %.4f, %.4f, %.4f]\n, true_beta);这里我引入了一个重要的实操检查rank(XtX) size(XtX, 1)。这个条件判断XX是否满秩。如果不满秩直接使用inv函数求逆会得到NaN或极不稳定的结果。在这种情况下使用pinv伪逆基于奇异值分解是一种更稳健的数值解法但它给出的是一种最小二乘解需要结合业务理解来解读。在数学建模中遇到这种情况你更应该去检查数据是否存在多重共线性而不是简单地用pinv绕过。3.3 模型预测、残差与拟合优度计算有了系数我们就可以进行预测并计算关键的评估指标。% 4. 模型预测与残差计算 Y_hat X_design * beta_hat; % 预测值 residuals Y - Y_hat; % 残差 % 5. 计算拟合优度指标 Y_mean mean(Y); % 总平方和 SST SST sum((Y - Y_mean).^2); % 回归平方和 SSR SSR sum((Y_hat - Y_mean).^2); % 误差平方和 SSE SSE sum(residuals.^2); % 验证 SST SSR SSE (在数值计算允许的误差内) fprintf(\nSST %.4f, SSR SSE %.4f\n, SST, SSRSSE); % 决定系数 R-squared R_squared SSR / SST; % 调整后的 R-squared n length(Y); k p; % 自变量个数不含截距 adj_R_squared 1 - (SSE/(n-k-1)) / (SST/(n-1)); % 残差标准误 RSE / 均方根误差 RMSE RMSE sqrt(SSE / (n-k-1)); fprintf(\n--- 模型拟合优度 ---\n); fprintf(R-squared: %.4f\n, R_squared); fprintf(Adjusted R-squared: %.4f\n, adj_R_squared); fprintf(均方根误差 RMSE: %.4f\n, RMSE); fprintf(残差和 (应为接近0): %.6f\n, sum(residuals));注意理论上最小二乘估计保证残差和为零。但在实际数值计算中由于浮点数精度问题sum(residuals)可能是一个极小的数如1e-14而不是绝对的0。这是一个很好的检查点如果这个值很大说明你的计算过程可能有误。3.4 统计推断F检验与t检验的实现这是自编程中最体现价值的部分让我们看清统计检验的每一个细节。% 6. 模型的显著性检验 (F检验) % 回归均方 MSR MSR SSR / k; % 残差均方 MSE MSE SSE / (n - k - 1); % F统计量 F_stat MSR / MSE; % F检验的p值 p_value_F 1 - fcdf(F_stat, k, n-k-1); % fcdf是F分布的累积分布函数 fprintf(\n--- 模型显著性检验 (F检验) ---\n); fprintf(F统计量: %.4f\n, F_stat); fprintf(自由度 (回归, 残差): (%d, %d)\n, k, n-k-1); fprintf(F检验的p值: %.6f\n, p_value_F); if p_value_F 0.05 fprintf(结论: 在0.05显著性水平下拒绝原假设模型整体显著。\n); else fprintf(结论: 在0.05显著性水平下无法拒绝原假设模型整体不显著。\n); end % 7. 回归系数的显著性检验 (t检验) % 计算系数估计的方差-协方差矩阵 cov_beta MSE * inv(XtX); % 这里假设XtX可逆否则用pinv(XtX) % 提取系数的标准误 se_beta sqrt(diag(cov_beta)); % 计算t统计量 t_stats beta_hat ./ se_beta; % 计算每个系数对应的p值 (双尾检验) p_values_t 2 * (1 - tcdf(abs(t_stats), n-k-1)); % tcdf是t分布的累积分布函数 fprintf(\n--- 回归系数显著性检验 (t检验) ---\n); fprintf(%10s %10s %10s %10s %12s\n, 系数, 估计值, 标准误, t统计量, p值); fprintf(%10s %10.4f %10.4f %10.4f %12.6f\n, beta0, beta_hat(1), se_beta(1), t_stats(1), p_values_t(1)); for i 1:p fprintf(beta%d(X%d) %10.4f %10.4f %10.4f %12.6f\n, i, i, beta_hat(i1), se_beta(i1), t_stats(i1), p_values_t(i1)); end在这段代码中fcdf和tcdf是MATLAB统计工具箱中的函数用于计算F分布和t分布的累积概率。即使我们自编程核心算法使用这些基础的分布函数也是合理且高效的。关键在于我们知道了传入的参数F_statdf1df2是如何计算出来的。3.5 与MATLAB内置函数对比验证为了确保我们的自编程结果是正确的最好的方法就是与MATLAB内置的权威函数进行对比。% 8. 使用MATLAB内置函数进行验证 % 使用 fitlm 函数 (需要Statistics and Machine Learning Toolbox) if exist(fitlm, file) 2 % 将X_raw作为表格变量传入fitlm会自动添加截距项 tbl array2table(X_raw, VariableNames, {X1, X2, X3}); tbl.Y Y; mdl fitlm(tbl, Y ~ X1 X2 X3); fprintf(\n 与MATLAB fitlm函数对比 \n); fprintf(\n1. 系数对比:\n); disp(自编程结果:); disp(beta_hat); disp(fitlm结果 (Coefficients.Estimate):); disp(mdl.Coefficients.Estimate); fprintf(\n2. 拟合优度对比:\n); fprintf(自编程 - R²: %.6f, Adj-R²: %.6f, RMSE: %.6f\n, R_squared, adj_R_squared, RMSE); fprintf(fitlm - R²: %.6f, Adj-R²: %.6f, RMSE: %.6f\n, mdl.Rsquared.Ordinary, mdl.Rsquared.Adjusted, mdl.RMSE); fprintf(\n3. 整体F检验对比:\n); fprintf(自编程 - F: %.6f, p-value: %.6f\n, F_stat, p_value_F); fprintf(fitlm - F: %.6f, p-value: %.6f\n, mdl.ModelFitVsNullModel.Fstat, mdl.ModelFitVsNullModel.Pvalue); fprintf(\n4. 系数t检验对比 (以beta1为例):\n); fprintf(自编程 - t: %.6f, p-value: %.6f\n, t_stats(2), p_values_t(2)); fprintf(fitlm - t: %.6f, p-value: %.6f\n, mdl.Coefficients.tStat(2), mdl.Coefficients.pValue(2)); else fprintf(\n未检测到Statistics and Machine Learning Toolbox跳过fitlm对比。\n); % 可以使用 regress 函数对比 [b, bint, r, rint, stats] regress(Y, X_design); fprintf(使用regress函数对比系数:\n); disp(自编程结果:); disp(beta_hat); disp(regress结果:); disp(b); fprintf(regress返回的stats向量 [R², F, p-value, 误差方差估计]:\n); disp(stats); end运行这段对比代码如果你的自编程结果与fitlm或regress的输出在数值上高度一致可能在小数点后第10位有细微差异源于浮点数计算那么恭喜你你的自编程实现是完全正确的这个对比过程不仅能验证代码更能给你巨大的信心。4. 超越基础自编程如何帮你解决实际问题如果你只是调用fitlm那么当结果不如预期时你的调试手段非常有限。而自编程赋予了你“透视”整个建模过程的能力让你能主动诊断和解决复杂问题。4.1 诊断与处理多重共线性多重共线性是多元回归中的常见病。它不会影响模型的预测能力但会使系数的估计值方差变大变得非常不稳定难以解释。fitlm会给出警告但自编程能让你更深入地理解它。诊断方法方差膨胀因子方差膨胀因子VIF是衡量共线性严重程度的常用指标。对于第j个自变量其 VIF 等于以该自变量为因变量对其他所有自变量进行回归所得到的 R² 的函数VIF_j 1 / (1 - R²_j)。VIF 大于 5 或 10 通常被认为存在较严重的共线性。% 计算方差膨胀因子 (VIF) fprintf(\n--- 多重共线性诊断方差膨胀因子(VIF) ---\n); vifs zeros(p, 1); for j 1:p % 将第j个自变量作为因变量其余自变量及截距项作为自变量 X_other X_design; X_other(:, j1) []; % 删除第j列自变量保留截距项列 Y_this X_design(:, j1); % 第j个自变量 % 使用我们自编的回归函数这里简单调用核心计算部分 XtX_other X_other * X_other; XtY_other X_other * Y_this; if rank(XtX_other) size(XtX_other,1) b_other inv(XtX_other) * XtY_other; else b_other pinv(XtX_other) * XtY_other; end Y_hat_other X_other * b_other; SST_other sum((Y_this - mean(Y_this)).^2); SSE_other sum((Y_this - Y_hat_other).^2); R2_j 1 - SSE_other / SST_other; vifs(j) 1 / (1 - R2_j); fprintf(自变量 X%d 的 VIF: %.4f\n, j, vifs(j)); end在我们模拟的数据中X2和X3被设计为有相关性你可能会看到X2或X3的 VIF 值显著高于X1。这就是共线性的信号。处理方法剔除变量如果某个高VIF的变量在业务上不重要可以直接剔除。主成分回归PCR利用我们自编程的框架可以轻松尝试。先对X_raw进行主成分分析PCA得到互不相关的主成分得分然后用这些主成分作为新的自变量进行回归。这能彻底消除共线性但缺点是主成分的解释性变差。岭回归Ridge Regression这是处理共线性的经典方法。它在最小二乘法的损失函数中加入了一个对系数大小的惩罚项L2正则化公式变为β_hat_ridge (XX λI)⁻¹ XY其中λ是惩罚系数I是单位阵。通过自编程实现岭回归并观察系数路径随λ变化的情况是理解正则化的绝佳方式。4.2 残差分析验证模型假设线性回归模型的有效性建立在几个关键假设上误差项独立、同方差、正态分布。这些假设是否成立需要通过残差分析来检验。自编程让你能自由地绘制和分析残差图。% 残差分析绘图 figure(Position, [100, 100, 1200, 800]); % 1. 残差与拟合值图 (检查同方差性) subplot(2,3,1); plot(Y_hat, residuals, o); hold on; plot([min(Y_hat), max(Y_hat)], [0,0], r--, LineWidth, 1.5); % 添加y0参考线 xlabel(拟合值 \^Y); ylabel(残差 e); title(残差 vs. 拟合值); grid on; % 理想情况残差随机均匀分布在0线两侧无明显趋势或漏斗形状。 % 2. 残差的正态概率图 (QQ图) subplot(2,3,2); normplot(residuals); title(残差正态概率图 (QQ图)); % 理想情况点大致分布在一条直线上。 % 3. 残差序列图 (检查独立性) subplot(2,3,3); plot(1:n, residuals, o-); hold on; plot([1, n], [0,0], r--, LineWidth, 1.5); xlabel(观测序号); ylabel(残差 e); title(残差序列图); grid on; % 理想情况残差随机波动无明显的周期性或趋势性。 % 4. 残差直方图 subplot(2,3,4); histogram(residuals, 15, Normalization, pdf); hold on; % 叠加正态分布曲线 mu mean(residuals); sigma std(residuals); x linspace(min(residuals), max(residuals), 100); y normpdf(x, mu, sigma); plot(x, y, r-, LineWidth, 2); xlabel(残差); ylabel(概率密度); title(残差分布直方图); legend(残差分布, 正态分布拟合); grid on; % 5. 残差与各自变量的关系图 (检查线性假设) for i 1:min(p, 3) % 只画前三个自变量 subplot(2,3,4i); plot(X_raw(:,i), residuals, o); hold on; plot([min(X_raw(:,i)), max(X_raw(:,i))], [0,0], r--, LineWidth, 1.5); xlabel(sprintf(自变量 X%d, i)); ylabel(残差 e); title(sprintf(残差 vs. X%d, i)); grid on; end sgtitle(多元线性回归残差分析图);通过观察这些图你可以判断残差-拟合值图如果残差随拟合值增大而扩散或收敛漏斗形则违反同方差假设可能需要进行变量变换如对数变换或使用加权最小二乘法。QQ图如果点严重偏离直线特别是两端则误差可能非正态。对于大样本量中心极限定理通常能保证推断的稳健性但对于小样本需谨慎。残差序列图如果残差呈现明显的趋势或周期说明误差项可能自相关常见于时间序列数据。这会影响标准误的估计可能需要使用时间序列模型。4.3 模型优化与特征工程尝试自编程的最大优势是灵活性。你可以轻松地尝试各种特征工程并立即评估其效果。例如尝试加入交互项或多项式项% 假设我们认为 X1 和 X2 可能存在交互效应 X_raw_with_interaction [X_raw, X_raw(:,1) .* X_raw(:,2)]; % 添加交互项 X1*X2 % 或者尝试加入 X1 的平方项 % X_raw_with_poly [X_raw, X_raw(:,1).^2]; % 然后用我们自编的函数重新拟合模型 % 只需修改数据准备部分后续计算代码完全复用 X_design_new [ones(n,1), X_raw_with_interaction]; % ... (重复3.2到3.4的计算步骤) % 计算新的R², Adj-R², F, p-value等 % 对比新模型和旧模型看Adj-R²是否提升新加入的项是否显著(t检验)。通过这种快速的迭代和对比你可以基于数据证据而非直觉来决定最终的特征组合。这是建模过程中最具创造性的部分而自编程为你提供了实现这种创造性的完整工具链。5. 封装与复用构建你自己的回归函数库将上述代码模块化封装成函数是工程实践的必然步骤。这不仅能让你在未来的项目中快速调用更是对知识体系的巩固。function [beta, stats] my_linregress(X, Y) % MY_LINREGRESS 自编多元线性回归函数 % 输入 % X: n x p 矩阵n个样本p个特征不含截距项 % Y: n x 1 向量因变量 % 输出 % beta: (p1) x 1 向量回归系数 [beta0; beta1; ...; betap] % stats: 结构体包含模型统计量 % .Y_hat: 预测值 % .residuals: 残差 % .R2: 决定系数 % .R2_adj: 调整决定系数 % .RMSE: 均方根误差 % .F_stat: F检验统计量 % .F_pval: F检验p值 % .t_stats: t检验统计量向量 % .t_pvals: t检验p值向量 % .cov_beta: 系数协方差矩阵 % .se_beta: 系数标准误向量 [n, p] size(X); % 1. 构建设计矩阵 X_design [ones(n,1), X]; % 2. 求解系数 (使用伪逆以增强稳定性) XtX X_design * X_design; XtY X_design * Y; beta pinv(XtX) * XtY; % 使用pinv替代inv % 3. 预测与残差 Y_hat X_design * beta; residuals Y - Y_hat; % 4. 拟合优度 Y_mean mean(Y); SST sum((Y - Y_mean).^2); SSR sum((Y_hat - Y_mean).^2); SSE sum(residuals.^2); R2 SSR / SST; R2_adj 1 - (SSE/(n-p-1)) / (SST/(n-1)); RMSE sqrt(SSE / (n-p-1)); % 5. F检验 MSR SSR / p; MSE_val SSE / (n-p-1); F_stat MSR / MSE_val; F_pval 1 - fcdf(F_stat, p, n-p-1); % 6. t检验 cov_beta MSE_val * pinv(XtX); % 使用与系数估计一致的伪逆 se_beta sqrt(diag(cov_beta)); t_stats beta ./ se_beta; t_pvals 2 * (1 - tcdf(abs(t_stats), n-p-1)); % 7. 打包输出 stats struct(); stats.Y_hat Y_hat; stats.residuals residuals; stats.R2 R2; stats.R2_adj R2_adj; stats.RMSE RMSE; stats.F_stat F_stat; stats.F_pval F_pval; stats.t_stats t_stats; stats.t_pvals t_pvals; stats.cov_beta cov_beta; stats.se_beta se_beta; % 8. 简单的结果打印 fprintf(\n 自编线性回归结果 \n); fprintf(样本数 n %d, 自变量数 p %d\n, n, p); fprintf(R² %.4f, Adj-R² %.4f, RMSE %.4f\n, R2, R2_adj, RMSE); fprintf(模型F检验: F(%d, %d) %.4f, p %.6f\n, p, n-p-1, F_stat, F_pval); fprintf(\n系数估计与检验:\n); fprintf(%8s %12s %12s %12s %12s\n, 变量, 系数估计, 标准误, t值, p值); fprintf(%8s %12.4f %12.4f %12.4f %12.6f\n, 截距, beta(1), se_beta(1), t_stats(1), t_pvals(1)); for i 1:p fprintf(X%-7d %12.4f %12.4f %12.4f %12.6f\n, i, beta(i1), se_beta(i1), t_stats(i1), t_pvals(i1)); end fprintf(\n); end将这个函数保存为my_linregress.m文件。以后在任何项目中你都可以像使用fitlm一样使用它但你对它的内部逻辑了如指掌。你还可以在此基础上继续扩展比如添加VIF计算、岭回归、逐步回归等功能逐步构建起属于你自己的、功能强大的建模工具箱。走到这里你已经不再是那个只会调用fitlm的建模新手了。你亲手搭建了多元线性回归的每一个部件理解了从数据输入到统计推断的完整链条。下次面对数据你将有底气选择最合适的工具并有能力在模型“出错”时深入其内部进行诊断和修复。这才是数学建模和数据分析的真正能力所在——不是记住几个函数名而是掌握驱动这些函数背后的思想与原理。