“`markdown

Python 数据分析：Numpy Polyfit 曲线拟合详解

在数据分析领域，我们经常会遇到需要从一组散点数据中找出其潜在趋势或规律的场景。曲线拟合就是这样一种强大的技术，它允许我们构建一个数学函数（曲线）来最好地描述这些数据点的行为。Python 的 numpy 库提供了 numpy.polyfit 函数，这是进行多项式曲线拟合的便捷且高效的工具。本文将深入探讨 numpy.polyfit 的工作原理、用法、参数选择以及一些高级应用。

1. 什么是曲线拟合？

曲线拟合是指选择一个函数模型（例如，直线、抛物线或其他多项式函数），并调整其参数，使其尽可能地与给定的数据点吻合。这个拟合出来的函数可以用于：

• 趋势预测： 根据现有数据预测未来值。
• 数据平滑： 消除数据中的噪声，揭示底层模式。
• 模型构建： 理解变量之间的关系。
• 插值与外推： 在已知数据点之间或之外估计值。

2. numpy.polyfit 简介

numpy.polyfit 是 NumPy 库中用于执行最小二乘多项式拟合的函数。它能够找到一个指定次数的多项式，使得该多项式曲线与给定数据点的距离平方和最小。

函数签名

python
numpy.polyfit(x, y, deg, rcond=None, full=False, w=None, cov=False)

主要参数解释

• x：一维数组，表示数据点的 x 坐标。
• y：一维数组，表示数据点的 y 坐标，与 x 的长度相同。
• deg：整数，表示拟合多项式的次数。例如，deg=1 表示拟合一条直线 (ax + b)，deg=2 表示拟合一个二次多项式 (ax^2 + bx + c)。
• w (可选)：与 y 相同长度的数组，表示每个数据点的权重。如果提供，拟合将最小化加权平方误差。这在某些数据点比其他数据点更可靠时非常有用。
• full (可选)：布尔值。如果为 True，除了系数外，还会返回残差和一些其他诊断信息。
• cov (可选)：布尔值。如果为 True，并且 full 为 False，则返回系数的协方差矩阵，可用于估计系数的不确定性。

返回值

numpy.polyfit 主要返回多项式系数。对于 deg 次多项式 p(x) = c_0 * x^deg + c_1 * x^(deg-1) + ... + c_deg，它将返回一个长度为 deg + 1 的数组 [c_0, c_1, ..., c_deg]。请注意，系数是按照降幂顺序排列的。

3. numpy.polyfit 工作原理：最小二乘法

numpy.polyfit 的核心是最小二乘法 (Least Squares Method)。给定一组数据点 (x_i, y_i) 和一个指定次数 deg 的多项式 P(x) = a_n * x^n + ... + a_1 * x + a_0，最小二乘法的目标是找到一组系数 a_n, ..., a_0，使得所有数据点到多项式曲线的垂直距离的平方和最小。

数学上，就是要最小化以下残差平方和 (Sum of Squared Residuals, SSR)：

$$ \text{SSR} = \sum_{i=1}^{N} (y_i – P(x_i))^2 $$

通过对每个系数求偏导并令其为零，可以得到一组线性方程组，解出这些方程组即可得到最佳拟合多项式的系数。numpy.polyfit 在底层高效地完成了这一计算。

4. numpy.polyfit 的基本使用与可视化

让我们通过一个简单的例子来演示如何使用 numpy.polyfit 并可视化拟合结果。

“`python
import numpy as np
import matplotlib.pyplot as plt

1. 生成一些示例数据

模拟一个二次曲线加上一些随机噪声

np.random.seed(42) # 为了结果可复现
x = np.linspace(0, 10, 50)
y_true = 0.5 * x**2 – 3 * x + 5
y_noisy = y_true + np.random.normal(0, 5, x.shape) # 添加噪声

2. 使用 polyfit 进行多项式拟合

拟合一个二次多项式 (deg=2)

coefficients = np.polyfit(x, y_noisy, deg=2)
print(“拟合的多项式系数 (降幂):”, coefficients)

3. 构建拟合多项式函数

numpy.poly1d 可以很方便地根据系数构建多项式函数

p = np.poly1d(coefficients)
print(“拟合的多项式函数:”, p)

4. 计算拟合曲线的y值

y_fitted = p(x)

5. 可视化结果

plt.figure(figsize=(10, 6))
plt.scatter(x, y_noisy, label=’噪声数据’, alpha=0.7)
plt.plot(x, y_true, color=’green’, linestyle=’–‘, label=’真实曲线’)
plt.plot(x, y_fitted, color=’red’, label=f’拟合曲线 (deg={2}): {p}’)
plt.title(‘NumPy Polyfit 曲线拟合示例’)
plt.xlabel(‘X’)
plt.ylabel(‘Y’)
plt.legend()
plt.grid(True)
plt.show()
“`

运行上述代码，你将看到一个散点图，以及真实曲线和 polyfit 拟合出的二次多项式曲线。拟合曲线会尽可能地穿过噪声数据点，并接近真实曲线。

5. 如何选择多项式次数 (deg)？

选择合适的多项式次数 deg 是曲线拟合中的关键一步。不当的 deg 会导致：

• 欠拟合 (Underfitting)：如果 deg 太小，模型过于简单，无法捕捉数据的真实模式。拟合曲线会显得过于平坦，不能很好地通过数据点。
• 过拟合 (Overfitting)：如果 deg 太大（特别是当 deg 接近或超过数据点数量时），模型会变得过于复杂，它不仅拟合了数据的真实模式，还拟合了数据中的随机噪声。这会导致模型在训练数据上表现良好，但在新数据上表现很差。过拟合的曲线通常会呈现出剧烈的波动。

选择 deg 的策略

1. 领域知识：如果对数据背后物理过程有了解，可能知道函数大致的形式（例如，线性关系、二次关系等）。
2. 可视化检查：这是最直观的方法。尝试不同的 deg 值，并绘制拟合曲线。观察哪条曲线既平滑又很好地捕捉了数据的整体趋势，而没有过度跟随噪声。
3. 残差分析：计算 y_noisy - y_fitted 作为残差。理想情况下，残差应随机分布在零附近，没有明显的模式。如果残差呈现出某种模式（例如，U形或倒U形），可能意味着多项式次数过低。
4. 交叉验证 (Cross-validation)：将数据分为训练集和验证集。在训练集上拟合不同次数的多项式，然后在验证集上评估它们的性能（例如，使用均方误差）。选择在验证集上表现最好的模型。
5. 赤池信息准则 (AIC) 或贝叶斯信息准则 (BIC)：这些统计量可以帮助在模型复杂性和拟合优度之间取得平衡，选择更优的模型。

示例：不同 deg 的影响

“`python

… (保留上面生成数据的代码)

degrees = [1, 2, 5, 10] # 尝试不同次数的多项式

plt.figure(figsize=(12, 8))
plt.scatter(x, y_noisy, label=’噪声数据’, alpha=0.7, zorder=10)
plt.plot(x, y_true, color=’green’, linestyle=’–‘, label=’真实曲线’, linewidth=2)

for deg in degrees:
coefficients = np.polyfit(x, y_noisy, deg=deg)
p = np.poly1d(coefficients)
y_fitted = p(x)
plt.plot(x, y_fitted, label=f’拟合曲线 (deg={deg})’, linestyle=’-‘)

plt.title(‘不同多项式次数对拟合结果的影响’)
plt.xlabel(‘X’)
plt.ylabel(‘Y’)
plt.legend()
plt.grid(True)
plt.ylim(min(y_noisy) – 10, max(y_noisy) + 10) # 调整Y轴范围
plt.show()
“`

从图中可以看到：
* deg=1 (直线) 明显欠拟合，未能捕捉二次趋势。
* deg=2 (二次) 拟合效果最好，因为它与真实曲线的次数相同。
* deg=5 开始出现轻微过拟合的迹象，曲线开始试图穿过更多的噪声点。
* deg=10 严重过拟合，曲线波动剧烈，在数据点之间表现出不自然的振荡。

6. 高级用法与考虑

加权拟合 (w 参数)

在某些情况下，数据的测量精度可能不同。例如，某些数据点可能比其他数据点更可靠。w 参数允许我们为每个数据点指定一个权重，使得 polyfit 在拟合时更看重权重高的数据点。

“`python

假设前25个数据点更可靠，赋予更高权重

weights = np.ones_like(x)
weights[:25] = 2.0 # 前半部分数据权重更高

coefficients_weighted = np.polyfit(x, y_noisy, deg=2, w=weights)
p_weighted = np.poly1d(coefficients_weighted)
y_fitted_weighted = p_weighted(x)

plt.figure(figsize=(10, 6))
plt.scatter(x, y_noisy, label=’噪声数据’, alpha=0.7)
plt.plot(x, y_true, color=’green’, linestyle=’–‘, label=’真实曲线’)
plt.plot(x, y_fitted, color=’red’, label=’未加权拟合曲线 (deg=2)’)
plt.plot(x, y_fitted_weighted, color=’purple’, linestyle=’:’, label=’加权拟合曲线 (deg=2)’)
plt.title(‘加权拟合与未加权拟合对比’)
plt.xlabel(‘X’)
plt.ylabel(‘Y’)
plt.legend()
plt.grid(True)
plt.show()
“`

加权拟合会使曲线更倾向于靠近权重较大的数据点。

返回残差与协方差 (full 和 cov 参数)

• full=True: 返回更多诊断信息，包括残差平方和、秩、奇异值以及条件数。这对于评估拟合质量很有用。

  python
coefficients_full, residuals, rank, singular_values, rcond = np.polyfit(x, y_noisy, deg=2, full=True)
print("\n拟合残差平方和:", residuals[0]) # 第一个元素是残差平方和

• cov=True: 返回系数的协方差矩阵。协方差矩阵的对角线元素是各个系数的方差，其平方根是标准误差，可用于构建置信区间。

  “`python
  coefficients_cov, cov_matrix = np.polyfit(x, y_noisy, deg=2, cov=True)
  print(“\n系数的协方差矩阵:\n”, cov_matrix)

  系数标准误差

  std_errors = np.sqrt(np.diag(cov_matrix))
  print(“系数的标准误差:”, std_errors)
  “`

  协方差信息对于量化拟合参数的不确定性非常重要。

7. numpy.polyfit 的局限性与替代方案

numpy.polyfit 在处理多项式关系时非常强大和高效。然而，它也有其局限性：

• 仅限于多项式函数：polyfit 只能拟合多项式。如果数据的潜在关系是非多项式（例如指数函数、对数函数、三角函数等），则 polyfit 无法直接拟合。
• 高次多项式的数值稳定性：当多项式次数很高时，拟合过程可能会出现数值不稳定问题，导致结果不准确或振荡。

替代方案

对于非多项式曲线拟合，或者需要更灵活地定义拟合函数，可以使用 scipy.optimize.curve_fit：

• scipy.optimize.curve_fit：这是 SciPy 库中一个更通用的非线性最小二乘拟合工具。它允许用户自定义任何形式的函数来拟合数据。你只需提供一个接受 x 和模型参数作为输入的函数，以及初始参数猜测值。

  “`python
  from scipy.optimize import curve_fit

  定义一个指数函数模型

  def exponential_func(x, a, b, c):
  return a * np.exp(b * x) + c

  假设有一些指数增长的数据

  x_exp = np.linspace(0, 4, 50)
  y_exp_true = 2 * np.exp(0.5 * x_exp) + 1
  y_exp_noisy = y_exp_true + np.random.normal(0, 0.5, x_exp.shape)

  使用 curve_fit 拟合

  提供一个初始猜测值 (p0) 非常重要

  params, covariance = curve_fit(exponential_func, x_exp, y_exp_noisy, p0=[1, 0.1, 1])

  y_fitted_exp = exponential_func(x_exp, *params)

  plt.figure(figsize=(10, 6))
  plt.scatter(x_exp, y_exp_noisy, label=’噪声数据’)
  plt.plot(x_exp, y_exp_true, color=’green’, linestyle=’–‘, label=’真实指数曲线’)
  plt.plot(x_exp, y_fitted_exp, color=’red’, label=f’拟合指数曲线’)
  plt.title(‘SciPy curve_fit 指数函数拟合示例’)
  plt.xlabel(‘X’)
  plt.ylabel(‘Y’)
  plt.legend()
  plt.grid(True)
  plt.show()
  “`

8. 总结

numpy.polyfit 是 Python 中进行多项式曲线拟合的基石。它简单易用，基于强大的最小二乘法原理，能够有效地从噪声数据中提取出多项式趋势。理解其参数（尤其是 deg 的选择）和工作原理对于成功进行数据分析至关重要。当遇到非多项式关系时，scipy.optimize.curve_fit 提供了一个更通用的解决方案。掌握这些工具，你将能够更深入地理解和利用数据中的隐藏模式。
“`