1. 从一次失败的曲线拟合说起前阵子我帮一个做材料分析的朋友处理一组实验数据。他测了不同温度下某种材料的电阻率想找出两者之间的数学关系。数据点散落在图上看起来像是一条直线但又没那么“直”。他随手用Excel画了条趋势线拿了个R²值就准备写进报告里。我一看这线画得有点“任性”——它似乎更照顾了几个数值偏大的点而对中间密集的小数据点“视而不见”。结果就是用这条线去预测新温度下的电阻率偏差可能会很大。这其实就是最朴素、也最容易出问题的“目测拟合”。我们真正需要的是一条能“公平”对待所有数据点的线让每个数据点到这条线的垂直距离我们称之为“残差”的平方和最小。这个方法就是最小二乘法。它不是什么新潮的算法而是两百多年前高斯和勒让德玩转天文学数据时就奠定的基石至今仍是数据分析、机器学习、工程优化的核心工具。为什么是“平方和”最小而不是绝对值和或者四次方和这背后有深刻的概率统计原理简单说在误差服从正态分布的假设下它给出的是最可能接近真实情况的估计但更直观的原因是平方项求导方便能导出一个干净利落的解析解。而Matlab几乎是实现这个想法的“标准答案”。它内置的矩阵运算和“反斜杠”运算符让最小二乘求解变得像做除法一样简单。但“简单”背后藏着很多新手甚至老手都会忽略的细节你的数据适合用直线拟合吗拟合出的系数物理意义是什么那个R²值到底可信度有多高Matlab给出的那一大堆统计信息又该怎么解读这篇文章我就以这个电阻率-温度数据的案例为引子带你从零开始在Matlab里亲手实现最小二乘法。我们不只满足于调通一行代码更要拆开“黑箱”看看Matlab到底帮我们算了什么并在这个过程中避开那些教科书里不提、但实践中一定会踩的坑。2. 最小二乘法的数学内核不止是“背公式”在打开Matlab之前我们必须先搞清楚我们要计算的对象是什么。很多人对最小二乘法的印象停留在y kx b这个公式上但这只是特例。它的核心思想是线性模型注意是“参数的线性”而非“变量的线性”。2.1 问题的一般化表述假设我们有n组观测数据(x_i, y_i), i1,2,...,n。我们认为y和x之间存在一种关系可以用一组基函数φ_j(x)的线性组合来近似描述y ≈ β_1*φ_1(x) β_2*φ_2(x) ... β_p*φ_p(x)这里的β_1, β_2, ..., β_p就是我们要求解的未知系数。φ_j(x)可以是任意函数。例如线性拟合φ_1(x)1,φ_2(x)x。模型为y β_1 β_2*x。多项式拟合φ_1(x)1,φ_2(x)x,φ_3(x)x², ...。模型为y β_1 β_2*x β_3*x² ...。非线性关系的线性化比如指数衰减y a*exp(b*x)两边取对数得ln(y) ln(a) b*x令Yln(y),β_1ln(a),β_2b就化为了线性模型Y β_1 β_2*x。将n组数据代入我们会得到一个线性方程组用矩阵表示最为清晰| y1 | | φ_1(x1) φ_2(x1) ... φ_p(x1) | | β_1 | | ε_1 | | y2 | | φ_1(x2) φ_2(x2) ... φ_p(x2) | * | β_2 | | ε_2 | | .. | | ... | | ... | | ... | | yn | | φ_1(xn) φ_2(xn) ... φ_p(xn) | | β_p | | ε_n |简写为Y Xβ ε其中Y是n×1的观测值向量。X是n×p的设计矩阵每一列对应一个基函数在所有x_i上的取值。它的构造是建模的关键一步。β是p×1的待求系数向量。ε是n×1的误差向量代表模型无法解释的部分。2.2 最小二乘解与法方程“最小二乘”的目标就是找到一组系数β使得所有误差的平方和S Σ ε_i² εᵀε达到最小。将ε Y - Xβ代入得到目标函数S(β) (Y - Xβ)ᵀ (Y - Xβ)这是一个关于β的二次函数。为了求其最小值我们令其对β的梯度为零∇S(β) 0。经过推导这是理解最小二乘法的关键一步我们得到著名的正规方程(XᵀX) β XᵀY这个方程被称为“法方程”。只要设计矩阵X是列满秩的即各列线性无关XᵀX就是一个p×p的可逆方阵。那么最小二乘解β_hat就可以直接写出β_hat (XᵀX)⁻¹ XᵀY这就是最小二乘法的“心脏”。所有软件包括Matlab最终都是在以某种方式求解这个方程。理解了这个你就知道了Matlab里那个神奇的“\”运算符反斜杠在背后做了什么。注意这里隐藏着一个巨大的实践陷阱。直接计算(XᵀX)⁻¹ XᵀY在数学上是正确的但在数值计算上可能是不稳定的。如果X的列之间存在近似线性相关称为“多重共线性”或者数据量纲差异巨大XᵀX的条件数会变得非常大求逆会放大舍入误差导致结果极不可靠。Matlab的“\”运算符以及pinv函数使用了更稳定的数值算法如QR分解、SVD分解来绕过直接求逆这是专业工具的价值所在。我们自己编程实现时如果数据规模不大、性质良好可以用公式法但对于严肃的数据分析强烈建议使用稳定的数值方法。3. Matlab实战从“傻瓜式”操作到自定义实现现在让我们回到Matlab用三种由浅入深的方式来解决开头的电阻率拟合问题。假设我们有以下数据 温度T(℃): [20, 40, 60, 80, 100, 120] 电阻率ρ(Ω·m): [1.02, 1.15, 1.33, 1.51, 1.75, 2.01]3.1 方法一使用内置拟合工具最快上手对于快速可视化探索Matlab的图形化工具是无敌的。T [20, 40, 60, 80, 100, 120]; rho [1.02, 1.15, 1.33, 1.51, 1.75, 2.01]; figure; plot(T, rho, bo, MarkerFaceColor, b, MarkerSize, 8); xlabel(温度 T (℃)); ylabel(电阻率 \rho (\Omega\cdot m)); grid on;画好散点图后在图形窗口的菜单栏点击“工具” - “基本拟合”。在弹出的窗口中勾选“线性”拟合并勾选“显示方程”和“绘制残差图”。你会立刻得到拟合直线ρ 0.0099*T 0.828以及残差分布图。优点直观、快速适合数据初探。缺点过程不透明无法自动化难以进行复杂的后续统计分析如计算置信区间。3.2 方法二使用polyfit函数最常用对于多项式拟合polyfit是专用函数语法简洁。% 进行一次线性拟合 p polyfit(T, rho, 1); % 第三个参数‘1’代表一次多项式直线 % p 是一个向量p(1)是斜率kp(2)是截距b k p(1); b p(2); fprintf(拟合直线: rho %.4f * T %.4f\n, k, b); % 用拟合出的系数计算预测值 rho_fit polyval(p, T); % 计算R² SS_res sum((rho - rho_fit).^2); % 残差平方和 SS_tot sum((rho - mean(rho)).^2); % 总平方和 R2 1 - SS_res / SS_tot; fprintf(R-squared (R²) %.4f\n, R2); % 绘图 figure; plot(T, rho, bo, MarkerFaceColor, b, MarkerSize, 8); hold on; plot(T, rho_fit, r-, LineWidth, 2); legend(原始数据, 最小二乘拟合, Location, northwest); xlabel(温度 T (℃)); ylabel(电阻率 \rho (\Omega\cdot m)); title(sprintf(线性拟合: \\rho %.4fT %.4f (R^2%.4f), k, b, R2)); grid on;polyfit内部也是通过构造范德蒙德矩阵一种特殊的设计矩阵并求解最小二乘问题来实现的。对于更高次的多项式只需改变第三个参数即可。实操心得polyfit虽然方便但在拟合高次多项式如5次以上时极易出现数值不稳定问题因为范德蒙德矩阵是著名的病态矩阵。此时考虑使用正交多项式如polyfit的‘centering and scaling’选项即polyfit(T, rho, n, Scale, on)或直接使用更稳定的\运算符配合精心设计的设计矩阵。3.3 方法三使用“反斜杠”运算符\最通用、最核心这是Matlab解决线性最小二乘问题的“瑞士军刀”。它直接对应我们第二节中的矩阵方程Xβ Y。% 步骤1构造设计矩阵 X % 对于线性模型 y β1 β2*x设计矩阵第一列是全1对应截距第二列是x值 X [ones(length(T), 1), T]; % 注意将行向量转置为列向量 Y rho; % 步骤2使用反斜杠求解 β X \ Y beta X \ Y; % 这就是最小二乘解 b beta(1); % 截距 k beta(2); % 斜率 fprintf(使用 \\ 运算符拟合: rho %.4f * T %.4f\n, k, b); % 验证与polyfit结果一致应完全一致 rho_fit_backslash X * beta;这行beta X \ Y是精髓。当方程组Xβ Y超定方程数多于未知数时\自动计算最小二乘解。它内部会根据矩阵X的具体情况智能选择最高效、最稳定的数值算法如QR分解。3.4 方法四手动实现矩阵运算理解原理为了彻底搞懂我们可以手动实现第二节中的正规方程解法。再次警告此方法数值稳定性较差仅用于教学理解。% 手动计算最小二乘解beta (X*X)^(-1) * (X*Y) X [ones(length(T), 1), T]; Y rho; % 计算 XX 和 XY XTX X * X; XTY X * Y; % 求逆并求解这是不稳定的根源 beta_manual inv(XTX) * XTY; % 不推荐在实际中使用inv直接求逆 fprintf(手动矩阵求逆拟合: rho %.4f * T %.4f\n, beta_manual(2), beta_manual(1)); % 更稳定的手动实现使用QR分解 [Q, R] qr(X, 0); % 经济型QR分解R是上三角阵 beta_qr R \ (Q * Y); % 解上三角方程组数值上稳定得多 fprintf(手动QR分解拟合: rho %.4f * T %.4f\n, beta_qr(2), beta_qr(1));通过对比你会发现beta_manual、beta_qr和之前beta的结果在理论上应该一致。但在数据条件恶劣时beta_manual可能会产生显著误差。4. 结果评估与深入分析你的拟合“好”吗得到拟合系数只是第一步。一个负责任的建模者必须评估这个模型的可靠性。Matlab的统计工具箱提供了强大的工具但我们也可以手动计算关键指标。4.1 核心统计量计算除了R²我们还需要关注残差residuals Y - X*beta。残差应该随机分布在0附近没有明显的模式如趋势或周期性。绘制残差图是检验模型假设如误差同方差、独立性的重要手段。均方根误差RMSE sqrt(mean(residuals.^2))。它反映了模型预测的平均误差大小与因变量Y有相同量纲更直观。系数标准误衡量系数估计的精度。这需要计算误差方差估计和(XᵀX)⁻¹的对角线元素。% 接续方法三的代码 n length(Y); % 样本数 p size(X, 2); % 参数个数包括截距 % 计算残差和误差方差估计 residuals Y - X * beta; sigma2_hat (residuals * residuals) / (n - p); % 无偏估计 % 计算系数协方差矩阵 (XX)^(-1) * sigma2_hat % 使用更稳定的方式计算 (XX)的逆cov_beta sigma2_hat * inv(X*X); [Q,R] qr(X, 0); Rinv inv(R); % R是方阵且上三角求逆稳定且高效 cov_beta sigma2_hat * (Rinv * Rinv); % 系数标准误 协方差矩阵对角线的平方根 se_beta sqrt(diag(cov_beta)); fprintf(截距 b %.4f ± %.4f\n, beta(1), se_beta(1)); fprintf(斜率 k %.4f ± %.4f\n, beta(2), se_beta(2)); % t统计量和p值检验系数是否显著不为0 t_stat beta ./ se_beta; df n - p; % 自由度 p_value 2 * (1 - tcdf(abs(t_stat), df)); % 双尾检验 fprintf(斜率k的t统计量 %.4f, p值 %.4e\n, t_stat(2), p_value(2));如果斜率的p值远小于0.05例如 0.01我们通常认为温度对电阻率有显著影响。4.2 使用fitlm函数进行专业回归分析对于需要全面回归诊断的场景统计工具箱中的fitlm函数是终极武器。% 将数据转换为表fitlm推荐格式 tbl table(T, rho, VariableNames, {Temperature, Resistivity}); % 拟合线性模型 mdl fitlm(tbl, Resistivity ~ Temperature); % 公式表示Resistivity 由 Temperature 解释 % 显示完整的回归结果摘要 disp(mdl); % 获取关键信息 coefficients mdl.Coefficients; % 系数表包含Estimate, SE, tStat, pValue R2 mdl.Rsquared.Ordinary; RMSE mdl.RMSE; fprintf(\n模型摘要:\n); fprintf(R² %.4f\n, R2); fprintf(调整后R² %.4f\n, mdl.Rsquared.Adjusted); fprintf(RMSE %.4f\n, RMSE); % 绘制诊断图 figure; plotResiduals(mdl, fitted); % 残差 vs. 拟合值图检查同方差性 figure; plotDiagnostics(mdl, cookd); % Cook距离识别强影响点fitlm的输出包含了专业统计软件级别的信息如调整R²、F统计量、DW检验自相关等并提供了丰富的绘图函数进行模型诊断。踩坑实录我曾用一组存在异方差误差方差随X增大而增大的数据做线性拟合polyfit和\都给出了“漂亮”的直线和R²。但fitlm的残差图清晰地显示出漏斗形状警告我模型假设不成立。此时盲目使用基于同方差假设的标准误和p值会导致错误推断。解决方案可能是对Y变量进行变换如取对数或使用加权最小二乘法。这个例子说明只看拟合方程和R²是远远不够的模型诊断至关重要。5. 进阶话题与常见陷阱掌握了基础线性拟合后我们可以应对更复杂的场景。5.1 加权最小二乘法当不同数据点的测量精度不同时我们应该给更精确的数据点更高的权重。这就是加权最小二乘。假设我们有一个权重向量w对应的加权最小二乘解为β_wls (XᵀWX)⁻¹ XᵀWY其中W是以w为对角元素的对角阵。 在Matlab中可以用lscov函数轻松实现weights [1, 1, 1, 0.5, 0.5, 0.3]; % 假设后三个数据点测量不确定性更大 beta_wls lscov(X, Y, weights);fitlm也支持通过‘Weights’参数指定权重。5.2 处理带有约束的拟合有时我们的模型参数需要满足某些约束。例如物理上电阻率在绝对零度时应趋近于一个本征值ρ0我们可以强制截距为正值。这是一个线性等式约束可以表示为Aeq * β beq。对于更复杂的约束需要使用优化工具箱的lsqlin函数。% 例如约束截距 0.8 A [-1, 0]; % -β1 -0.8 等价于 β1 0.8 b -0.8; beta_con lsqlin(X, Y, A, b); % 在约束下求解最小二乘5.3 模型选择与过拟合用高阶多项式去拟合数据R²会越来越高但模型会变得毫无预测能力这就是过拟合。对于我们的6个数据点用5次多项式可以完美穿过所有点R²1但这显然是荒谬的。p5 polyfit(T, rho, 5); rho_fit5 polyval(p5, T); R2_5 1 - sum((rho - rho_fit5).^2) / sum((rho - mean(rho)).^2);如何避免看调整R²fitlm输出的Rsquared.Adjusted考虑了参数个数惩罚复杂模型。交叉验证将数据分成训练集和测试集用训练集拟合用测试集计算预测误差。观察残差一个好的模型其残差应该是随机的。如果残差呈现明显的系统性模式说明模型缺失了关键结构。赤池信息准则使用aicbic函数比较不同模型的AIC或BIC值值越小越好。5.4 与相关热词“ttest”和“ttest2”的联系在回归分析中我们对系数进行的t检验fitlm输出的tStat和pValue其原理与热词中提到的ttest单样本t检验和ttest2双样本t检验同源。回归中的t检验是检验某个自变量的系数是否显著不为零即该变量是否对因变量有解释力。而ttest是检验单个样本的均值是否等于某个理论值ttest2是检验两个独立样本的均值是否相等。它们都是基于t统计量的假设检验只是应用场景和原假设不同。理解这一点就能把回归分析和基础的差异检验联系起来。6. 从拟合到应用预测与置信区间拟合模型的最终目的是预测。Matlab可以方便地计算预测值及其置信区间。% 使用 fitlm 模型 mdl 进行预测 new_T [25; 150]; % 想要预测的温度点包括一个外推点(150℃) [rho_pred, rho_ci] predict(mdl, new_T); % 默认返回95%置信区间 fprintf(在T25℃时预测电阻率 %.4f95%% CI: [%.4f, %.4f]\n, rho_pred(1), rho_ci(1,1), rho_ci(1,2)); fprintf(在T150℃时预测电阻率 %.4f95%% CI: [%.4f, %.4f]\n, rho_pred(2), rho_ci(2,1), rho_ci(2,2)); % 绘制拟合线及置信区间带 figure; plot(mdl); hold on; plot(new_T, rho_pred, r*, MarkerSize, 10); % 标出预测点 xlabel(温度 T (℃)); ylabel(电阻率 \rho (\Omega\cdot m)); title(线性回归拟合与预测区间);重要提示从输出可以看到对150℃的预测区间远宽于25℃。这是因为外推预测的不确定性远大于内插预测。我们的模型只在20-120℃的数据范围内得到验证将其应用于范围之外需格外谨慎置信区间会急剧变宽以反映这种不确定性。永远不要轻信模型在外推区域的结果。绕了一大圈我们回到了最初那个电阻率拟合的问题。现在我们不仅能用Matlab画出一条拟合直线更能说清楚这条线是怎么来的矩阵运算与数值分解它的可靠性如何统计检验与诊断以及用它做预测时需要注意什么置信区间与外推风险。最小二乘法不是点一下鼠标就完事的魔法而是一套从模型构建、参数估计、到统计推断的完整方法论。Matlab提供了从快捷操作到底层实现的完整工具链理解每一层工具背后的原理才能让你从“会用软件”变成“会做分析”。下次当你再看到散点图时希望你能本能地思考该用什么模型设计矩阵如何构造残差是否随机系数是否显著这才是数据工作者的核心素养。