我很难弄清楚 numpy polyfit 中协方差矩阵的缩放比例。
在文档中，我读到从未缩放到缩放协方差矩阵的比例因子为
chi2 / sqrt(N - DOF).
在下面附加的代码中，似乎比例因子实际上是
chi2 / DOF
这是我的代码
# Generate synthetically the data
# True parameters
import numpy as np
true_slope = 3
true_intercept = 7
x_data = np.linspace(-5, 5, 30)
# The y-data will have a noise term, to simulate imperfect observations
sigma = 1
y_data = true_slope * np.linspace(-5, 5, 30) + true_intercept
y_obs = y_data + np.random.normal(loc=0.0, scale=sigma, size=x_data.size)
# Here I generate artificially some unequal uncertainties
# (even if there is no reason for them to be so)
y_uncertainties = sigma * np.random.normal(loc=1.0, scale=0.5*sigma, size=x_data.size)
# Make the fit
popt, pcov = np.polyfit(x_data, y_obs, 1, w=1/y_uncertainties, cov='unscaled')
popt, pcov_scaled = np.polyfit(x_data, y_obs, 1, w=1/y_uncertainties, cov=True)
my_scale_factor = np.sum((y_obs - popt[0] * x_data - popt[1])**2 / y_uncertainties**2)
/ (len(y_obs)-2)
scale_factor = pcov_scaled[0,0] / pcov[0,0]
如果我运行代码，我看到实际的比例因子是 chi2/DOF，而不是文档中报告的值。这是真的还是我错过了什么？
我还有一个问题。在不确定性呈正态分布的情况下，为什么建议仅使用 y 数据误差的倒数而不是 y 数据误差的倒数作为权重？
编辑以添加代码运行生成的数据
x_data = array([-5. , -4.65517241, -4.31034483, -3.96551724, -3.62068966,
-3.27586207, -2.93103448, -2.5862069 , -2.24137931, -1.89655172,
-1.55172414, -1.20689655, -0.86206897, -0.51724138, -0.17241379,
0.17241379, 0.51724138, 0.86206897, 1.20689655, 1.55172414,
1.89655172, 2.24137931, 2.5862069 , 2.93103448, 3.27586207,
3.62068966, 3.96551724, 4.31034483, 4.65517241, 5. ])
y_obs = array([-7.27819725, -8.41939411, -3.9089926 , -5.24622589, -3.78747379,
-1.92898727, -1.375255 , -1.84388812, -0.37092441, 0.27572306,
2.57470918, 3.860485 , 4.62580789, 5.34147103, 6.68231985,
7.38242258, 8.28346559, 9.46008873, 10.69300274, 12.46051285,
13.35049975, 13.28279961, 14.31604781, 16.8226239 , 16.81708308,
18.64342284, 19.37375515, 19.6714002 , 20.13700708, 22.72327533])
y_uncertainties = array([ 0.63543112, 1.07608924, 0.83603265, -0.03442888, -0.07049299,
1.30864191, 1.36015322, 1.42125414, 1.04099854, 1.20556608,
0.43749964, 1.635056 , 1.00627014, 0.40512511, 1.19638787,
1.26230966, 0.68253139, 0.98055035, 1.01512232, 1.83910276,
0.96763007, 0.57373151, 1.69358475, 0.62068133, 0.70030971,
0.34648312, 1.85234844, 1.18687269, 1.23841579, 1.19741206])
有了这些数据，我得到了scale_factor = 1.6534129347542432，my_scale_factor = 1.653412934754234和文档中报告的"名义"比例因子，即
nominal_scale_factor = np.sum((y_obs - popt[0] * x_data - popt[1])**2 /
y_uncertainties**2) / np.sqrt(len(y_obs) - len(y_obs) + 2)
具有价值nominal_scale_factor = 32.73590595145554
我的 numpy 版本是1.18.5 3.7.7 (default, May 6 2020, 11:45:54) [MSC v.1916 64 bit (AMD64)]
关于numpy.polyfit文档：
默认情况下，协方差按 chi2/sqrt(N-dof) 缩放，即，权重被假定为不可靠，除非在相对意义上，并且所有内容都缩放，使得简化的 chi2 是单位。
这看起来像一个文档错误。协方差的正确比例因子chi_square/(N-M)其中M是拟合参数的数量，N-M是自由度的数量。看起来np.polyfit实现了正确，因为my_scale_factor和scale_factor是一致的。
关于为什么不是"y数据误差的逆平方"的问题：多项式拟合或更一般地说，最小二乘拟合涉及求解p向量
A @ p = y
其中A是y中N数据点和p中M元素的(N, M)矩阵，A中的每一列是在相应x值处计算的多项式项。
该解决方案最大限度地减少了
(SUM_j A[i, j] p[j] - y[i])^2
SUM -----------------------------
i sigma_y[i]^2
在计算上，最便宜的计算方法是将A中的每一行和每个y值乘以相应的1/sigma_y，然后取A@p=y方程的标准最小二乘解。通过让用户提供反误差，它可以避免拟合例程处理除以零问题和慢速平方根运算。
关于第一部分，我打开了一个 Github 问题
https://github.com/numpy/numpy/issues/16842
该线程的结论是文档是错误的，但函数行为正确。
文档应更新为
默认情况下，协方差按chi2/dof缩放，即，除非在相对意义上，否则权重被假定为不可靠，并且所有内容都缩放，使得简化的 chi2 是单位。