简介这份教学课件围绕线性参数的最小二乘法处理展开面向误差理论、测量数据处理相关课程的学习者也可供需要掌握数据拟合方法的工程技术人员参考。内容先阐明最小二乘原理——通过使残差平方和最小化从带随机误差的测量数据中找出最可信赖的参数值随后系统讲解线性测量方程组、残差方程组与正规方程组的建立和求解讨论等权与不等权下的标准差估计及相关系数计算。课件通过例5-1演示电容器精密测定中电容量和标准偏差的完整求解流程例5-2则展示新增串联电容测量后非线性参数的处理方式包括泰勒展开线性化、线性化残差方程组构建及迭代收敛过程。此外还涉及组合测量的基本概念帮助理解多未知量联合估计问题。资源为1个pptx教案文件约306KB共24页章节划分清晰重点难点标注明确。已有60人学习适合课堂教学辅助、课后巩固及考前系统复习使用。1. 最小二乘法公式背得再多线性参数模型才是这个处理方案的真正主角线性参数的最小二乘法处理这个标题很容易让人误以为只讲直线拟合。实际上它覆盖的是更大一类问题只要待估参数以一次方形式进入模型最小二乘法公式就能写成统一的矩阵形式一次解出所有系数。这类学习教案通常会把推导讲得很细但真正做传感器标定、热电偶分度表拟合、光谱校正或应变标定时你缺的往往不是公式而是公式背后的数值稳定性、权处理和残差检验。这篇就按工程落地顺序把原理、代码、参数和坑一次讲全新手可以照着跑熟手可以直接跳到第 3 章以后看边界条件。2. 最小二乘法原理与正规方程从偏导为零到能跑的 Python 代码2.1 线性参数并不等价于直线先弄清楚待求参数怎么进入模型线性参数指的是模型关于待求参数线性不是指曲线必须是一条直线。例如 y a bx 关于 a、b 线性y a·x b·cos(x) 关于 a、b 也是线性因为 a 和 b 都只出现一次不乘对方不取指数不进三角函数。把 x 和 cos(x) 看成两个特征列模型就能写成 y A·θ其中 A 的每一列是一个基函数θ 是待求系数向量。反过来y a·e^(bx) 虽然看起来简单但 b 在指数里求偏导后得到的方程组关于 b 是非线性的不能用下面这套正规方程直接解。判断一个拟合问题能不能用线性最小二乘法就看目标函数对每个待求参数的偏导里是否还含有该参数本身。这份讲义如果只给了 y a bx 的推导你在实际套用时至少要再多问自己一句我的模型能写成 y Aθ 吗能就继续往下走不能就换非线性拟合工具。常见的线性参数模型包括多项式拟合 y a0 a1·x a2·x²、多变量线性回归、傅里叶级数拟合等。这些模型看起来形态差异很大落到矩阵里都是同一件事把已知数据组装成设计矩阵 A然后求解 θ。理解了这一层最小二乘法公式就不再是一条直线的专属公式而是一套通用的参数估计框架。2.2 最小二乘法公式为什么偏导为零对应的是最优解以一元直线模型为例目标函数是所有点的残差平方和J(a, b) Σ(yi - a - b·xi)²要求 J 最小对 a 和 b 分别求偏导并令其为零。对 a 偏导Σ(yi - a - b·xi) 0对 b 偏导Σxi·(yi - a - b·xi) 0整理后得到两个方程n·a Σxi·b ΣyiΣxi·a Σxi²·b Σxi·yi这就是最小二乘法的正规方程。写成矩阵形式会更通用(AᵀA)·θ Aᵀy其中 A 的第一列是 1第二列是 xθ [a, b]ᵀ。正规方程的由来是对残差平方和求梯度并置零它给出的解是所有线性无偏估计中方差最小的那一个这是最小二乘法公式在统计学里的核心价值也是它被大量工程场景选为默认算法的原因。需要特别注意的是正规方程在理论上成立数值上却不一定稳妥。AᵀA 的条件数是 A 条件数的平方。当 x 的动态范围很大比如从 0.001 到 100000AᵀA 会变得病态求解时微小舍入误差可能被放大。工程上我几乎不用显式的逆矩阵解而是用 QR 分解或奇异值分解把问题交给数值稳定的求解器处理。教程推导里可以写 θ (AᵀA)⁻¹Aᵀy落地代码里要换成其他解法。2.3 用 Python 把最小二乘公式落地三步能跑的代码先准备一组带噪声的线性数据用 numpy 的最小二乘求解器跑一遍这也是后续所有讨论的基座。import numpy as np # 构造带噪声的线性数据 rng np.random.default_rng(42) x np.linspace(0, 10, 30) y_true 2.5 * x 1.2 y y_true rng.normal(0, 1.0, sizex.shape) # 组装设计矩阵第一列全是 1对应截距第二列是 x对应斜率 A np.vstack([x, np.ones_like(x)]).T # lstsq 内部走 SVD 路径比 inv(A.T A) 数值稳得多 theta, residuals, rank, s np.linalg.lstsq(A, y, rcondNone) slope, intercept theta print(f斜率 {slope:.4f}, 截距 {intercept:.4f})这段代码有三个关键点。第一设计矩阵 A 的列顺序决定了参数顺序第一列是常数项第二列是 x因此 theta[0] 是斜率theta[1] 是截距。第二rcondNone 表示小于最大奇异值乘以机器精度的奇异值被截断这是 numpy 的默认推荐行为不要改为 0否则病态矩阵会直接报错或给出天文数字。第三np.linalg.lstsq 返回的 rank 和 s 能帮你判断矩阵是否满秩如果 rank 小于 A 的列数说明基函数之间存在线性相关后面会专门讲这个坑。再看一个对照写法方便你理解正规方程和 SVD 的差别# 正规方程适合教学演示实际使用注意条件数 A_T_A A.T A A_T_y A.T y theta_ne np.linalg.solve(A_T_A, A_T_y) # SVD 手动实现看清楚 lstsq 内部在做什么 u, s, vt np.linalg.svd(A, full_matricesFalse) theta_svd vt.T (np.linalg.inv(np.diag(s)) (u.T y)) print(f正规方程 {theta_ne}) print(fSVD 手动版 {theta_svd})两种写法在 x 范围不大时结果几乎一致但一旦 x 从 0 到 1e6正规方程的舍入误差会比 SVD 大好几个量级。我的习惯是只要特征列的量纲差超过三个数量级一律用 lstsq不要让 inv 出现在正式代码里。3. 加权最小二乘与数据预处理测量精度不同时OLS 会翻车3.1 权重矩阵的物理含义最小二乘里可以给每个点不同的信任票普通最小二乘假设每个点的测量误差方差相同但工程数据几乎都违反这个假设。比如用不同量程的传感器拼接测一条曲线低量程段精度高高量程段噪声大或者同一组实验里某些点重复测量了 5 次另一些只测了 1 次。这种情况下方差小的点应当对拟合结果有更大的发言权。加权最小二乘的做法是给每个残差除以它的标准差 σi目标函数变成J Σ((yi - A·θi) / σi)²写成矩阵形式权重矩阵 W 是对角阵对角线元素为 1/σi²θ_w (Aᵀ·W·A)⁻¹ · Aᵀ·W·y权重矩阵不是玄学它来自你对每个数据点不确定度的先验估计。σi 可以是传感器手册给出的重复性指标也可以是同一条件多次测量的样本标准差。没有这些信息时可以先做一次普通最小二乘用残差分布估计方差再重新加权这就是后面要讲的迭代重加权。这里就要提到这类讲义和真实工程的差距了。PPT 里推导最小二乘法原理时通常默认所有点等精度因为那会让公式简洁很多。但你在现场拟合标定数据时不等精度是常态。如果死活不用权重那些噪声特别大的点会像锚一样把拟合曲线拖向它们结果就是斜率偏低、残差集中在低精度段。想避免这种翻车第一步就是检查你的数据点来源凡是有明确不同精度分组的先建权重。3.2 必调参数条件数、归一化与奇异值截断怎么看填入代码之前先学会看三个诊断量。第一个是条件数用 np.linalg.cond(A) 或 np.linalg.cond(A.T A) 计算。条件数越大输入数据的微小扰动对解的影响越大。实战经验大致可以这样判断条件数在 100 以内是良态问题从 100 到 1e6 要警惕噪声放大效应超过 1e12 基本可以认为解已经不可信这时候拟合出的系数改一位小数都会剧烈跳动。条件数范围状态建议 100良态常规求解即可100 ~ 1e6中等病态优先用 lstsq考虑归一化 1e12严重病态必须归一化或换基函数第二个是奇异值数组 s。np.linalg.lstsq 返回的 s 是按从大到小排列的奇异值。如果最后几个奇异值比最大的小 8 个量级以上说明设计矩阵里的特征列几乎是线性相关的解虽然能出来但极不稳定。此时不要急着增大正则项先检查基函数是否重复比如同时放入了 x 和 2x 这样的列。第三个是归一化。处理大动态范围 x 数据的常见做法是把每一列特征减去均值再除以标准差求解完成后再把系数换算回原始尺度。这么做不是因为算法要求而是因为浮点数在 1e-6 和 1e6 混合运算时有效位数会被浪费归一化能给求解器留出更多有效数字空间。注意归一化不是只针对 x多项式拟合里 x²、x³ 列也要一起处理否则条件数依然爆炸。3.3 加权最小二乘落地代码从 OLS 到 WLS 只改一步接第 2.3 节的数据人为构造不同的测量精度然后用权重矩阵重新求解。# 假设每个点的测量标准差已知前 15 个点精度高后 15 个点噪声大 sigma np.where(x 5, 0.3, 2.0) W np.diag(1.0 / sigma**2) # 用 solve 解加权正规方程A.T W A 是实对称正定矩阵 theta_wls np.linalg.solve(A.T W A, A.T W y) # 对比普通最小二乘结果 theta_ols np.linalg.lstsq(A, y, rcondNone)[0] print(fOLS : {theta_ols}) print(fWLS : {theta_wls})这段代码的核心只有一处改动在 AᵀA 和 Aᵀy 中间插入了 W。权重矩阵 W 的对角线是 1/σ²σ 越小权重越大这符合物理直觉。需要说明的是W 的数据类型是 float64 就够了不要用 float32否则系数出来了最后几位全是噪声。如果你没有每个点的 σ 先验值可以用迭代重加权最小二乘。做法是先算普通最小二乘残差然后用 Huber 函数把残差大的点降权再重新拟合重复两三次。下面是一个精简版实现def huber_weight(r, c1.345): # 残差绝对值小于 c 的点权重为 1大于 c 的点按 c/|r| 降权 return np.where(np.abs(r) c, 1.0, c / np.abs(r)) theta np.linalg.lstsq(A, y, rcondNone)[0] for _ in range(3): r y - A theta w huber_weight(r) W np.diag(w) theta np.linalg.solve(A.T W A, A.T W y)迭代重加权对野值的抵抗力比一次性剔除数据点强很多因为它是连续降权而不是硬删保留了部分信息同时也避免了因删除点导致样本量过少的连锁反应。这一招在处理传感器偶发跳变时特别管用也是我在做标定数据处理时最常用的后悔药之一。需要留意的是 c 取 1.345 来自统计学经典经验值如果数据量很小可以适当提高到 2 或 3降权力度更温和。4. 用残差和判定系数验证拟合质量R² 高不一定是模型好4.1 残差是自检报告模式比大小更值得注意残差是 y 减去 A·θ 的结果。很多人只关心 R² 高不高其实残差图才是更早暴露问题的窗口。正规的做法是把残差对 x 画散点图。理想的残差应该在零附近随机分布没有明显的弯曲、喇叭口或分层结构。如果残差呈 U 型说明模型缺了二次项如果残差从小到大张开发散说明数据存在异方差对应的是本文第 3 章的加权需求如果残差在某个区间突然全部同号那多半是那个区间存在系统误差比如仪器在这个量程段没校准。残差的绝对大小也有意义。计算残差标准差 s sqrt(Σr²/(n-k))其中 n 是样本数k 是参数个数。这个 s 是对测量噪声 σ 的估计可以用于后续计算参数标准误。最小二乘算法本身只给出一组系数不给系数的不确定度你需要靠残差标准差和设计矩阵去补参数协方差矩阵近似为 s²·(AᵀA)⁻¹。这一步在标定证书里尤其重要客户要的不只是“斜率是多少”而是“斜率的 95% 置信区间是多少”。教学中可以只讲系数求解但工程报告里没有不确定度评估的拟合结果基本等于半成品。4.2 判定系数 R² 的边界什么情况下它骗人R² 1 - RSS / TSS其中 RSS 是残差平方和TSS 是 y 相对均值的平方和。R² 越接近 1说明模型解释的方差比例越大。但这里有一个常被忽略的陷阱R² 只在线性最小二乘、且模型含截距项时才有完整的方差分解含义。对于加权最小二乘要使用加权后的 R²对于非线性模型R² 的分子分母含义会变得模糊不能直接和线性拟合结果比大小。R² 高并不代表模型正确有几个典型场景会骗人。样本量很少时比如只有 5 个点拟合 4 参数多项式R² 几乎必然接近 1但这只是过拟合的假象x 取值范围很窄时即使 y 和 x 之间是明显曲线关系直线拟合的 R² 也可能高达 0.98数据包含系统误差时比如时间序列存在漂移但没有被建模R² 依然可以很高残差图却能清楚看到非随机分布。因此我的判断顺序是先看残差图再算参数置信区间最后才看 R²。R² 适合用来比较同一批数据上不同候选模型的解释力不适合用来证明一个模型是否物理上成立。针对参数数量多的模型还应该看自由度修正的 R²adj_R² 1 - (1 - R²) · (n - 1) / (n - k - 1)修正后的 R² 会在增加无用参数时下降比普通 R² 更能反映模型复杂度带来的惩罚。4.3 留一交叉验证与 AIC用代码确定模型该不该加项当你要决定用 2 阶还是 5 阶多项式时光靠 R² 不够因为每加一项 R² 都只升不降。这时我给你一个不依赖额外工具包的留一交叉验证实现def loo_cv(A, y): n len(y) err 0.0 for i in range(n): # 每次留出一个样本 mask np.ones(n, dtypebool) mask[i] False theta_i, _, _, _ np.linalg.lstsq(A[mask], y[mask], rcondNone) pred A[i] theta_i err (y[i] - pred) ** 2 return err / n # 比较不同阶数的多项式模型 for degree in [1, 2, 3, 5]: A_poly np.vstack([x**d for d in range(degree 1)]).T mse loo_cv(A_poly, y) print(fdegree{degree}, LOO MSE{mse:.4f})留一交叉验证的核心思想是每次拿 n-1 个点拟合预测被留下的那一个点累计所有预测误差。它比直接在训练集上看 R² 更贴近真实泛化能力也是判断阶数选择有没有翻车的硬指标。LOO MSE 最小的阶数通常就是你该选的阶数。另一个常见准则是 AIC 赤池信息量准则用下面的公式计算aic n·ln(rss / n) 2k其中 k 是参数个数。AIC 在 RSS 和参数数量之间取平衡适合在阶数相近时做快速筛选。注意 AIC 只能用于比较同一数据集上的模型跨数据集比较没有意义。留一交叉验证和 AIC 的说法略有不同前者看预测误差后者看信息损失两者结论不一致时我更信任留一交叉验证的结果因为它是直接对预测能力做测试。5. 线性最小二乘避坑5 条血泪经验与排查方法5.1 条件数爆炸解出来的系数十万八千里现象跑完拟合斜率大得离谱加上微小的数据扰动系数符号都能翻转。原因特征列之间量纲差异太大比如 x 在 0.001 附近而 x² 在 1000 附近AᵀA 的条件数轻松突破 1e12浮点舍入误差被无限放大。解决对每一列特征做标准化求解后还原。还原方法设 x_scaled (x - mean) / std拟合得到 a_s、b_s 是对标准化后的 x 的系数原始斜率 a a_s / std截距 b b_s - a_s · mean / std。一个简单的检查方式是打印 np.linalg.cond(A)见到 1e12 以上就立刻停止手写正规方程改用 lstsq 或先归一化。5.2 高阶多项式过拟合边缘振荡与异常大系数现象拟合阶数加到 7 甚至 10训练集上 R² 接近 0.999但拟合曲线在两个端点剧烈上下摆动系数出现几十上百的数值。原因高阶多项式参数太多模型开始记忆噪声而不是学习趋势正交多项式没有铺开时还会引入数值相关。解决先用留一交叉验证选阶数再加正则化。岭回归是最小二乘的直接扩展把目标函数改为 Σ(yi - Aθ)² λ‖θ‖²λ 用交叉验证选择。低于 5 阶的多项式通常不会碰到这个坑阶数越高越要克制。5.3 指数数据硬套线性模型残差呈 U 型现象数据本身是指数增长却用 y a bx 去拟合残差图呈明显的碗形拟合值在两端都低于观测值。原因模型形式错误把非线性模型误判成线性参数模型此时最小二乘公式再精确也救不了。解决先画散点图看趋势指数型数据换用取对数再线性化或者直接用非线性最小二乘。这里要特别注意一个误用解耦取对数后 y ln(y)模型变为 y ln(A) b·x这时候待求参数是 ln(A) 和 b它们关于新因变量 y 是线性的但拟合得到的截距要取指数还原成 A。还原过程会放大误差所以 A 的置信区间需要用误差传递公式重新算不要直接把拟合截距的标准误当成 A 的标准误。5.4 强制截距为零斜率的额外偏差现象明明数据在 x0 时 y 不为零强行设置过零点后斜率明显偏离正常值残差在零点附近同号。原因线性回归里截距项实际上吸收了数据的基准偏移强制截距为零等于对模型施加了错误的约束所有偏差都被斜率项承接。解决不要轻易删除截距列。先保留截距拟合再看截距的置信区间是否包含零。若置信区间窄且包含零才考虑去掉截距项重新拟合若区间不包含零说明截距显著不能删。这个坑在化学分析的标准曲线数据里特别常见因为很多人误以为空白为零就代表截距必须是零实际上仪器的基线漂移和基质效应都会产生非零信号。5.5 忽略自变量误差最小二乘解向零收缩现象同一个数据集把 x 和 y 互换后拟合得到的斜率乘积不等于 1。原因普通最小二乘只考虑因变量 y 的噪声如果 x 也有不可忽略的误差斜率估计会偏向零这叫回归稀释效应。解决当 x 是人为设定的标准值且误差可忽略时常规最小二乘没问题当 x 和 y 都来自测量且有相当噪声时改用正交回归。正交回归的实现在 numpy 里只需一次 SVDdef total_least_squares(x, y): A np.vstack([x - x.mean(), y - y.mean()]).T _, _, vt np.linalg.svd(A, full_matricesFalse) v vt[-1] # 最小奇异值对应的方向即拟合直线的法向量 slope -v[0] / v[1] intercept y.mean() - slope * x.mean() return slope, intercept这段代码背后的思路是不再只最小化竖向残差而是最小化点到直线的垂直距离。x 和 y 噪声相近时正交回归才是更合理的估计。千万别拿普通最小二乘的结果硬解释成“x 和 y 的斜率关系”那是两个不同问题的答案。6. 从静态拟合到递推最小二乘在线标定与最后一条验证技巧6.1 递推最小二乘遗忘因子 λ 如何调工程里经常遇到数据是分批到的比如设备每天产出新标定点不想把历史数据全部重算一遍这时可以用递推最小二乘 RLS。它的核心是保存参数向量 θ 和协方差矩阵 P每来一个新数据点就更新一次K P_old·x_new / (λ x_newᵀ·P_old·x_new)θ_new θ_old K·(y_new - x_newᵀ·θ_old)P_new (I - K·x_newᵀ)·P_old / λλ 在 0 到 1 之间叫遗忘因子。λ1 表示所有历史数据同等重要适合静态系统λ0.98~0.99 用于参数缓慢漂移的场合比如温度变化造成的传感器灵敏度衰减λ0.9 适应更快但对噪声更敏感。下面是一段可直接运行的 RLS 更新函数def rls_update(theta, P, x_new, y_new, lam0.98): # x_new 是当前点的设计向量比如 [x, 1] x np.asarray(x_new, dtypefloat) # 新息预测误差 innov y_new - float(x theta) # 增益向量 K (P x) / (lam float(x P x)) theta_new theta K * innov P_new (P - np.outer(K, x) P) / lam return theta_new, P_new注意初始化 P 为一个较大的单位阵乘以一个常数比如 1000·I表示初始时对参数一无所知。每轮调用 rls_update 之前先把 x_new 按和原始拟合一样的预处理方式标准化否则第 5 章的条件数雷区会在递推版本里重新引爆。我一般会在递推第 10 个点之后才开始信任输出前几个点的参数跳动属于正常收敛过程。6.2 我的验证习惯合成数据反演与参数还原写完最小二乘代码后我习惯先做一次合成数据反演测试设定一组已知 θ生成无噪声的 y再叠加高斯噪声跑完算法后检查估计值是否落在真实值的 3 倍标准误以内。这个测试能一次性暴露三件事设计矩阵组装有没有错、预处理还原逻辑有没有漏、求解器有没有数值问题。如果估计值和真值偏差超过 3 倍标准误不要去调算法先检查代码。另一个容易忽略的细节是中间结果精度。拟合计算里全程用 float64不要因为显示方便把 x 或矩阵打印出来手动重新录入任何一次手工抄写都等于把有效位数砍到 6 位。保存结果时多用几个有效数字只在最终报告阶段才四舍五入。我现在拿到一批数据的第一件事不是求解而是先画散点图、看条件数、看残差模式这三步做完再谈系数。这个习惯帮我挡掉了至少三成无意义的拟合事故。希望帮到你。本文还有配套的精品资源点击获取
