1. 从“差不多”到“刚刚好”为什么我们需要曲线拟合做数据分析、信号处理或者搞工程建模的朋友肯定都遇到过这种场景你手里有一堆实验数据点横七竖八地散落在坐标图上像一群迷路的蚂蚁。你心里清楚这些数据背后应该藏着某种规律比如一个线性关系、一个指数衰减或者一个更复杂的函数。你的任务就是找到一根“线”能最好地穿过这些点或者至少能最贴切地描述它们背后的趋势。这根“线”就是我们要的“拟合曲线”。这个过程就叫曲线拟合。它绝不是简单地画一根看起来顺眼的线而是一个有严格数学定义的优化过程。那么问题来了什么叫“最好”什么叫“最贴切”对于同一组数据我画一根高一点的线你说离某些点远了我画一根低一点的线他又说离另一些点远了。公说公有理婆说婆有理我们需要一个客观、统一的评判标准。这个标准就是最小二乘原理。它可以说是数据拟合领域里最经典、最基础也几乎是应用最广泛的方法。它的核心思想非常直观甚至有点“朴素”我找的那条曲线应该让所有数据点到这条曲线的垂直距离的平方和最小。为什么是“平方和”而不是直接的距离和这里有两个关键原因。第一距离有正有负直接相加可能会相互抵消比如一个点在曲线上方正距离一个点在下方负距离一加和可能变成0这显然不能反映整体的偏离程度。第二使用平方可以放大那些偏离较远的点的影响让拟合曲线对“异常点”更敏感同时也让整个数学问题变得“友好”——平方函数是光滑可导的这为我们后续用微积分工具求解极值点铺平了道路。所以最小二乘拟合本质上就是在解决一个最优化问题在某一类候选曲线比如所有一次函数y ax b中寻找一组特定的参数a和b使得根据这组参数计算出的预测值与所有实际观测值之间的误差平方和达到全局最小。这个思想由大名鼎鼎的高斯和勒让德在两百多年前分别独立提出并应用于天体轨道计算至今仍是工程和科学研究的基石。2. 最小二乘的数学内核误差、模型与目标函数要真正理解最小二乘我们不能只停留在“让距离平方和最小”这个口号上必须拆开看看它的数学骨架。这个过程也是我们建立任何定量分析模型的标准思路。2.1 定义误差观测值与预测值的差距假设我们有一组观测数据共有n个点。第i个点的坐标是(x_i, y_i)。这里x_i是自变量比如时间、温度、压力y_i是因变量比如位移、电阻、产量是我们实际测量得到的值。现在我们猜测y和x之间存在某种函数关系我们用一个带参数的函数f(x; θ)来表示它。这里的θ代表了一组待确定的参数。例如对于线性模型f(x; a, b) a * x b那么参数θ就是[a, b]。对于每一个数据点(x_i, y_i)我们用模型预测出的y值是ŷ_i f(x_i; θ)。那么这个点的误差或称残差e_i就定义为e_i y_i - ŷ_i y_i - f(x_i; θ)这个误差e_i可正可负代表了模型在该点预测的偏差。2.2 构建目标函数误差平方和SSE如果只有一个误差我们很容易判断模型好坏。但现在有n个误差e_1, e_2, ..., e_n。我们需要一个单一的标量来整体衡量模型在所有数据点上的表现。最小二乘法选择的衡量标准就是误差平方和Sum of Squared Errors, SSES(θ) Σ_{i1}^{n} e_i^2 Σ_{i1}^{n} [y_i - f(x_i; θ)]^2这个S(θ)就是我们的目标函数。它依赖于模型参数θ。不同的参数θ会给出不同的预测值ŷ_i从而计算出不同的S(θ)。2.3 优化求解寻找最优参数最小二乘拟合的目标就是找到一组参数θ*使得目标函数S(θ)的值达到最小θ* argmin_{θ} S(θ)这是一个标准的无约束优化问题。对于许多常见的模型特别是线性模型我们可以通过数学方法直接求出这个最优解的解析表达式。为什么这个方法如此强大从统计学的视角看最小二乘估计在满足一系列理想假设如误差独立、同方差、均值为零时是所有无偏估计中方差最小的即最佳线性无偏估计。从计算的角度看平方项导致目标函数是凸函数对于线性模型这意味着通常只有一个全局最小值我们可以稳定地求解。从感性的角度看它惩罚大的误差比惩罚小的误差更严厉这符合我们“不希望模型在任何一点上偏离太远”的直觉。3. 线性最小二乘手把手推导与几何意义线性拟合是最简单也最常用的情况。我们的模型是y a * x b。此时参数θ [a, b]目标函数为S(a, b) Σ (y_i - (a*x_i b))^2我们的任务是找到a和b最小化S(a, b)。根据微积分函数在极值点处对各个自变量的偏导数应为零。这引出了著名的正规方程组。3.1 正规方程组的推导分别对a和b求偏导并令其等于0∂S/∂a -2 * Σ [x_i * (y_i - a*x_i - b)] 0∂S/∂b -2 * Σ (y_i - a*x_i - b) 0整理后得到正规方程组a * Σ x_i^2 b * Σ x_i Σ (x_i * y_i) a * Σ x_i b * n Σ y_i这是一个关于a和b的二元一次方程组。解这个方程组就得到了最小二乘估计值a (n * Σ(x_i*y_i) - Σx_i * Σy_i) / (n * Σ(x_i^2) - (Σx_i)^2) b (Σy_i * Σ(x_i^2) - Σx_i * Σ(x_i*y_i)) / (n * Σ(x_i^2) - (Σx_i)^2)或者用均值的概念表示会更简洁。令x̄ (Σx_i)/n,ȳ (Σy_i)/n则a Σ[(x_i - x̄)(y_i - ȳ)] / Σ[(x_i - x̄)^2] b ȳ - a * x̄这个形式揭示了a实际上是x和y的协方差除以x的方差非常直观。3.2 线性最小二乘的几何透视我们可以从线性代数的角度获得更深刻的理解。把n个数据点的y值写成一个列向量Y [y_1, y_2, ..., y_n]^T。 我们的线性模型ŷ_i a*x_i b可以改写为Ŷ a * X b * 1其中X [x_1, x_2, ..., x_n]^T1是一个全为1的n维列向量。令矩阵A [X, 1]参数向量β [a, b]^T。那么模型可以写成矩阵形式Ŷ A β我们的目标是最小化真实向量Y与预测向量Ŷ之间的欧氏距离的平方即||Y - Aβ||^2。在几何上Aβ代表的是由矩阵A的列向量即X和1所张成的列空间中的一个向量。Y是空间中的一个点。最小二乘解β*所对应的Ŷ* Aβ*正是Y在这个列空间上的正交投影。注意这个几何解释是理解最小二乘乃至更广义线性模型的关键。它意味着我们寻找的拟合直线其预测值向量Ŷ是真实数据向量Y在由“自变量”和“常数项”所构成平面上的“影子”。误差向量e Y - Ŷ垂直于这个平面。这也解释了为什么正规方程(A^T A) β A^T Y成立——它本质上是要求误差向量与列空间的所有基向量即A的列都正交。4. 超越线性非线性最小二乘与模型线性化现实世界的关系远非总是线性的。可能是指数增长y a * e^{bx}可能是幂律关系y a * x^b也可能是多项式y a_0 a_1*x a_2*x^2 ...。这些都属于非线性最小二乘问题因为待估参数θ与预测值f(x; θ)之间的关系是非线性的。4.1 多项式拟合披着非线性外衣的线性问题多项式拟合y Σ_{j0}^{m} a_j * x^j是一个特例。虽然y关于x是非线性的但关于参数a_j却是线性的我们可以令φ_j(x) x^j那么模型就变成了y a_0*φ_0(x) a_1*φ_1(x) ... a_m*φ_m(x)。这被称为线性于参数的模型。对于这类模型我们完全可以套用线性最小二乘的框架。只需将设计矩阵A中的列从[X, 1]扩展为[1, X, X.^2, ..., X.^m]其中X.^k表示对向量X的每个元素求k次幂。然后求解正规方程(A^T A) β A^T Y其中β [a_0, a_1, ..., a_m]^T即可。在 MATLAB 或 Python (NumPy) 中这通常只需一两行代码。实操心得多项式阶数m的选择这里有一个经典的陷阱过拟合。阶数m越高曲线越“柔软”能更精确地穿过每一个数据点甚至让 SSE 降为 0但这样的曲线往往震荡剧烈失去了揭示底层规律的能力对新数据的预测能力极差。我个人的经验法则是先可视化画出散点图观察大致趋势。线性二次饱和增长从低阶开始优先尝试m1线性m2二次。很多时候简单的模型更稳健。交叉验证如果有足够数据将数据分为训练集和测试集。用训练集拟合不同阶数的模型在测试集上计算误差。选择测试集误差最小的模型。观察系数如果高阶项的系数绝对值非常小或者其置信区间包含0通常可以考虑去掉该项。一个实用警告尽量不要让多项式阶数m超过数据点数量n的十分之一对于m接近n的情况结果基本是灾难性的。4.2 真正的非线性拟合迭代优化与线性化技巧对于像y a * e^{bx}或y a / (b x)这类参数非线性的模型目标函数S(θ)关于θ是非线性的可能有很多局部极小值。我们无法直接解出像正规方程那样的解析解必须借助迭代优化算法例如高斯-牛顿法一种利用目标函数近似为二次型的迭代方法收敛速度快但需要计算雅可比矩阵且对初始值敏感。列文伯格-马夸尔特法高斯-牛顿法的改进版通过引入阻尼因子在梯度下降和高斯-牛顿法之间自适应切换更稳定、更常用。SciPy 和 MATLAB 中的lsqcurvefit、curve_fit等函数默认或常用此算法。信任域反射法另一种稳健的迭代算法。实操中的关键参数初始值对于非线性拟合提供一个好的参数初始猜测θ0至关重要它直接决定了算法能否收敛到全局最优以及收敛的速度。我常用的策略是物理意义法如果参数有物理意义如衰减率、饱和值根据数据范围和经验给出粗略估计。线性化法这是最实用的技巧之一。以指数模型y a * e^{bx}为例两边取自然对数ln(y) ln(a) b*x。令Y ln(y)A ln(a) 则模型变为Y A b*x 这是一个关于A和b的线性模型我们可以先用线性最小二乘拟合(x, ln(y))的数据得到A和b的初始估计然后反推出a exp(A)。这个方法对于幂律、饱和增长等许多可线性化模型都有效。网格搜索法如果参数范围大致可知可以在一个粗糙的网格上计算S(θ)选择使S(θ)最小的点作为初始值。5. 实战避坑从MATLAB到Python的完整流程与陷阱理论懂了不实践等于零。我们以 MATLAB 和 Python (SciPy) 为例走一遍完整的拟合流程并指出每一步可能遇到的坑。5.1 数据准备与可视化第一步就错了后面全白费拿到数据千万别急着fit或curve_fit。第一步永远是可视化。% MATLAB 示例 % 假设已有数据向量 x_data 和 y_data figure; scatter(x_data, y_data, 50, filled, DisplayName, 原始数据); xlabel(自变量 X); ylabel(因变量 Y); title(数据散点图); grid on; legend; % 仔细观察点的分布是线性有弯曲有异常点# Python (Matplotlib) 示例 import matplotlib.pyplot as plt import numpy as np # 假设已有数据数组 x_data 和 y_data plt.figure(figsize(8, 6)) plt.scatter(x_data, y_data, s50, label原始数据, alpha0.7) plt.xlabel(自变量 X) plt.ylabel(因变量 Y) plt.title(数据散点图) plt.grid(True, linestyle--, alpha0.5) plt.legend() plt.show()这一步的陷阱异常值图中是否有一两个点离群索居它们会严重扭曲最小二乘的结果因为平方项放大了大误差的影响。需要判断是测量错误应剔除还是重要现象需用稳健回归等方法。数据密度不均如果数据在某个x区间非常密集在另一个区间非常稀疏最小二乘拟合的结果会被密集区域“主导”。需要考虑是否进行数据加权或变换。异方差性误差的方差是否随着x变化如果数据在x较大时波动范围明显变大这就违反了“同方差”假设标准最小二乘的统计性质会变差。5.2 模型选择与拟合执行根据可视化结果选择模型。假设我们决定用二次多项式y p0 p1*x p2*x^2拟合。% MATLAB 多项式拟合 (二次) p polyfit(x_data, y_data, 2); % p 是系数向量从高次到低次: p(1)*x^2 p(2)*x p(3) % 生成拟合曲线上的点 x_fit linspace(min(x_data), max(x_data), 100); y_fit polyval(p, x_fit); % 画图对比 hold on; plot(x_fit, y_fit, r-, LineWidth, 2, DisplayName, 二次多项式拟合); hold off;# Python 多项式拟合 (使用 NumPy) import numpy as np p np.polyfit(x_data, y_data, deg2) # p 是系数向量从高次到低次: p[0]*x^2 p[1]*x p[2] # 生成拟合曲线上的点 x_fit np.linspace(np.min(x_data), np.max(x_data), 100) y_fit np.polyval(p, x_fit) # 画图对比 (接续之前的画图代码) plt.plot(x_fit, y_fit, r-, linewidth2, label二次多项式拟合) plt.legend() plt.show()对于非线性模型例如指数衰减y a * exp(-b*x) c% MATLAB 非线性拟合 (使用 fit 函数或 lsqcurvefit) % 定义模型函数句柄 modelfun (p, x) p(1) * exp(-p(2)*x) p(3); % 初始猜测值 p0 [a_guess, b_guess, c_guess] p0 [max(y_data), 0.1, min(y_data)]; % 使用 lsqcurvefit options optimoptions(lsqcurvefit, Display, iter); % 显示迭代过程 [p_opt, resnorm] lsqcurvefit(modelfun, p0, x_data, y_data, [], [], options); % p_opt 是最优参数 y_fit_nl modelfun(p_opt, x_fit);# Python 非线性拟合 (使用 SciPy) from scipy.optimize import curve_fit import numpy as np # 定义模型函数 def exp_decay(x, a, b, c): return a * np.exp(-b * x) c # 初始猜测值 p0 [a_guess, b_guess, c_guess] p0 [np.max(y_data), 0.1, np.min(y_data)] # 执行拟合 p_opt, p_cov curve_fit(exp_decay, x_data, y_data, p0p0) # p_opt 是最优参数 p_cov 是参数的协方差矩阵可用于计算标准差 y_fit_nl exp_decay(x_fit, *p_opt)这一步的陷阱初始值导致收敛到局部最优对于非线性拟合如果结果不合理如预测曲线完全偏离数据首先怀疑初始值。尝试不同的p0或使用前面提到的线性化技巧获取初始值。参数物理意义与约束有时参数应有物理范围如衰减率b 0。curve_fit和lsqcurvefit都支持设置参数的上下界 (bounds)。算法不收敛可能模型过于复杂或数据噪声太大。可以尝试简化模型增加最大迭代次数 (maxfev)或调整优化算法参数。5.3 结果评估与诊断拟合好≠模型好拟合出一条曲线只是开始评估其质量才是重点。绝不能只看“拟合曲线和散点图看起来挺近”。1. 残差分析这是诊断模型缺陷的最有力工具。残差e_i y_i - ŷ_i应该随机分布在0附近不应有任何明显的模式。# Python 残差图 y_pred exp_decay(x_data, *p_opt) # 或用 polyval 计算预测值 residuals y_data - y_pred fig, axes plt.subplots(1, 2, figsize(12, 4)) # 残差 vs. 预测值 axes[0].scatter(y_pred, residuals, alpha0.7) axes[0].axhline(y0, colorr, linestyle--) axes[0].set_xlabel(预测值 ŷ) axes[0].set_ylabel(残差) axes[0].set_title(残差 vs. 预测值) axes[0].grid(True, alpha0.3) # 残差 vs. 自变量 axes[1].scatter(x_data, residuals, alpha0.7) axes[1].axhline(y0, colorr, linestyle--) axes[1].set_xlabel(自变量 X) axes[1].set_ylabel(残差) axes[1].set_title(残差 vs. 自变量) axes[1].grid(True, alpha0.3) plt.tight_layout() plt.show()理想的残差图点随机、均匀地分布在水平线y0两侧无明显趋势、漏斗形或弯曲。有问题的残差图漏斗形残差随预测值增大而发散提示异方差需要考虑加权最小二乘或对y做变换如取对数。U型或倒U型曲线残差与预测值/自变量呈系统性弯曲提示模型缺失了重要项如该用二次但用了线性。周期性波动残差呈现周期性提示数据可能存在未被模型捕捉的周期成分。2. 量化指标R-squared (决定系数)表示模型解释的数据变异性的比例。越接近1越好。但注意对于非线性模型其定义和解释与线性模型不同需谨慎使用。调整后的 R-squared考虑了模型复杂度参数个数用于比较不同复杂度模型。均方根误差RMSE sqrt(SSE / n)。它与y有相同量纲可以直观理解为“平均预测误差有多大”。参数的标准误与置信区间从协方差矩阵p_cov可以计算参数的标准差进而评估参数估计的精度。如果某个参数的置信区间包含0意味着该参数可能不显著对于线性模型。3. 过拟合检验针对多项式或复杂模型将数据随机分为训练集如70%和测试集30%。用训练集拟合模型然后计算模型在测试集上的 RMSE。如果训练集 RMSE 远小于测试集 RMSE说明模型过拟合了训练数据的噪声泛化能力差。6. 系统辨识中的应用从数据到动态模型在控制工程和系统辨识领域最小二乘法是构建动态系统数学模型如差分方程、传递函数的基石。这里的“曲线”变成了系统输入输出数据在时间序列上展现的动态关系。假设我们有一个离散时间系统我们认为其输出y(k)与过去若干时刻的输入u(k-1), u(k-2)...和自身过去的输出y(k-1), y(k-2)...有关这可以用一个线性差分方程描述自回归外生模型y(k) a1*y(k-1) ... ana*y(k-na) b1*u(k-1) ... bnb*u(k-nb) e(k)其中e(k)是白噪声。我们可以将其重写为y(k) -a1*y(k-1) - ... - ana*y(k-na) b1*u(k-1) ... bnb*u(k-nb) e(k) φ(k)^T * θ e(k)其中φ(k) [-y(k-1), ..., -y(k-na), u(k-1), ..., u(k-nb)]^T是回归向量θ [a1, ..., ana, b1, ..., bnb]^T是待辨识的参数向量。对于从k1到kN的N组观测数据我们可以构建Y [y(1), y(2), ..., y(N)]^TΦ [φ(1)^T; φ(2)^T; ...; φ(N)^T]这是一个N x (nanb)的矩阵E [e(1), e(2), ..., e(N)]^T于是系统方程可以写成矩阵形式Y Φθ E。这正是一个标准的线性最小二乘问题目标是最小化误差平方和E^T E。其最小二乘解为θ_LS (Φ^T Φ)^{-1} Φ^T Y这个公式与线性回归的正规方程解在形式上完全一致。通过采集系统的输入输出数据构造出Φ和Y我们就可以利用这个公式一次性估计出模型的所有参数θ。在MATLAB中的系统辨识工具箱和Python的 SciPy 或 SysIdentPy 库中都有现成的函数来实现这个过程。例如在较新版本的MATLAB中可以使用arx函数或tfest函数在Python中可以自定义构建Φ矩阵并使用np.linalg.lstsq求解。系统辨识中的特殊考量持续激励输入信号u(k)需要足够“丰富”如包含多种频率成分的伪随机二进制序列才能激励出系统的所有动态模式使矩阵Φ^T Φ可逆满秩。模型阶次选择na和nb的选择至关重要。阶次太低模型欠拟合无法捕捉系统动态阶次太高模型过拟合会去拟合噪声。通常结合残差检验残差是否像白噪声和信息准则如AIC、BIC来综合确定。闭环辨识如果数据来自闭环控制系统直接使用上述最小二乘可能会产生有偏估计需要采用其他方法如辅助变量法、预报误差法等。从静态的曲线拟合到动态的系统辨识最小二乘原理提供了一套统一、强大的框架将我们面对杂乱数据时的直觉——“找一条最贴近所有点的线”——转化为了严谨的数学优化问题。理解其背后的误差定义、目标函数构建和优化求解过程不仅能让你在调用polyfit或curve_fit时心里有底更能让你在面对更复杂的建模任务时知道如何将问题“翻译”成最小二乘能够解决的形式。这或许就是这个古老原理历经两百年而不衰的魅力所在。
1. 从“差不多”到“刚刚好”为什么我们需要曲线拟合做数据分析、信号处理或者搞工程建模的朋友肯定都遇到过这种场景你手里有一堆实验数据点横七竖八地散落在坐标图上像一群迷路的蚂蚁。你心里清楚这些数据背后应该藏着某种规律比如一个线性关系、一个指数衰减或者一个更复杂的函数。你的任务就是找到一根“线”能最好地穿过这些点或者至少能最贴切地描述它们背后的趋势。这根“线”就是我们要的“拟合曲线”。这个过程就叫曲线拟合。它绝不是简单地画一根看起来顺眼的线而是一个有严格数学定义的优化过程。那么问题来了什么叫“最好”什么叫“最贴切”对于同一组数据我画一根高一点的线你说离某些点远了我画一根低一点的线他又说离另一些点远了。公说公有理婆说婆有理我们需要一个客观、统一的评判标准。这个标准就是最小二乘原理。它可以说是数据拟合领域里最经典、最基础也几乎是应用最广泛的方法。它的核心思想非常直观甚至有点“朴素”我找的那条曲线应该让所有数据点到这条曲线的垂直距离的平方和最小。为什么是“平方和”而不是直接的距离和这里有两个关键原因。第一距离有正有负直接相加可能会相互抵消比如一个点在曲线上方正距离一个点在下方负距离一加和可能变成0这显然不能反映整体的偏离程度。第二使用平方可以放大那些偏离较远的点的影响让拟合曲线对“异常点”更敏感同时也让整个数学问题变得“友好”——平方函数是光滑可导的这为我们后续用微积分工具求解极值点铺平了道路。所以最小二乘拟合本质上就是在解决一个最优化问题在某一类候选曲线比如所有一次函数y ax b中寻找一组特定的参数a和b使得根据这组参数计算出的预测值与所有实际观测值之间的误差平方和达到全局最小。这个思想由大名鼎鼎的高斯和勒让德在两百多年前分别独立提出并应用于天体轨道计算至今仍是工程和科学研究的基石。2. 最小二乘的数学内核误差、模型与目标函数要真正理解最小二乘我们不能只停留在“让距离平方和最小”这个口号上必须拆开看看它的数学骨架。这个过程也是我们建立任何定量分析模型的标准思路。2.1 定义误差观测值与预测值的差距假设我们有一组观测数据共有n个点。第i个点的坐标是(x_i, y_i)。这里x_i是自变量比如时间、温度、压力y_i是因变量比如位移、电阻、产量是我们实际测量得到的值。现在我们猜测y和x之间存在某种函数关系我们用一个带参数的函数f(x; θ)来表示它。这里的θ代表了一组待确定的参数。例如对于线性模型f(x; a, b) a * x b那么参数θ就是[a, b]。对于每一个数据点(x_i, y_i)我们用模型预测出的y值是ŷ_i f(x_i; θ)。那么这个点的误差或称残差e_i就定义为e_i y_i - ŷ_i y_i - f(x_i; θ)这个误差e_i可正可负代表了模型在该点预测的偏差。2.2 构建目标函数误差平方和SSE如果只有一个误差我们很容易判断模型好坏。但现在有n个误差e_1, e_2, ..., e_n。我们需要一个单一的标量来整体衡量模型在所有数据点上的表现。最小二乘法选择的衡量标准就是误差平方和Sum of Squared Errors, SSES(θ) Σ_{i1}^{n} e_i^2 Σ_{i1}^{n} [y_i - f(x_i; θ)]^2这个S(θ)就是我们的目标函数。它依赖于模型参数θ。不同的参数θ会给出不同的预测值ŷ_i从而计算出不同的S(θ)。2.3 优化求解寻找最优参数最小二乘拟合的目标就是找到一组参数θ*使得目标函数S(θ)的值达到最小θ* argmin_{θ} S(θ)这是一个标准的无约束优化问题。对于许多常见的模型特别是线性模型我们可以通过数学方法直接求出这个最优解的解析表达式。为什么这个方法如此强大从统计学的视角看最小二乘估计在满足一系列理想假设如误差独立、同方差、均值为零时是所有无偏估计中方差最小的即最佳线性无偏估计。从计算的角度看平方项导致目标函数是凸函数对于线性模型这意味着通常只有一个全局最小值我们可以稳定地求解。从感性的角度看它惩罚大的误差比惩罚小的误差更严厉这符合我们“不希望模型在任何一点上偏离太远”的直觉。3. 线性最小二乘手把手推导与几何意义线性拟合是最简单也最常用的情况。我们的模型是y a * x b。此时参数θ [a, b]目标函数为S(a, b) Σ (y_i - (a*x_i b))^2我们的任务是找到a和b最小化S(a, b)。根据微积分函数在极值点处对各个自变量的偏导数应为零。这引出了著名的正规方程组。3.1 正规方程组的推导分别对a和b求偏导并令其等于0∂S/∂a -2 * Σ [x_i * (y_i - a*x_i - b)] 0∂S/∂b -2 * Σ (y_i - a*x_i - b) 0整理后得到正规方程组a * Σ x_i^2 b * Σ x_i Σ (x_i * y_i) a * Σ x_i b * n Σ y_i这是一个关于a和b的二元一次方程组。解这个方程组就得到了最小二乘估计值a (n * Σ(x_i*y_i) - Σx_i * Σy_i) / (n * Σ(x_i^2) - (Σx_i)^2) b (Σy_i * Σ(x_i^2) - Σx_i * Σ(x_i*y_i)) / (n * Σ(x_i^2) - (Σx_i)^2)或者用均值的概念表示会更简洁。令x̄ (Σx_i)/n,ȳ (Σy_i)/n则a Σ[(x_i - x̄)(y_i - ȳ)] / Σ[(x_i - x̄)^2] b ȳ - a * x̄这个形式揭示了a实际上是x和y的协方差除以x的方差非常直观。3.2 线性最小二乘的几何透视我们可以从线性代数的角度获得更深刻的理解。把n个数据点的y值写成一个列向量Y [y_1, y_2, ..., y_n]^T。 我们的线性模型ŷ_i a*x_i b可以改写为Ŷ a * X b * 1其中X [x_1, x_2, ..., x_n]^T1是一个全为1的n维列向量。令矩阵A [X, 1]参数向量β [a, b]^T。那么模型可以写成矩阵形式Ŷ A β我们的目标是最小化真实向量Y与预测向量Ŷ之间的欧氏距离的平方即||Y - Aβ||^2。在几何上Aβ代表的是由矩阵A的列向量即X和1所张成的列空间中的一个向量。Y是空间中的一个点。最小二乘解β*所对应的Ŷ* Aβ*正是Y在这个列空间上的正交投影。注意这个几何解释是理解最小二乘乃至更广义线性模型的关键。它意味着我们寻找的拟合直线其预测值向量Ŷ是真实数据向量Y在由“自变量”和“常数项”所构成平面上的“影子”。误差向量e Y - Ŷ垂直于这个平面。这也解释了为什么正规方程(A^T A) β A^T Y成立——它本质上是要求误差向量与列空间的所有基向量即A的列都正交。4. 超越线性非线性最小二乘与模型线性化现实世界的关系远非总是线性的。可能是指数增长y a * e^{bx}可能是幂律关系y a * x^b也可能是多项式y a_0 a_1*x a_2*x^2 ...。这些都属于非线性最小二乘问题因为待估参数θ与预测值f(x; θ)之间的关系是非线性的。4.1 多项式拟合披着非线性外衣的线性问题多项式拟合y Σ_{j0}^{m} a_j * x^j是一个特例。虽然y关于x是非线性的但关于参数a_j却是线性的我们可以令φ_j(x) x^j那么模型就变成了y a_0*φ_0(x) a_1*φ_1(x) ... a_m*φ_m(x)。这被称为线性于参数的模型。对于这类模型我们完全可以套用线性最小二乘的框架。只需将设计矩阵A中的列从[X, 1]扩展为[1, X, X.^2, ..., X.^m]其中X.^k表示对向量X的每个元素求k次幂。然后求解正规方程(A^T A) β A^T Y其中β [a_0, a_1, ..., a_m]^T即可。在 MATLAB 或 Python (NumPy) 中这通常只需一两行代码。实操心得多项式阶数m的选择这里有一个经典的陷阱过拟合。阶数m越高曲线越“柔软”能更精确地穿过每一个数据点甚至让 SSE 降为 0但这样的曲线往往震荡剧烈失去了揭示底层规律的能力对新数据的预测能力极差。我个人的经验法则是先可视化画出散点图观察大致趋势。线性二次饱和增长从低阶开始优先尝试m1线性m2二次。很多时候简单的模型更稳健。交叉验证如果有足够数据将数据分为训练集和测试集。用训练集拟合不同阶数的模型在测试集上计算误差。选择测试集误差最小的模型。观察系数如果高阶项的系数绝对值非常小或者其置信区间包含0通常可以考虑去掉该项。一个实用警告尽量不要让多项式阶数m超过数据点数量n的十分之一对于m接近n的情况结果基本是灾难性的。4.2 真正的非线性拟合迭代优化与线性化技巧对于像y a * e^{bx}或y a / (b x)这类参数非线性的模型目标函数S(θ)关于θ是非线性的可能有很多局部极小值。我们无法直接解出像正规方程那样的解析解必须借助迭代优化算法例如高斯-牛顿法一种利用目标函数近似为二次型的迭代方法收敛速度快但需要计算雅可比矩阵且对初始值敏感。列文伯格-马夸尔特法高斯-牛顿法的改进版通过引入阻尼因子在梯度下降和高斯-牛顿法之间自适应切换更稳定、更常用。SciPy 和 MATLAB 中的lsqcurvefit、curve_fit等函数默认或常用此算法。信任域反射法另一种稳健的迭代算法。实操中的关键参数初始值对于非线性拟合提供一个好的参数初始猜测θ0至关重要它直接决定了算法能否收敛到全局最优以及收敛的速度。我常用的策略是物理意义法如果参数有物理意义如衰减率、饱和值根据数据范围和经验给出粗略估计。线性化法这是最实用的技巧之一。以指数模型y a * e^{bx}为例两边取自然对数ln(y) ln(a) b*x。令Y ln(y)A ln(a) 则模型变为Y A b*x 这是一个关于A和b的线性模型我们可以先用线性最小二乘拟合(x, ln(y))的数据得到A和b的初始估计然后反推出a exp(A)。这个方法对于幂律、饱和增长等许多可线性化模型都有效。网格搜索法如果参数范围大致可知可以在一个粗糙的网格上计算S(θ)选择使S(θ)最小的点作为初始值。5. 实战避坑从MATLAB到Python的完整流程与陷阱理论懂了不实践等于零。我们以 MATLAB 和 Python (SciPy) 为例走一遍完整的拟合流程并指出每一步可能遇到的坑。5.1 数据准备与可视化第一步就错了后面全白费拿到数据千万别急着fit或curve_fit。第一步永远是可视化。% MATLAB 示例 % 假设已有数据向量 x_data 和 y_data figure; scatter(x_data, y_data, 50, filled, DisplayName, 原始数据); xlabel(自变量 X); ylabel(因变量 Y); title(数据散点图); grid on; legend; % 仔细观察点的分布是线性有弯曲有异常点# Python (Matplotlib) 示例 import matplotlib.pyplot as plt import numpy as np # 假设已有数据数组 x_data 和 y_data plt.figure(figsize(8, 6)) plt.scatter(x_data, y_data, s50, label原始数据, alpha0.7) plt.xlabel(自变量 X) plt.ylabel(因变量 Y) plt.title(数据散点图) plt.grid(True, linestyle--, alpha0.5) plt.legend() plt.show()这一步的陷阱异常值图中是否有一两个点离群索居它们会严重扭曲最小二乘的结果因为平方项放大了大误差的影响。需要判断是测量错误应剔除还是重要现象需用稳健回归等方法。数据密度不均如果数据在某个x区间非常密集在另一个区间非常稀疏最小二乘拟合的结果会被密集区域“主导”。需要考虑是否进行数据加权或变换。异方差性误差的方差是否随着x变化如果数据在x较大时波动范围明显变大这就违反了“同方差”假设标准最小二乘的统计性质会变差。5.2 模型选择与拟合执行根据可视化结果选择模型。假设我们决定用二次多项式y p0 p1*x p2*x^2拟合。% MATLAB 多项式拟合 (二次) p polyfit(x_data, y_data, 2); % p 是系数向量从高次到低次: p(1)*x^2 p(2)*x p(3) % 生成拟合曲线上的点 x_fit linspace(min(x_data), max(x_data), 100); y_fit polyval(p, x_fit); % 画图对比 hold on; plot(x_fit, y_fit, r-, LineWidth, 2, DisplayName, 二次多项式拟合); hold off;# Python 多项式拟合 (使用 NumPy) import numpy as np p np.polyfit(x_data, y_data, deg2) # p 是系数向量从高次到低次: p[0]*x^2 p[1]*x p[2] # 生成拟合曲线上的点 x_fit np.linspace(np.min(x_data), np.max(x_data), 100) y_fit np.polyval(p, x_fit) # 画图对比 (接续之前的画图代码) plt.plot(x_fit, y_fit, r-, linewidth2, label二次多项式拟合) plt.legend() plt.show()对于非线性模型例如指数衰减y a * exp(-b*x) c% MATLAB 非线性拟合 (使用 fit 函数或 lsqcurvefit) % 定义模型函数句柄 modelfun (p, x) p(1) * exp(-p(2)*x) p(3); % 初始猜测值 p0 [a_guess, b_guess, c_guess] p0 [max(y_data), 0.1, min(y_data)]; % 使用 lsqcurvefit options optimoptions(lsqcurvefit, Display, iter); % 显示迭代过程 [p_opt, resnorm] lsqcurvefit(modelfun, p0, x_data, y_data, [], [], options); % p_opt 是最优参数 y_fit_nl modelfun(p_opt, x_fit);# Python 非线性拟合 (使用 SciPy) from scipy.optimize import curve_fit import numpy as np # 定义模型函数 def exp_decay(x, a, b, c): return a * np.exp(-b * x) c # 初始猜测值 p0 [a_guess, b_guess, c_guess] p0 [np.max(y_data), 0.1, np.min(y_data)] # 执行拟合 p_opt, p_cov curve_fit(exp_decay, x_data, y_data, p0p0) # p_opt 是最优参数 p_cov 是参数的协方差矩阵可用于计算标准差 y_fit_nl exp_decay(x_fit, *p_opt)这一步的陷阱初始值导致收敛到局部最优对于非线性拟合如果结果不合理如预测曲线完全偏离数据首先怀疑初始值。尝试不同的p0或使用前面提到的线性化技巧获取初始值。参数物理意义与约束有时参数应有物理范围如衰减率b 0。curve_fit和lsqcurvefit都支持设置参数的上下界 (bounds)。算法不收敛可能模型过于复杂或数据噪声太大。可以尝试简化模型增加最大迭代次数 (maxfev)或调整优化算法参数。5.3 结果评估与诊断拟合好≠模型好拟合出一条曲线只是开始评估其质量才是重点。绝不能只看“拟合曲线和散点图看起来挺近”。1. 残差分析这是诊断模型缺陷的最有力工具。残差e_i y_i - ŷ_i应该随机分布在0附近不应有任何明显的模式。# Python 残差图 y_pred exp_decay(x_data, *p_opt) # 或用 polyval 计算预测值 residuals y_data - y_pred fig, axes plt.subplots(1, 2, figsize(12, 4)) # 残差 vs. 预测值 axes[0].scatter(y_pred, residuals, alpha0.7) axes[0].axhline(y0, colorr, linestyle--) axes[0].set_xlabel(预测值 ŷ) axes[0].set_ylabel(残差) axes[0].set_title(残差 vs. 预测值) axes[0].grid(True, alpha0.3) # 残差 vs. 自变量 axes[1].scatter(x_data, residuals, alpha0.7) axes[1].axhline(y0, colorr, linestyle--) axes[1].set_xlabel(自变量 X) axes[1].set_ylabel(残差) axes[1].set_title(残差 vs. 自变量) axes[1].grid(True, alpha0.3) plt.tight_layout() plt.show()理想的残差图点随机、均匀地分布在水平线y0两侧无明显趋势、漏斗形或弯曲。有问题的残差图漏斗形残差随预测值增大而发散提示异方差需要考虑加权最小二乘或对y做变换如取对数。U型或倒U型曲线残差与预测值/自变量呈系统性弯曲提示模型缺失了重要项如该用二次但用了线性。周期性波动残差呈现周期性提示数据可能存在未被模型捕捉的周期成分。2. 量化指标R-squared (决定系数)表示模型解释的数据变异性的比例。越接近1越好。但注意对于非线性模型其定义和解释与线性模型不同需谨慎使用。调整后的 R-squared考虑了模型复杂度参数个数用于比较不同复杂度模型。均方根误差RMSE sqrt(SSE / n)。它与y有相同量纲可以直观理解为“平均预测误差有多大”。参数的标准误与置信区间从协方差矩阵p_cov可以计算参数的标准差进而评估参数估计的精度。如果某个参数的置信区间包含0意味着该参数可能不显著对于线性模型。3. 过拟合检验针对多项式或复杂模型将数据随机分为训练集如70%和测试集30%。用训练集拟合模型然后计算模型在测试集上的 RMSE。如果训练集 RMSE 远小于测试集 RMSE说明模型过拟合了训练数据的噪声泛化能力差。6. 系统辨识中的应用从数据到动态模型在控制工程和系统辨识领域最小二乘法是构建动态系统数学模型如差分方程、传递函数的基石。这里的“曲线”变成了系统输入输出数据在时间序列上展现的动态关系。假设我们有一个离散时间系统我们认为其输出y(k)与过去若干时刻的输入u(k-1), u(k-2)...和自身过去的输出y(k-1), y(k-2)...有关这可以用一个线性差分方程描述自回归外生模型y(k) a1*y(k-1) ... ana*y(k-na) b1*u(k-1) ... bnb*u(k-nb) e(k)其中e(k)是白噪声。我们可以将其重写为y(k) -a1*y(k-1) - ... - ana*y(k-na) b1*u(k-1) ... bnb*u(k-nb) e(k) φ(k)^T * θ e(k)其中φ(k) [-y(k-1), ..., -y(k-na), u(k-1), ..., u(k-nb)]^T是回归向量θ [a1, ..., ana, b1, ..., bnb]^T是待辨识的参数向量。对于从k1到kN的N组观测数据我们可以构建Y [y(1), y(2), ..., y(N)]^TΦ [φ(1)^T; φ(2)^T; ...; φ(N)^T]这是一个N x (nanb)的矩阵E [e(1), e(2), ..., e(N)]^T于是系统方程可以写成矩阵形式Y Φθ E。这正是一个标准的线性最小二乘问题目标是最小化误差平方和E^T E。其最小二乘解为θ_LS (Φ^T Φ)^{-1} Φ^T Y这个公式与线性回归的正规方程解在形式上完全一致。通过采集系统的输入输出数据构造出Φ和Y我们就可以利用这个公式一次性估计出模型的所有参数θ。在MATLAB中的系统辨识工具箱和Python的 SciPy 或 SysIdentPy 库中都有现成的函数来实现这个过程。例如在较新版本的MATLAB中可以使用arx函数或tfest函数在Python中可以自定义构建Φ矩阵并使用np.linalg.lstsq求解。系统辨识中的特殊考量持续激励输入信号u(k)需要足够“丰富”如包含多种频率成分的伪随机二进制序列才能激励出系统的所有动态模式使矩阵Φ^T Φ可逆满秩。模型阶次选择na和nb的选择至关重要。阶次太低模型欠拟合无法捕捉系统动态阶次太高模型过拟合会去拟合噪声。通常结合残差检验残差是否像白噪声和信息准则如AIC、BIC来综合确定。闭环辨识如果数据来自闭环控制系统直接使用上述最小二乘可能会产生有偏估计需要采用其他方法如辅助变量法、预报误差法等。从静态的曲线拟合到动态的系统辨识最小二乘原理提供了一套统一、强大的框架将我们面对杂乱数据时的直觉——“找一条最贴近所有点的线”——转化为了严谨的数学优化问题。理解其背后的误差定义、目标函数构建和优化求解过程不仅能让你在调用polyfit或curve_fit时心里有底更能让你在面对更复杂的建模任务时知道如何将问题“翻译”成最小二乘能够解决的形式。这或许就是这个古老原理历经两百年而不衰的魅力所在。
复杂度计算叶子节点与度1/2节点)
)
