简介最小二乘法是数据拟合与回归分析的核心原理其本质是通过最小化误差平方和来确定模型参数。线性模型可通过正规方程直接求解如numpy.polyfit非线性模型则依赖scipy.optimize.curve_fit进行迭代优化。该方法在传感器标定、光谱分析、信号处理等工程场景中应用广泛但实际落地时模型选择、初值设置、权重分配与数值稳定性都会显著影响结果。本文从10个真实案例出发覆盖直线、多项式、对数、指数、幂律、高斯峰、洛伦兹峰、正弦、加权拟合、分段线性与圆拟合演示如何使用python生态工具完成稳健的数据拟合并总结常见坑点与排查思路。 最小二乘这个名词学统计的、做测绘的、搞算法的、写硬件的都绕不开。我在实际项目里前前后后用了起码几十种变体最大感受是它入门确实不难但想用得稳、选得对背后有非常多门道。这篇内容整理了10个我在真实工作中反复用过的最小二乘拟合法案例从一条直线拟合到带误差棒的加权拟合从高斯峰到分段线性再到圆拟合基本把工程和科研里最常见的场景都覆盖了。适合刚接触数据拟合的同学照着抄也适合已经写过polyfit、但总觉得拟合结果心里没底的同行对照着排查。先把结论放在前面最小二乘不是调用一个函数那么简单模型形式的选择、初值的给法、数据预处理的方式、权重怎么设置每一步都在决定拟合结果靠不靠谱。这篇就按案例一个个说每个都给思路、给代码、给踩坑经验。1. 最小二乘到底在做什么1.1 一句话理解最小二乘最小二乘的核心思想特别朴素找一组模型参数让模型预测值和实际观测值之间的误差平方和最小。误差平方和长这样S Σ (y_i - f(x_i, θ))²其中f是模型函数θ是待求参数y_i是观测值。为什么用平方而不是绝对值两个原因一是平方函数处处可导求极值方便二是平方对大误差的惩罚更重这既能让整体拟合更紧贴数据主体也让离群点对结果的影响被放大。后面这点既是优点也是坑后面案例里面讲。线性最小二乘有解析解也就是直接算出参数不需要迭代。设模型为y Xθ那么解就是θ (XᵀX)⁻¹Xᵀy这就是正规方程。实际再用np.linalg.lstsq或者np.polyfit都能直接得到结果。非线性模型没有解析解就要靠数值迭代比如高斯牛顿法、莱文伯格-马夸特法也就是scipy.optimize.curve_fit背后干的事。1.2 为什么工程里不能只会polyfit很多初学者拿到一组点就np.polyfit(x, y, 2)也不看模型背后的物理意义结果就是拟合曲线画出来挺好看但是参数解释不了换一组数据立刻报废。最小二乘的真正难点其实有两个。第一是模型选择也就是f(x, θ)应该长成什么样这需要你对数据来源有基本判断比如衰减曲线选指数、峰形选高斯或洛伦兹、周期信号选正弦。第二是数值稳定性高次多项式、量级差异很大的数据、相关性很强的参数都会让最小二乘解变得非常脆弱。这篇文章的10个案例我特意按线性模型到非线性模型的梯度排了一遍前面5个大多可以先转成线性再求解后面5个基本得靠curve_fit或专用方法。你把这个顺序走完脑子里就会建立一张什么数据用什么模型的索引表。2. 十个案例上线性、多项式、对数、指数、幂律2.1 案例一直线拟合与传感器标定先来最简单的。传感器标定是最典型的直线拟合场景采集ADC读数x和标准温度y然后想建立y a x b的关系用来把之后的ADC读数转换成温度。这种场景直接上np.polyfit就行但有一件事我会特别提醒先画散点图确认线性关系真的存在再去做拟合。否则一个弯曲的数据带硬用直线拟合斜率截距虽然能算出来但没有意义。import numpy as np import matplotlib.pyplot as plt # 模拟一组传感器标定数据真实温度 20~80°CADC读数含噪声 x_adc np.linspace(200, 800, 20) true_a 0.1 true_b 5.0 y_temp true_a * x_adc true_b np.random.normal(0, 0.8, sizex_adc.shape) # 一次多项式拟合 a, b np.polyfit(x_adc, y_temp, 1) # 拟合优度计算R² y_pred a * x_adc b ss_res np.sum((y_temp - y_pred) ** 2) ss_tot np.sum((y_temp - np.mean(y_temp)) ** 2) r_squared 1 - ss_res / ss_tot print(f拟合斜率 a {a:.4f}, 截距 b {b:.4f}, R² {r_squared:.4f})这里截距b在传感器标定里通常表示零位偏置斜率a表示灵敏度。如果b偏大说明传感器零点有漂移可能要做硬件校准。这种物理解读才是拟合的价值所在而不只是画条线过去。实操心得polyfit默认返回的是降幂排列的系数也就是最高次在前。如果你习惯按a, b这样读很容易搞反。我平时只用一次多项式所以更习惯直接写np.polyfit(x, y, 1)然后自己知道返回值是a, b但到三次以上的时候用np.poly1d来包一层调用起来更不容易弄混。2.2 案例二多项式拟合与应力应变曲线多项式拟合是polyfit最直接的延展但也是最容易用错的地方。材料拉伸实验的应力应变曲线常常在弹性段之后出现非线性段用一个三次多项式y a x³ b x² c x d去近似在工程上很常见。deg选多少是这里最大的问题。deg太低欠拟合曲线形状抓不住。deg太高过拟合拟合曲线会把噪声都吃进去Runge现象一出现曲线两端剧烈震荡完全没法用。# 模拟应力应变数据 strain np.linspace(0, 0.1, 30) stress 2000 * strain**3 - 50 * strain**2 800 * strain np.random.normal(0, 2, sizestrain.shape) # 比较3次和8次多项式拟合 coeff3 np.polyfit(strain, stress, 3) coeff8 np.polyfit(strain, stress, 8) p3 np.poly1d(coeff3) p8 np.poly1d(coeff8) strain_smooth np.linspace(0, 0.1, 300) plt.scatter(strain, stress, labeldata) plt.plot(strain_smooth, p3(strain_smooth), labeldeg3) plt.plot(strain_smooth, p8(strain_smooth), labeldeg8) plt.legend() plt.show()这里的核心经验是多项式拟合的阶数能低就别高一般工程数据不超过5次。判断阶数是否合适不要只看R²要去看残差。如果残差呈现明显的弯曲趋势说明模型还有结构没拟合进去如果残差近乎随机那这个阶数就够了。另外当自变量的数值很大比如10000、20000且次数高时正规方程会病态。解决办法是对x做归一化也就是(x - x_mean) / x_std再拟合或者用np.polynomial.Polynomial.fit这个函数它内部会自动处理domain和window的问题。2.3 案例三对数拟合与浓度标定对数形式y a b ln x在化学分析和生物实验里很常见比如某些光谱测定的吸光度与浓度的关系或者心理物理学里的韦伯-费希纳定律。对数拟合本身不是非线性因为令t ln x之后模型就变成了y a b t完全可以用线性最小二乘来做。我这里建议直接用Scipy的curve_fit也是可以的但对于这种能线性化的模型转成线性回归来解其实更快、更稳定。from scipy.optimize import curve_fit def log_model(x, a, b): return a b * np.log(x) x_data np.array([1, 2, 5, 10, 20, 50, 100, 200]) y_data 1.5 0.8 * np.log(x_data) np.random.normal(0, 0.1, sizex_data.shape) # 方式一转成线性拟合 t np.log(x_data) b, a np.polyfit(t, y_data, 1) # 方式二直接非线性拟合 params, _ curve_fit(log_model, x_data, y_data) print(f线性化: a{a:.3f}, b{b:.3f}) print(fcurve_fit: a{params[0]:.3f}, b{params[1]:.3f})这里要注意的是x必须严格大于0因为ln 0和ln负数都没有定义。如果数据里有x接近0的点就非常麻烦你得考虑是不是要换成ln(x c)这种带偏移的形式或者换模型。2.4 案例四指数衰减拟合指数模型y a e^(bx)出现的频率远比想象中高RC电路的放电曲线、荧光寿命衰减、细菌生长曲线、药物在体内的代谢过程全是这个形式。指数拟合很容易踩一个思维定式的坑很多人一上来就两边取对数转线性ln y ln a bx然后用polyfit。这在数据干净、噪声很小的时候没问题但一旦噪声变大取对数后小值部分的噪声会被严重放大导致拟合结果偏向于数值小的区域整体效果反而差。我的建议是如果噪声不大可以用线性化快速估计初始值如果噪声明显直接上curve_fit做非线性拟合。初值用线性化结果来给这样又稳又准。def exp_model(x, a, b): return a * np.exp(b * x) x_data np.linspace(0, 5, 50) true_a 2.5 true_b -0.6 y_data true_a * np.exp(true_b * x_data) np.random.normal(0, 0.1, sizex_data.shape) # 先用对数线性化选初值 log_y np.log(y_data[y_data 0]) x_pos x_data[y_data 0] b0, log_a0 np.polyfit(x_pos, log_y, 1) a0 np.exp(log_a0) # 再去做非线性拟合 params, cov curve_fit(exp_model, x_data, y_data, p0[a0, b0]) a_fit, b_fit params print(f拟合结果: a{a_fit:.3f}, b{b_fit:.3f})curve_fit返回的cov是参数协方差矩阵对角线开根号就是参数的标准差。我强烈建议你每次都打印一下这个不确定性哪怕只是大概看看量级。如果标准差比参数本身还大那说明数据对参数没有约束力模型可能选错了或者数据范围不足。2.5 案例五幂律拟合幂律形式y a x^b在物理学里到处都是异速生长、湍流标度律、地震频度-震级关系、城市规模分布、网站流量分布。它的拟合方式和指数拟合的逻辑相似但取对数时使用的是log-logln y ln a b ln x也就是说如果数据在双对数坐标下呈直线那就符合幂律关系。def power_model(x, a, b): return a * np.power(x, b) x_data np.logspace(0, 3, 30) true_a 1.2 true_b 1.5 y_data true_a * np.power(x_data, true_b) np.random.normal(0, 3, sizex_data.shape) # log-log线性化 log_x np.log(x_data) log_y np.log(y_data) b0, log_a0 np.polyfit(log_x, log_y, 1) a0 np.exp(log_a0) params, _ curve_fit(power_model, x_data, y_data, p0[a0, b0], maxfev5000) print(fa{params[0]:.3f}, b{params[1]:.3f})幂律拟合在实际数据上有个很隐蔽的坑数据范围太小的话指数b的估计极不稳定。比如你只有x从1到10的数据想判断是不是幂律根本区分不出幂律和指数因为两者的曲线在窄范围内都可以被另一种形式近似。所以做幂律分析至少要覆盖两个数量级的x范围否则结论不可靠。3. 十个案例下高斯、洛伦兹、正弦、加权、分段线性与圆拟合3.1 案例六高斯峰拟合光谱分析、色谱分析、XRD衍射峰分析里峰形拟合是最常见的需求。高斯峰模型长这样y A exp(-(x - μ)² / (2σ²)) B其中A是峰高μ是峰位置σ是峰宽B是基线。这种模型用curve_fit做非线性最小二乘。初值给得好不好决定了非线性拟合是2步收敛还是200步都不收敛。我的固定套路是峰高A用max(y) - min(y)来估计峰位置μ用np.argmax(y)对应的x值σ用半高全宽估计一下基线B用数据左右端点的平均值。def gaussian(x, a, mu, sigma, b): return a * np.exp(-(x - mu)**2 / (2 * sigma**2)) b x_data np.linspace(0, 10, 200) true_mu 5.0 y_data 3.0 * np.exp(-(x_data - true_mu)**2 / (2 * 0.8**2)) 1.0 np.random.normal(0, 0.05, sizex_data.shape) # 初值估计 A0 np.max(y_data) - np.min(y_data) mu0 x_data[np.argmax(y_data)] sigma0 0.5 # 半经验给个值 B0 np.mean([y_data[0], y_data[-1]]) params, cov curve_fit(gaussian, x_data, y_data, p0[A0, mu0, sigma0, B0]) print(f峰高{params[0]:.3f}, 峰位{params[1]:.3f}, σ{params[2]:.3f}, 基线{params[3]:.3f})这里的经验是峰形数据一定要先做基线扣除或者把基线B作为参数一起拟合不然峰高和峰宽都会被基线带偏。另一个很容易犯的错是把高斯拟合和正态分布概率密度拟合当成一回事模型本身是一样但A的含义完全不同写代码的时候别搞混了。3.2 案例七洛伦兹峰拟合洛伦兹线型是另一类常见峰形公式为y A / (1 ((x - x0) / γ)²) B它和高斯的差别在于高斯峰衰减快尾巴薄洛伦兹峰有厚尾巴在机械谐振、原子光谱、Raman光谱里更常见。两种线形放到一起光靠眼睛很难分辨需要用拟合结果去比较残差和不确定度。def lorentzian(x, a, x0, gamma, b): return a / (1 ((x - x0) / gamma)**2) b x_data np.linspace(0, 10, 300) y_data 5.0 / (1 ((x_data - 5.0) / 0.5)**2) 0.5 np.random.normal(0, 0.05, sizex_data.shape) A0 np.max(y_data) - np.min(y_data) x00 x_data[np.argmax(y_data)] gamma0 0.5 B0 np.mean([y_data[0], y_data[-1]]) params, _ curve_fit(lorentzian, x_data, y_data, p0[A0, x00, gamma0, B0]) print(f峰高{params[0]:.3f}, 峰位{params[1]:.3f}, γ{params[2]:.3f})实际工作中同一样品可能会用高斯和洛伦兹各拟合一次看哪个残差更小。这里有个经验如果两种模型的R²差别很小优先选物理上更合理的那个不要被纯数值结果带着走。3.3 案例八正弦信号拟合正弦拟合是周期性信号处理里的常见需求模型是y A sin(2π f t φ) C这里的频率f是整个模型最敏感的参数初值稍微偏一点拟合就很容易掉进局部最优。我处理这种问题的习惯是先用FFT粗估频率再把它作为初值传给curve_fit这样成功率能提升非常多。def sin_model(t, A, f, phi, C): return A * np.sin(2 * np.pi * f * t phi) C # 模拟信号 t np.linspace(0, 1, 500) f_true 50.0 y_data 2.0 * np.sin(2 * np.pi * f_true * t 0.3) 1.0 np.random.normal(0, 0.1, sizet.shape) # FFT粗估频率 dt t[1] - t[0] fft_vals np.fft.rfft(y_data - np.mean(y_data)) freqs np.fft.rfftfreq(len(t), dt) f0 freqs[np.argmax(np.abs(fft_vals))] # 拟合 A0 np.ptp(y_data) / 2 C0 np.mean(y_data) phi0 0.0 params, _ curve_fit(sin_model, t, y_data, p0[A0, f0, phi0, C0], maxfev10000) print(f频率{params[1]:.3f} Hz)这里要注意FFT频域分辨率是1/TT是信号总时长。如果信号只有几个周期FFT初值的误差可能较大但哪怕粗略估计也比随机给一个初值强得多。曲线拟合在频率接近真实值之后收敛很快通常几十次迭代就能搞定。3.4 案例九带误差棒的加权最小二乘前面所有案例都默认每个点的权重一样。但实验数据通常每个点有不同误差比如同一位置重复测量了多次标准差不同。这时候就得上加权最小二乘权重取1/σ²让误差小的点主导拟合误差大的点少参与。在curve_fit里直接用sigma参数传入每个点的标准差即可默认是None即等权重。x_data np.linspace(1, 10, 10) y_true 2.0 1.5 * x_data # 不同点给不同的噪声水平 sigma np.array([0.1, 0.15, 0.2, 0.3, 0.5, 0.8, 1.0, 1.2, 1.5, 2.0]) y_data y_true np.random.normal(0, sigma) # 加权拟合sigma传入标准差 params_w, _ curve_fit(lambda x, a, b: a b * x, x_data, y_data, sigmasigma, absolute_sigmaTrue) # 不加权拟合 params_nw, _ curve_fit(lambda x, a, b: a b * x, x_data, y_data) print(f加权拟合: a{params_w[0]:.3f}, b{params_w[1]:.3f}) print(f不加权拟合: a{params_nw[0]:.3f}, b{params_nw[1]:.3f})absolute_sigma这个参数很重要。如果sigma只代表相对权重不关心绝对量级设False如果sigma是真实的标准差估计设True协方差矩阵才能给出正确的参数误差。我见过很多人在这上面栽跟头拟合参数没问题但误差低估了一两个量级。手动实现加权最小二乘也不复杂把每个点都除以σ变成y_i/σ_i (1/σ_i) * f(x_i, θ)再用普通最小二乘解。本质上加权就是数据变换。3.5 案例十分段线性拟合与圆拟合最后两个在工业里很常用我合在一个小节讲。分段线性拟合适合那种在不同区间呈现不同线性趋势的数据比如材料在弹塑性转变前后的曲线、物体在启动和匀速阶段的位置-时间曲线。断点位置可以手工指定也可以作为参数参与拟合。手工指定简单但不够灵活。用curve_fit把断点也当参数问题是目标函数在断点处不可导优化器容易卡住或发散。一个更稳的方案是用pwlf库它专门处理分段线性拟合断点数量可配断点位置自动优化# 需要先安装pip install pwlf import pwlf x_data np.linspace(0, 10, 100) y_data np.piecewise(x_data, [x_data 4, x_data 4], [lambda x: 1.0 * x 0.5, lambda x: -0.8 * x 7.7]) y_data np.random.normal(0, 0.2, sizex_data.shape) my_pwlf pwlf.PiecewiseLinFit(x_data, y_data) breaks my_pwlf.fit(2) # 2段1个断点 print(断点位置:, breaks)圆拟合是机器视觉里定位圆孔、检测工件圆度的经典需求。最小二乘拟合圆的核心思路是把圆方程x² y² Dx Ey F 0变成线性最小二乘问题解出D、E、F之后圆心为(-D/2, -E/2)半径为√((D²E²)/4 - F)。# 模拟圆弧上的点 theta np.linspace(0, np.pi, 50) cx_true, cy_true, r_true 3.0, 4.0, 2.0 x_data cx_true r_true * np.cos(theta) np.random.normal(0, 0.05, sizetheta.shape) y_data cy_true r_true * np.sin(theta) np.random.normal(0, 0.05, sizetheta.shape) # 代数拟合 A_mat np.vstack([x_data**2 y_data**2, x_data, y_data, np.ones_like(x_data)]).T coeffs, _, _, _ np.linalg.lstsq(A_mat, -np.ones_like(x_data), rcondNone) D, E, F coeffs[1], coeffs[2], coeffs[3] cx -D / 2 cy -E / 2 r np.sqrt((D**2 E**2) / 4 - F) print(f圆心: ({cx:.3f}, {cy:.3f}), 半径: {r:.3f})这里有个值得关注的细节代数圆拟合在圆弧覆盖角度很小时比如不到90度半径估计偏差极大会系统性偏小。如果数据只能覆盖一小段圆弧建议改用几何距离迭代方法比如Taubin算法或Pratt算法。实际视觉项目中圆孔通常能看到大半个圆代数拟合够用但如果只能看到一小段弧形轮廓一定要额外评估圆参数的稳定性。4. 实操中最容易翻车的四个环节4.1 模型选错了拟合再漂亮也没用很多拟合问题本质不是数值问题而是模型识别问题。一组衰减数据你用指数模型还是幂律模型在窄数据范围里拟合出来的曲线可能几乎重合但外推结果天差地别。我的习惯是先看物理背景从机制上去推断模型然后在双对数、半对数坐标下画图看数据在什么坐标下更接近直线最后再用拟合残差验证。举一个真实的例子某个热敏电阻的阻值-温度关系数据有人用三次多项式拟合在数据范围内误差很小但外推到低温区间完全不合理。改用Steinhart-Hart方程本身就是基于物理机制的对数模型之后虽然数据范围内拟合误差稍大但外推行为合理得多。提示拟合的目的不是把现有数据点全部踩在脚下而是让模型能够解释数据背后真实的生成机制。这一点想不清楚后面全是白忙。4.2 初值不给curve_fit原地发脾气curve_fit的p0如果不给默认全部初值是1。这对于很多模型来说完全不在收敛域内结果就是要么迭代超时要么收敛到一个局部最优解曲线歪得离谱。我自己的固定工作流是三步第一步用图形或者简单统计粗估参数。比如高斯峰峰高就是max-min峰位就是argmax对应的x指数衰减就先取对数线性拟合得到a、b。第二步把粗估值作为p0传给curve_fit。第三步把拟合结果再作为初值跑一次看两次的答案是否一致。如果两次结果有明显出入说明参数空间里存在多个局部极小模型本身对参数分辨力不足需要加约束或者换模型。4.3 高次多项式怎么防病态高次多项式拟合最大的敌人是矩阵病态。当x值很大比如从几千到几万再来个6次多项式正规方程(XᵀX)⁻¹这一步几乎必然引起严重的数值误差。解决办法有三个方向。第一对x做标准化比如(x - mean(x)) / std(x)拟合完再还原系数。第二使用np.polynomial.Polynomial.fit它会自动处理domain和window。第三更高阶的情况直接用正交多项式比如Legendre多项式做基函数从根本上避免病态。在普通工程场景下第一种和第二种足够用。我强烈不建议硬上8次、10次多项式除非有明确的物理依据否则高阶多项式本身就是过拟合的代名词。4.4 残差和不确定度的理解R²是大多数人看拟合好坏的唯一指标但它有严重盲区。R²高不代表模型正确只代表模型在数据范围里拟合出了一部分方差。即使模型形式完全错误R²也可能高达0.95以上。正确的做法是画残差图横轴是x纵轴是y - y_pred。残差如果随机散布在零线附近说明模型捕捉到了所有结构如果残差还有明显的曲线趋势说明还有信号没有被模型解释。至于参数不确定度记住一句话拟合得到的参数不是点值是有分布的点估计。协方差矩阵对角线开根号给出标准误在写报告时最好带上比如a 2.35 ± 0.12。这样别人能直接判断这个参数的可信度。5. 工具选型什么场景用什么函数5.1 常用函数对比表在Python生态里最小二乘相关工具看似很多实际上每个函数都有明确的定位和边界。用错了工具不一定会出错但一定会写出很多多余代码或者在某些场景下踩到性能陷阱。工具适用场景注意点numpy.polyfit一次或低次多项式拟合系数顺序是高次在前高次需小心病态numpy.linalg.lstsq自定义线性模型多元回归万能线性最小二乘返回残余平方和scipy.optimize.curve_fit非线性模型通用拟合需要给初值支持sigma权重scipy.optimize.least_squares需要加边界约束、鲁棒损失比curve_fit更底层灵活性高statsmodels.OLS需要完整统计推断、p值置信区间适合论文和正式报告pwlf.PiecewiseLinFit分段线性拟合断点位置自动优化sklearn.linear_model.LinearRegression大数据量、多特征适合机器学习场景有predict方法这里要特别强调一下curve_fit和least_squares的关系curve_fit是least_squares的封装它内部已经帮你处理了雅可比矩阵、参数上下界、权重这些最常见的需求。只有当你要自定义损失函数、启用鲁棒核函数比如Huber损失时才需要直接上least_squares。5.2 我自己的选择习惯我在实际项目里通常遵循一条简单的原则能线性化就线性化不能用线性化就curve_fit再不行才上least_squares带约束。为什么我偏好线性化解法两个原因。第一线性最小二乘有唯一解没有初值问题不会因为初值没给对就发散。第二速度快。非线性拟合每次都在迭代几百个数据点还好几万个数据点的话耗时就上来了。对于只需要快速看一下趋势的脚本我直接用np.polyfit。对于需要严谨参数误差分析的正式实验我用curve_fit加absolute_sigmaTrue把协方差矩阵一起打出来。对于有物理边界约束的场景比如参数必须为正、频率必须在一定范围我在curve_fit里传bounds。这里有个很少人提的细节bounds如果传了curve_fit会自动换成trf算法收敛行为会变化所以设置了边界之后一定要再看一眼拟合曲线和残差防止收敛到边界附近的伪最优。注意很多时候拟合结果不理想不是工具的问题而是数据本身的问题。数据的采集范围、噪声水平、异常点这些因素决定了拟合结果的极限。工具只负责在给定数据的前提下取最优它无法变出数据里本就不存在的信息。最后分享一个我改不掉的习惯每做完一次拟合我一定把原始数据点和拟合曲线画在同一张图里并把残差图放在下方。这一步看起来多余却能最快暴露问题。哪怕数值指标全是绿的只要图形上发现拟合曲线在某一段系统性地偏离数据我都会停下来重新检查模型而不是急着把结果写进报告。因为拟合最大的魅力恰恰藏在那些看起来差不多但总觉得哪里不对劲的细节里。本文还有配套的精品资源点击获取