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

1. 从“调包”到“通透”：为什么你需要亲手实现多元线性回归

如果你正在学习数学建模，或者任何与数据分析、机器学习相关的课程，那么“多元线性回归”这个词对你来说一定不陌生。在MATLAB里，你很可能已经熟练地敲下fitlm或者regress这些函数，几秒钟后，模型系数、R²、p值等结果就整齐地呈现在眼前。这很方便，对吧？但问题也恰恰出在这里——这种“黑箱”式的操作，让你错过了理解模型最核心、最精妙的部分。

我见过太多同学，在建模比赛中，面对一个复杂的多元回归问题，只会调用内置函数。当模型结果不理想、出现多重共线性、或者需要对残差进行深入分析时，就完全束手无策了。他们知道“是什么”，但完全不明白“为什么”以及“怎么办”。这就像你学会了开车，但不知道引擎盖下发生了什么，一旦抛锚，就只能等待救援。

所以，今天我们不谈怎么用fitlm最快地跑出一个结果。我们要做的是，从最根本的矩阵运算出发，用MATLAB自编程，一步步“徒手”实现多元线性回归的整个流程。这个过程会让你彻底明白：

1. 那个神秘的回归系数β到底是怎么算出来的？
2. R²、调整R²这些评价指标背后的数学意义是什么？
3. 如何从零开始进行模型的显著性检验（F检验）和系数的显著性检验（t检验）？
4. 当遇到“设计矩阵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 × (p+1)的矩阵。注意，它多了一列，这是为了容纳截距项β₀。通常，我们会将这一列全部设为1，称为“全1列”或“截距列”。所以X看起来是这样的：

  X = [1, x₁₁, x₁₂, ..., x₁ₚ; 1, x₂₁, x₂₂, ..., x₂ₚ; ... ... ... ... ...; 1, xₙ₁, xₙ₂, ..., xₙₚ]

  第一列全是1，后面p列对应p个自变量的观测值。
• 系数向量 β：一个(p+1) × 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关于β求导，并令导数为零。经过推导（这里不展开微积分过程），我们可以得到著名的正规方程：(X'X) β = X'Y其中，X'是X的转置。如果X'X这个矩阵是可逆的（即X是列满秩的），那么我们就可以直接解出β的估计值：β_hat = (X'X)⁻¹ X'Y这就是整个多元线性回归最核心的公式。我们后续的MATLAB代码，将紧紧围绕这个公式展开。

注意：X'X的可逆性是一个关键前提。如果X的列之间存在严格的线性关系（例如，一个变量是另外两个变量的和），或者观测数n小于变量数p+1，那么X'X就是奇异矩阵，不可逆。这就是我们常说的“多重共线性”问题或“欠定”问题。在自编程时，我们必须处理这种异常情况。

2.3 模型评估：不止是R²

得到β_hat后，我们可以计算预测值Ŷ = X * β_hat和残差e = Y - Ŷ。如何评价这个模型的好坏呢？

1. 总平方和（SST）：SST = Σ(yᵢ - ȳ)² = (Y - ȳ)'(Y - ȳ)，其中ȳ是Y的均值。它代表了因变量的总波动。
2. 回归平方和（SSR）：SSR = Σ(ŷᵢ - ȳ)² = (Ŷ - ȳ)'(Ŷ - ȳ)。它代表了模型解释掉的那部分波动。
3. 误差平方和（SSE）：SSE = Σ(yᵢ - ŷᵢ)² = e'e。它代表了模型未能解释的波动。
4. 决定系数 R²：R² = SSR / SST = 1 - SSE/SST。它衡量了模型对数据波动的解释比例，介于0到1之间，越接近1越好。
5. 调整后 R²：当自变量增加时，R² 总会增加，即使这个变量无关紧要。调整R² 引入了惩罚项：Adj-R² = 1 - [(SSE/(n-p-1)) / (SST/(n-1))]。它更适用于比较不同自变量数量的模型。
6. 残差标准误（RSE）：RSE = 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) = σ² (X'X)⁻¹，其中σ²是误差项方差的估计，通常用MSE = SSE/(n-p-1)来估计。
  • 系数βⱼ的标准误（Standard Error）就是Cov(β_hat)矩阵第j个对角线元素的平方根：SE(βⱼ) = sqrt( MSE * ((X'X)⁻¹)[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 × (p+1)的设计矩阵，第一列全是1，用于估计截距项β₀。这是实现中非常容易忘记但至关重要的一步。

3.2 核心计算：求解正规方程与系数估计

接下来，我们根据公式β_hat = (X'X)⁻¹ X'Y来计算系数。

% 3. 核心：利用最小二乘法求解回归系数 % 计算 X'X 和 X'Y XtX = X_design' * X_design; XtY = X_design' * Y; % 检查 X'X 是否可逆（满秩） 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(i+1)); end fprintf('真实系数为: [%.4f, %.4f, %.4f, %.4f]\n', true_beta);

这里我引入了一个重要的实操检查：rank(XtX) < size(XtX, 1)。这个条件判断X'X是否满秩。如果不满秩，直接使用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, SSR+SSE); % 决定系数 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(i+1), se_beta(i+1), t_stats(i+1), p_values_t(i+1)); end

在这段代码中，fcdf和tcdf是MATLAB统计工具箱中的函数，用于计算F分布和t分布的累积概率。即使我们自编程核心算法，使用这些基础的分布函数也是合理且高效的。关键在于，我们知道了传入的参数（F_stat，df1，df2）是如何计算出来的。

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(:, j+1) = []; % 删除第j列自变量（保留截距项列） Y_this = X_design(:, j+1); % 第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。这就是共线性的信号。

处理方法：

1. 剔除变量：如果某个高VIF的变量在业务上不重要，可以直接剔除。
2. 主成分回归（PCR）：利用我们自编程的框架，可以轻松尝试。先对X_raw进行主成分分析（PCA），得到互不相关的主成分得分，然后用这些主成分作为新的自变量进行回归。这能彻底消除共线性，但缺点是主成分的解释性变差。
3. 岭回归（Ridge Regression）：这是处理共线性的经典方法。它在最小二乘法的损失函数中加入了一个对系数大小的惩罚项（L2正则化），公式变为：β_hat_ridge = (X'X + λI)⁻¹ X'Y，其中λ是惩罚系数，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); % 添加y=0参考线 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,4+i); 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: (p+1) 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(i+1), se_beta(i+1), t_stats(i+1), t_pvals(i+1)); end fprintf('==============================\n'); end

将这个函数保存为my_linregress.m文件。以后在任何项目中，你都可以像使用fitlm一样使用它，但你对它的内部逻辑了如指掌。你还可以在此基础上继续扩展，比如添加VIF计算、岭回归、逐步回归等功能，逐步构建起属于你自己的、功能强大的建模工具箱。

走到这里，你已经不再是那个只会调用fitlm的建模新手了。你亲手搭建了多元线性回归的每一个部件，理解了从数据输入到统计推断的完整链条。下次面对数据，你将有底气选择最合适的工具，并有能力在模型“出错”时，深入其内部进行诊断和修复。这才是数学建模和数据分析的真正能力所在——不是记住几个函数名，而是掌握驱动这些函数背后的思想与原理。