三次样条插值与多项式拟合：从数学原理到MATLAB/Python实战

1. 从“画线”到“造路”：为什么我们需要插值与拟合？

做数据分析或者工程计算的朋友，经常会遇到一个头疼的问题：手里只有一堆离散的数据点，像夜空里的星星，稀稀拉拉。但我们需要知道星星之间那片“黑暗”里，到底藏着什么。比如，你每隔一小时测一次室外温度，得到了12个数据点，但你想知道凌晨3点15分到底有多冷；又比如，你通过实验得到了一组材料在不同应力下的形变数据，但需要预测一个实验没做过的应力值对应的形变。

这时候，你就需要一种方法，根据已知的“星星”（数据点），去构建一条穿过或者贴近它们的、连续的“路径”（函数）。这条路径，就是我们理解数据内在规律、进行预测和分析的桥梁。而“插值”和“拟合”，就是搭建这座桥梁最核心的两种数学工具。它们都源于同一个需求，但思路和适用场景截然不同，选错了方法，就像在沼泽地里开跑车，再好的算法也跑不起来。

简单来说：

• 插值：追求“精确穿过”。构造的函数曲线必须严丝合缝地经过每一个已知数据点。它假设你的数据点本身是精确无误的，目标是在点与点之间进行“内插”，填补空白。就像用尺子画直线连接各个点，或者用更光滑的曲线把它们串起来。三次样条插值就是插值家族中的“明星成员”，因为它平衡了光滑性和计算复杂度。
• 拟合：追求“大势所趋”。承认数据可能存在误差（比如测量噪声），不要求曲线穿过每一个点，而是寻找一条从整体趋势上看最贴近所有数据点的曲线。它更关注数据背后隐藏的宏观规律。多项式拟合是最基础、最直观的拟合方法，试图用一个多项式函数来概括数据的走势。

那么，什么时候用插值，什么时候用拟合呢？我个人的经验是：如果你的数据点本身精度很高，近乎“真理”，并且你需要得到点与点之间任意位置的值，那就用插值。如果你的数据带有明显的噪声或误差，你更关心整体趋势和预测，那就用拟合。举个例子，你有一张精确的等高线地图上几个点的海拔（插值适用），和你有一群人在不同年龄的身高数据，但每个人的测量都有细微误差（拟合适用）。

接下来的内容，我将结合近十年在科学计算和工程分析中的实战经验，以MATLAB和Python这两个最主流的工具为例，手把手带你深入三次样条插值和多项式拟合的“内核”。我们不止讲怎么调用函数，更要讲清楚函数背后的数学逻辑、参数选择的门道，以及那些官方文档里不会写、但实践中一踩一个准的“坑”。

2. 光滑的“弹簧尺”：三次样条插值原理深度拆解

想象一下，你有一把特制的“弹簧尺”，它被固定在几个已知的点上（我们的数据点）。这把尺子非常柔软，但在固定点处又非常听话。我们希望这把尺子在连接这些点时，整体看起来非常光滑流畅，没有突兀的折角。这就是三次样条插值要解决的问题。

为什么是“三次”？为什么叫“样条”？这得从它的数学本质说起。早期的工程师或绘图员，会用一根有弹性的细木条或金属条（称为“样条”）固定在几个点上，让它自然弯曲形成一条光滑的曲线。数学上，这条曲线可以用分段的三次多项式来完美模拟。

2.1 核心思想：分段三次，全局光滑

三次样条插值的基本策略是：

1. 分段处理：将整个数据区间，按照已知数据点的x坐标，分割成若干个小区间。在每个小区间[x_i, x_{i+1}]上，用一个独立的三次多项式S_i(x)来连接左右两个端点。S_i(x) = a_i + b_i(x - x_i) + c_i(x - x_i)^2 + d_i(x - x_i)^3这里有四个未知系数(a_i, b_i, c_i, d_i)。
2. 连接条件：为了保证整条曲线看起来是“一根”而不是“几段拼凑的”，我们在每个内节点（非边界的数据点）x_i处，要求左右两段多项式满足三个“光滑”条件：
   • 函数值连续：左边多项式在x_i的值等于右边多项式在x_i的值。这保证了曲线不断开。S_{i-1}(x_i) = S_i(x_i) = y_i
   • 一阶导数连续：左边多项式在x_i的斜率（一阶导数）等于右边多项式在x_i的斜率。这保证了曲线没有尖角，是光滑的。S'_{i-1}(x_i) = S'_i(x_i)
   • 二阶导数连续：左边多项式在x_i的弯曲程度（二阶导数）等于右边多项式在x_i的弯曲程度。这保证了曲线的弯曲变化也是平滑的，没有突兀的转折。S''_{i-1}(x_i) = S''_i(x_i)
3. 边界条件：上面这些条件只确定了内部节点的关系，整条曲线的“起势”和“收势”还需要额外规定，这就是边界条件。最常见的两种是：
   • 自然样条 (Natural Spline)：规定曲线在两个端点处的二阶导数为零，即S''(x_0) = 0和S''(x_n) = 0。这好比让那根“弹簧尺”在两端处于自然放松、不受弯矩的状态。这样得到的曲线在端点附近可能会有些“平直”。
   • 固定斜率样条 (Clamped Spline)：直接指定曲线在两个端点处的一阶导数值，即S'(x_0) = f'_0和S'(x_n) = f'_n。如果你知道数据在边界的变化趋势（比如物理上的速度、梯度），用这个条件会更准确。
   • 非扭结样条 (Not-a-Knot)：这是一种更“聪明”的默认条件，它要求第一个和第二个内部节点处的三阶导数也连续。相当于去掉了第一个和最后一个内节点，让曲线在整体上更“自然”，通常是软件里的默认选项。

通过以上条件，我们可以列出一个关于所有分段多项式系数的线性方程组。求解这个方程组，就得到了唯一确定的三次样条函数。这个过程虽然复杂，但好在MATLAB和Python都为我们封装好了。

2.2 MATLAB实战：spline与csape的玄机

在MATLAB中，最常用的函数是spline和interp1（选择'spline'方法）。

基础用法：

% 假设我们有原始数据点 x = [0, 1, 2, 3, 4]; y = [0, 0.5, 2, 1.5, 0]; % 方法1：使用 spline 函数，它默认使用非扭结(not-a-knot)边界条件 pp = spline(x, y); % 返回一个结构体，包含了分段多项式的所有系数信息 % 计算在更细网格上的插值结果 xx = linspace(0, 4, 100); yy_spline = ppval(pp, xx); % 使用ppval计算插值 % 方法2：使用 interp1 函数，更直观 yy_interp = interp1(x, y, xx, 'spline'); % 绘图对比 figure; plot(x, y, 'ro', 'MarkerSize', 10, 'LineWidth', 2); hold on; plot(xx, yy_spline, 'b-', 'LineWidth', 1.5); plot(xx, yy_interp, 'g--', 'LineWidth', 1.5); legend('原始数据点', 'spline函数', 'interp1(spline)'); title('MATLAB 三次样条插值对比'); grid on;

你会发现yy_spline和yy_interp的结果几乎完全一致。interp1内部调用的就是spline算法。

进阶控制：csape函数如果你需要对边界条件进行精细控制，曲线拟合工具箱里的csape函数是更强大的选择。

% 使用 csape，需要安装曲线拟合工具箱 (Curve Fitting Toolbox) % 指定边界条件：第二类边界条件（自然样条，二阶导为零） pp_natural = csape(x, y, 'second'); % 'second' 表示指定二阶导，默认值为[0,0]即自然样条 % 指定边界条件：第一类边界条件（固定斜率） % 假设我们知道左端点斜率为1，右端点斜率为-1 clamped_conds = [1, -1]; % 第一个元素是左端点一阶导，第二个是右端点一阶导 pp_clamped = csape(x, y, 'complete', clamped_conds); % 计算并绘图 yy_natural = ppval(pp_natural, xx); yy_clamped = ppval(pp_clamped, xx); figure; plot(x, y, 'ko', 'MarkerSize', 10, 'LineWidth', 2); hold on; plot(xx, yy_spline, 'b-', 'DisplayName', 'spline (not-a-knot)'); plot(xx, yy_natural, 'r--', 'LineWidth', 2, 'DisplayName', 'csape (natural)'); plot(xx, yy_clamped, 'g:', 'LineWidth', 2, 'DisplayName', 'csape (clamped, slope=[1,-1])'); legend('show'); title('不同边界条件的三次样条插值对比'); grid on;

运行这段代码，你可以清晰地看到不同边界条件对曲线两端形态的影响。自然样条在端点处趋于平直，而固定斜率样条则严格按照你指定的方向“出发”和“到达”。

实操心得：对于大多数内插问题，MATLAB默认的spline（非扭结条件）效果已经非常好，且计算稳定。只有在你有明确的物理约束（如已知边界梯度）或数据在边界有特殊行为时，才需要考虑使用csape指定边界条件。盲目使用自然样条有时会在边界引入不希望的“平坦化”效应。

2.3 Python实战：SciPy库的CubicSpline与interp1d

在Python的科学计算栈中，SciPy库是处理这类问题的绝对主力。

基础用法：scipy.interpolate.CubicSpline这是最现代、接口最清晰的三次样条插值类。

import numpy as np import matplotlib.pyplot as plt from scipy.interpolate import CubicSpline # 原始数据点 x = np.array([0, 1, 2, 3, 4]) y = np.array([0, 0.5, 2, 1.5, 0]) # 创建CubicSpline对象 # bc_type参数指定边界条件： # 'not-a-knot' (默认), 'natural', 'clamped'等 cs_notaknot = CubicSpline(x, y, bc_type='not-a-knot') # 默认 cs_natural = CubicSpline(x, y, bc_type='natural') # 自然样条 # 对于‘clamped’，需要提供边界的一阶导数值 # 假设左端点导数为1，右端点为-1 cs_clamped = CubicSpline(x, y, bc_type=((1, -1), )) # 生成插值点 xx = np.linspace(0, 4, 100) yy_notaknot = cs_notaknot(xx) yy_natural = cs_natural(xx) yy_clamped = cs_clamped(xx) # 绘图 plt.figure(figsize=(10, 6)) plt.plot(x, y, 'ko', label='Data points', markersize=10) plt.plot(xx, yy_notaknot, 'b-', label='CubicSpline (not-a-knot)', linewidth=2) plt.plot(xx, yy_natural, 'r--', label='CubicSpline (natural)', linewidth=2) plt.plot(xx, yy_clamped, 'g:', label='CubicSpline (clamped, [1,-1])', linewidth=2) plt.legend() plt.grid(True, linestyle='--', alpha=0.7) plt.title('Python SciPy CubicSpline with Different Boundary Conditions') plt.xlabel('x') plt.ylabel('y') plt.show()

另一种选择：scipy.interpolate.interp1d这是一个更通用的插值函数工厂，可以通过kind参数选择多种插值方法，包括三次样条。

from scipy.interpolate import interp1d # 使用 interp1d 进行三次样条插值 # kind='cubic' 在scipy中通常指三次样条 f_cubic = interp1d(x, y, kind='cubic') # 注意：旧版本可能默认不是not-a-knot yy_interp1d = f_cubic(xx) # 对比 plt.figure(figsize=(10, 6)) plt.plot(x, y, 'ko', markersize=10, label='Data points') plt.plot(xx, cs_notaknot(xx), 'b-', label='CubicSpline (not-a-knot)', linewidth=2) plt.plot(xx, yy_interp1d, 'm-.', label='interp1d (kind="cubic")', linewidth=2) plt.legend() plt.grid(True, linestyle='--', alpha=0.7) plt.title('CubicSpline vs interp1d(kind="cubic")') plt.show()

注意：interp1d的kind='cubic'在较新的SciPy版本中指的是三次样条，但其内部实现和默认边界条件可能与CubicSpline略有不同，对于关键应用，建议使用CubicSpline因为它对边界条件的控制更明确。

高阶信息获取CubicSpline对象的一个巨大优势是，你可以轻松获取插值函数在各点的导数，这对于很多工程分析（如求速度、加速度）非常有用。

# 计算一阶导数（斜率）和二阶导数（曲率） dy_dx = cs_notaknot(xx, 1) # 参数1表示一阶导 d2y_dx2 = cs_notaknot(xx, 2) # 参数2表示二阶导 fig, axes = plt.subplots(1, 3, figsize=(15, 4)) axes[0].plot(xx, cs_notaknot(xx)) axes[0].set_title('Interpolation') axes[0].grid(True) axes[1].plot(xx, dy_dx) axes[1].set_title('First Derivative') axes[1].grid(True) axes[2].plot(xx, d2y_dx2) axes[2].set_title('Second Derivative') axes[2].grid(True) plt.tight_layout() plt.show()

3. 当插值“过犹不及”：龙格现象与样条的稳定性

在讨论多项式拟合之前，必须提一下高次多项式插值的一个著名陷阱——龙格现象 (Runge‘s Phenomenon)。这能帮你理解为什么我们更偏爱分段低次（如三次）的样条插值，而不是一个贯穿所有点的单一高次多项式。

考虑在区间[-1, 1]上，用等距节点对函数f(x) = 1 / (1 + 25x^2)进行插值。随着插值多项式次数的增加（即使用更多的点），在区间两端会出现剧烈的振荡，误差不但没有减小，反而爆炸式增长。

import numpy as np import matplotlib.pyplot as plt def runge(x): return 1 / (1 + 25 * x**2) x_true = np.linspace(-1, 1, 400) y_true = runge(x_true) plt.figure(figsize=(12, 8)) plt.plot(x_true, y_true, 'k-', label='True function: $1/(1+25x^2)$', linewidth=3) for n in [5, 10, 15]: # 使用5, 10, 15个等距节点 x_nodes = np.linspace(-1, 1, n) y_nodes = runge(x_nodes) # 使用numpy的polyfit进行多项式插值（拟合阶数为n-1） coeff = np.polyfit(x_nodes, y_nodes, n-1) # 注意：高次多项式拟合数值不稳定，这里仅为演示 p = np.poly1d(coeff) y_poly = p(x_true) plt.plot(x_nodes, y_nodes, 'o', markersize=8) plt.plot(x_true, y_poly, '--', linewidth=2, label=f'Poly Interp (n={n})') plt.legend(loc='best') plt.ylim(-2, 2) plt.grid(True, alpha=0.3) plt.title("Runge's Phenomenon: High-degree Polynomial Interpolation Fails") plt.show()

运行这段代码，你会看到当n=15时，多项式在区间两端已经疯狂振荡，完全偏离了真实函数。这就是龙格现象：对于某些函数，在等距节点上使用高次多项式插值，在区间边缘会产生巨大的误差。

为什么三次样条能避免这个问题？

1. 分段低次：样条插值将整个区间分成很多小段，每段只用三次多项式。三次多项式的行为非常温和，不容易产生剧烈振荡。
2. 局部性：修改一个数据点或一段区间，主要只影响相邻的几段，而不会像高次多项式那样“牵一发而动全身”，影响全局。
3. 光滑性：通过强制一阶、二阶导数连续，保证了拼接处的平滑过渡，避免了生硬的拐角。

核心教训：当你有很多数据点时（比如超过10个），绝对不要试图去构造一个穿过所有点的单一高次多项式来做插值。三次样条插值几乎是这种情况下默认的、安全的选择。它用牺牲一点点全局解析表达式的简洁性，换来了极高的数值稳定性和局部可控性。

4. 抓住主要矛盾：多项式拟合的本质与最小二乘法

现在，让我们把目光从“精确穿过”转向“大势所趋”。拟合承认数据有瑕疵，我们的目标是找到一条曲线，它能以最小的“总体偏差”来代表这组数据。最常用的准则就是最小二乘法 (Least Squares)：让所有数据点到拟合曲线的垂直距离的平方和最小。

为什么是平方和？而不是绝对值和？简单来说，平方项对大的误差惩罚更重（让拟合曲线更“厌恶”离群点），并且在数学上可导，能导出漂亮的解析解（线性方程组）。

4.1 数学原理：从目标函数到正规方程

假设我们想用一个m次多项式来拟合n个数据点(x_i, y_i)：P(x) = a_0 + a_1 x + a_2 x^2 + ... + a_m x^m

我们的目标是找到一组系数[a_0, a_1, ..., a_m]，使得损失函数L最小：L = Σ_{i=1}^{n} [y_i - P(x_i)]^2 = Σ_{i=1}^{n} [y_i - (a_0 + a_1 x_i + ... + a_m x_i^m)]^2

这是一个关于系数a_k的多元二次函数。为了求最小值，我们分别对每个系数a_k求偏导数，并令其等于0：∂L/∂a_k = 0, for k = 0, 1, ..., m

经过推导（具体过程涉及线性代数，这里略过），我们可以得到一个(m+1)阶的线性方程组，称为正规方程 (Normal Equations)：(X^T X) A = X^T Y其中：

• X是一个n x (m+1)的矩阵，其第i行是[1, x_i, x_i^2, ..., x_i^m]（称为设计矩阵）。
• Y是一个n x 1的列向量，包含所有的y_i。
• A是一个(m+1) x 1的列向量，即我们要求的系数[a_0, a_1, ..., a_m]^T。

求解这个线性方程组，就得到了最优的多项式系数。理论上，只要X^T X可逆（通常要求n > m且数据点不特殊），就有唯一解。

4.2 关键抉择：多项式次数m到底选几？

这是多项式拟合中最核心、也最容易出错的问题。次数太低，模型太简单，无法捕捉数据的复杂趋势，称为“欠拟合”。次数太高，模型会拼命去“贴合”每一个数据点，包括噪声点，导致曲线剧烈波动，预测新数据的能力变差，称为“过拟合”。

如何选择？没有银弹，但有以下实战策略：

1. 可视化观察：先把数据点画出来，肉眼观察大致趋势是直线、抛物线还是更复杂的曲线。这是最直观的第一步。
2. 交叉验证：将数据随机分成“训练集”和“验证集”。用训练集拟合不同次数的多项式，然后在验证集上计算误差（如均方误差MSE）。选择在验证集上误差最小的那个次数。这是更可靠的方法。
3. 观察误差下降：绘制“多项式次数”与“训练集误差”的关系图。误差会随着次数增加而持续下降（因为模型越来越复杂）。我们寻找那个“拐点”，即再增加次数，误差下降不再明显的点。
4. 奥卡姆剃刀原则：在效果相近的情况下，永远选择更简单的模型（次数更低）。简单模型的泛化能力通常更强。

一个经典误区：很多人认为拟合的R^2（决定系数）越接近1越好，所以会不断增加次数直到R^2接近1。这是完全错误的！R^2在训练集上会随着模型复杂度（次数）增加而单调增加，即使是对纯噪声数据进行高次拟合，R^2也能接近1。高R^2只说明模型“记住”了训练数据，不代表它“理解”了规律。必须依赖验证集误差或类似准则来判断。

4.3 MATLAB实战：polyfit与polyval的黄金组合

MATLAB中进行多项式拟合异常简单，核心就是polyfit函数。

% 生成带噪声的示例数据 x = linspace(0, 4*pi, 50); y_true = sin(x); % 真实函数 rng(42); % 设置随机种子，保证结果可复现 noise = 0.1 * randn(size(x)); % 加入高斯噪声 y_noisy = y_true + noise; % 尝试用不同次数的多项式拟合 degrees = [1, 3, 5, 10]; % 线性，三次，五次，十次 figure; plot(x, y_true, 'k-', 'LineWidth', 3, 'DisplayName', 'True: sin(x)'); hold on; plot(x, y_noisy, 'bo', 'MarkerSize', 6, 'DisplayName', 'Noisy Data'); for i = 1:length(degrees) m = degrees(i); % polyfit: 返回从高次到低次排列的系数 p_coeffs = polyfit(x, y_noisy, m); % polyval: 利用系数计算多项式值 y_fit = polyval(p_coeffs, x); % 计算拟合优度 R-squared (仅作参考，非决定性指标) y_mean = mean(y_noisy); ss_tot = sum((y_noisy - y_mean).^2); ss_res = sum((y_noisy - y_fit).^2); r2 = 1 - (ss_res / ss_tot); % 绘图 plot(x, y_fit, 'LineWidth', 1.5, 'DisplayName', sprintf('Fit (deg=%d, R^2=%.4f)', m, r2)); end hold off; legend('Location', 'best'); xlabel('x'); ylabel('y'); title('Polynomial Fitting with Different Degrees (MATLAB)'); grid on;

运行后观察：

• 1次（线性）：明显欠拟合，无法捕捉正弦波的波动。
• 3次/5次：开始能捕捉到主要的波动趋势，与真实函数形状接近。
• 10次：曲线开始出现不自然的扭动，特别是在数据稀疏或两端区域，它试图去穿过每一个噪声点，这就是过拟合的苗头。虽然它在训练数据点上的R^2可能最高。

评估过拟合：引入新数据为了更清楚地看到过拟合，我们可以看看拟合出的多项式在新数据（测试集）上的表现。

% 生成新的测试数据（同分布，但未参与拟合） x_test = linspace(0, 4*pi, 200); y_test_true = sin(x_test); y_test_noisy = y_test_true + 0.1 * randn(size(x_test)); % 计算不同次数模型在测试集上的均方误差 (MSE) mse_train = zeros(size(degrees)); mse_test = zeros(size(degrees)); for i = 1:length(degrees) m = degrees(i); p_coeffs = polyfit(x, y_noisy, m); % 用原始训练数据拟合 y_fit_train = polyval(p_coeffs, x); y_fit_test = polyval(p_coeffs, x_test); mse_train(i) = mean((y_noisy - y_fit_train).^2); mse_test(i) = mean((y_test_noisy - y_fit_test).^2); end figure; plot(degrees, mse_train, 's-', 'LineWidth', 2, 'MarkerSize', 10, 'DisplayName', 'Training MSE'); hold on; plot(degrees, mse_test, 'd-', 'LineWidth', 2, 'MarkerSize', 10, 'DisplayName', 'Test MSE'); xlabel('Polynomial Degree'); ylabel('Mean Squared Error (MSE)'); title('Overfitting Visualization: Training vs Test Error'); legend; grid on;

你会观察到，随着次数增加，训练误差持续下降（模型越来越复杂，记住数据的能力越强），但测试误差会先下降后上升。测试误差最低点对应的次数，通常就是泛化能力最好的模型复杂度。这就是为什么不能只看训练集表现。

4.4 Python实战：NumPy的polyfit与sklearn的流程

Python中同样有numpy.polyfit，其用法与MATLAB极其相似。

import numpy as np import matplotlib.pyplot as plt # 生成数据 np.random.seed(42) x = np.linspace(0, 4*np.pi, 50) y_true = np.sin(x) y_noisy = y_true + 0.1 * np.random.randn(len(x)) degrees = [1, 3, 5, 10] plt.figure(figsize=(10, 6)) plt.plot(x, y_true, 'k-', label='True: sin(x)', linewidth=3) plt.plot(x, y_noisy, 'bo', label='Noisy Data', alpha=0.6) for deg in degrees: # numpy.polyfit: 返回从高次到低次排列的系数 coeffs = np.polyfit(x, y_noisy, deg) # numpy.poly1d: 将系数数组转换为多项式函数对象 p = np.poly1d(coeffs) y_fit = p(x) # 计算 R-squared residuals = y_noisy - y_fit ss_res = np.sum(residuals**2) ss_tot = np.sum((y_noisy - np.mean(y_noisy))**2) r_squared = 1 - (ss_res / ss_tot) plt.plot(x, y_fit, linewidth=1.5, label=f'Fit (deg={deg}, $R^2$={r_squared:.4f})') plt.legend(loc='best') plt.grid(True, alpha=0.3) plt.title('Polynomial Fitting with NumPy (Different Degrees)') plt.xlabel('x') plt.ylabel('y') plt.show()

使用sklearn的管道化操作对于更复杂的机器学习流程（如包含标准化、交叉验证），scikit-learn提供了更模块化的接口。虽然多项式拟合本身简单，但用sklearn可以方便地集成到更广的流程中。

from sklearn.preprocessing import PolynomialFeatures from sklearn.linear_model import LinearRegression from sklearn.pipeline import make_pipeline from sklearn.metrics import mean_squared_error, r2_score # 划分训练集和测试集 from sklearn.model_selection import train_test_split x_reshaped = x.reshape(-1, 1) # sklearn要求特征矩阵是二维的 X_train, X_test, y_train, y_test = train_test_split(x_reshaped, y_noisy, test_size=0.3, random_state=42) degrees = [1, 3, 5, 10] train_errors = [] test_errors = [] plt.figure(figsize=(12, 8)) for i, deg in enumerate(degrees): # 创建管道：先构造多项式特征，再进行线性回归 model = make_pipeline(PolynomialFeatures(degree=deg), LinearRegression()) model.fit(X_train, y_train) y_train_pred = model.predict(X_train) y_test_pred = model.predict(X_test) mse_train = mean_squared_error(y_train, y_train_pred) mse_test = mean_squared_error(y_test, y_test_pred) train_errors.append(mse_train) test_errors.append(mse_test) # 为了绘图平滑，在密集点上预测 x_plot = np.linspace(0, 4*np.pi, 300).reshape(-1, 1) y_plot_pred = model.predict(x_plot) plt.subplot(2, 2, i+1) plt.scatter(X_train, y_train, color='blue', s=20, alpha=0.6, label='Train data') plt.scatter(X_test, y_test, color='red', s=20, alpha=0.6, label='Test data') plt.plot(x_plot, y_plot_pred, color='black', linewidth=2, label=f'Deg {deg}') plt.title(f'Degree {deg}\nTrain MSE: {mse_train:.4f}, Test MSE: {mse_test:.4f}') plt.legend(loc='upper right') plt.grid(True, alpha=0.3) plt.tight_layout() plt.show() # 绘制误差随次数变化的曲线 plt.figure(figsize=(8, 5)) plt.plot(degrees, train_errors, 's-', label='Training MSE', markersize=8) plt.plot(degrees, test_errors, 'd-', label='Test MSE', markersize=8) plt.xlabel('Polynomial Degree') plt.ylabel('Mean Squared Error') plt.title('Bias-Variance Trade-off: Training vs Test Error') plt.legend() plt.grid(True, alpha=0.3) plt.show()

使用sklearn的优势在于，PolynomialFeatures可以轻松地与其他预处理步骤（如标准化StandardScaler）和正则化模型（如Ridge,Lasso）结合，以应对过拟合问题，这是numpy.polyfit所不具备的。

5. 殊途同归？插值与拟合的边界与联合应用

经过前面的详细拆解，你可能觉得插值和拟合是两条平行线。但在实际工程中，它们的边界有时是模糊的，甚至可以协同工作。

场景一：从拟合到插值当你用多项式拟合一组数据后，得到的拟合函数本身就是一个连续函数，你当然可以用它来计算任意x对应的y值，这看起来像“插值”。但关键在于，拟合曲线通常不穿过原始数据点（除非次数足够高导致过拟合）。所以，对于数据点本身，拟合给出的是“估计值”，而插值给出的是“精确值”。

场景二：样条拟合我们之前讨论的三次样条是“插值样条”。还有一种叫“平滑样条”或“拟合样条”。它不再要求曲线必须穿过每一个点，而是在“拟合程度”（残差小）和“曲线光滑度”（二阶导数积分小）之间找一个平衡。这通过一个惩罚参数λ来控制：min Σ [y_i - S(x_i)]^2 + λ ∫ [S''(x)]^2 dx当λ=0，就是插值样条；当λ→∞，就趋向于一条直线（最光滑）。这其实是拟合思想在样条框架下的应用，非常适合处理带噪声的数据。

在MATLAB中，曲线拟合工具箱的spaps函数或csaps函数可以完成平滑样条拟合。在Python的SciPy中，scipy.interpolate.UnivariateSpline可以通过s参数控制平滑度。

场景三：先拟合趋势，再插值残差这是一个高级技巧。对于有明显趋势但又有局部波动的数据，可以先用一个低次多项式（或其它简单模型）拟合其“整体趋势”，然后对残差（原始数据减去趋势）进行样条插值。最后将趋势和插值后的残差相加。这样既能抓住大方向，又能用插值刻画局部细节。

# 示例：趋势拟合 + 残差插值 x = np.linspace(0, 10, 30) y_trend = 0.5 * x + 0.05 * x**2 # 真实的二次趋势 y_fluctuation = 2 * np.sin(2 * np.pi * x / 4) # 真实的周期性波动 y_true = y_trend + y_fluctuation y_noisy = y_true + 0.5 * np.random.randn(len(x)) # 加入噪声 # 1. 用二次多项式拟合趋势 coeffs_trend = np.polyfit(x, y_noisy, 2) p_trend = np.poly1d(coeffs_trend) y_trend_fitted = p_trend(x) # 2. 计算残差 residuals = y_noisy - y_trend_fitted # 3. 对残差进行三次样条插值（得到连续的函数） from scipy.interpolate import CubicSpline cs_residual = CubicSpline(x, residuals) # 4. 合成最终模型：在任何点x_new上，y = 趋势(x_new) + 残差插值(x_new) x_new = np.linspace(0, 10, 300) y_new_trend = p_trend(x_new) y_new_residual = cs_residual(x_new) y_new_total = y_new_trend + y_new_residual # 绘图 fig, axes = plt.subplots(2, 2, figsize=(12, 10)) axes[0, 0].plot(x, y_noisy, 'bo', label='Noisy Data') axes[0, 0].plot(x_new, p_trend(x_new), 'r-', linewidth=3, label='Fitted Trend (Poly deg=2)') axes[0, 0].set_title('Step 1: Fit the Trend') axes[0, 0].legend() axes[0, 0].grid(True) axes[0, 1].plot(x, residuals, 'go', label='Residuals') axes[0, 1].plot(x_new, y_new_residual, 'm-', label='Spline Interp of Residuals') axes[0, 1].axhline(y=0, color='k', linestyle='--', alpha=0.3) axes[0, 1].set_title('Step 2: Interpolate the Residuals') axes[0, 1].legend() axes[0, 1].grid(True) axes[1, 0].plot(x, y_noisy, 'bo', alpha=0.6, label='Noisy Data') axes[1, 0].plot(x_new, y_new_total, 'k-', linewidth=3, label='Combined Model') axes[1, 0].set_title('Step 3: Combined Result (Trend + Residual Spline)') axes[1, 0].legend() axes[1, 0].grid(True) # 对比：直接高次多项式拟合（可能过拟合） coeffs_direct = np.polyfit(x, y_noisy, 10) p_direct = np.poly1d(coeffs_direct) axes[1, 1].plot(x, y_noisy, 'bo', alpha=0.6, label='Noisy Data') axes[1, 1].plot(x_new, p_direct(x_new), 'r--', linewidth=2, label='Direct Poly Fit (deg=10)') axes[1, 1].plot(x_new, y_new_total, 'k-', linewidth=2, label='Trend+Spline') axes[1, 1].set_title('Comparison: Direct High-degree vs Our Method') axes[1, 1].legend() axes[1, 1].grid(True) plt.tight_layout() plt.show()

这种方法结合了拟合的“趋势捕捉能力”和插值的“局部细节刻画能力”，往往能获得比单一方法更稳健、更物理可解释的结果。