如何通过numpy.polynomial.polynomial.polyfit获取拟合参数协方差矩阵

numpy.polynomial.polynomial.polyfit确实没有内置的cov参数直接返回协方差矩阵，但可以通过纯NumPy的手动计算得到，完全不需要依赖SciPy。核心是利用最小二乘拟合的协方差公式：

协方差矩阵 = 残差方差 × (XᵀX)⁻¹

其中：

• X是设计矩阵，对应多项式拟合的基函数（比如直线拟合时，X的第一列全为1，第二列为自变量x）
• 残差方差 = 残差平方和 / (数据点数量 - 参数数量)

具体步骤与代码示例 #

以直线拟合（一次多项式）为例：

1. 生成模拟数据

import numpy as np

# 模拟带噪声的直线数据
x = np.linspace(0, 10, 50)
y_true = 2 * x + 3
y = y_true + np.random.normal(0, 1, size=x.shape)

2. 用numpy.polynomial.polynomial.polyfit拟合系数

from numpy.polynomial.polynomial import polyfit, polyval

# 拟合一次多项式（直线），得到系数[a0, a1]，对应y = a0 + a1*x
coeffs = polyfit(x, y, deg=1)
a0, a1 = coeffs

3. 手动计算协方差矩阵

n = len(x)
params_num = len(coeffs)  # 直线拟合是2个参数

# 计算拟合值与残差
y_fit = polyval(x, coeffs)
residuals = y - y_fit

# 计算残差方差（自由度校正）
residual_variance = np.sum(residuals**2) / (n - params_num)

# 构造设计矩阵X：第一列全1（对应常数项），第二列是x（对应一次项）
X = np.column_stack((np.ones(n), x))

# 计算(XᵀX)的逆矩阵
XTX_inv = np.linalg.inv(X.T @ X)

# 计算协方差矩阵
cov_matrix = residual_variance * XTX_inv

4. 验证结果（对比传统numpy.polyfit的cov=True输出）

# 用旧版polyfit得到协方差做对比
old_coeffs, old_cov = np.polyfit(x, y, deg=1, cov=True)
# 注意旧版polyfit返回的系数是[a1, a0]（高次在前），所以协方差矩阵的顺序也要对应调整
old_cov_adjusted = old_cov[::-1, ::-1]

print("手动计算的协方差矩阵：")
print(cov_matrix)
print("\n旧版polyfit的协方差矩阵（调整顺序后）：")
print(old_cov_adjusted)

扩展到更高次多项式 #

如果是拟合更高次的多项式，只需要修改设计矩阵X的构造：比如二次多项式y = a0 + a1*x + a2*x²，X的列就是[1, x, x²]，其余计算逻辑完全一致。

为什么这个方法可行？ #

numpy.polynomial.polynomial.polyfit本质上也是用最小二乘法求解，所以手动推导的协方差计算逻辑和旧版numpy.polyfit完全一致，只是新库把协方差计算的步骤暴露给用户，而非封装成参数。这种手动计算的方式也能让大一物理系学生更直观理解最小二乘的统计意义，适合教学场景。

内容的提问来源于stack exchange，提问作者Ben