python - numpy.polyfit 提供有用的拟合，但协方差矩阵无限大

Question

我正在尝试将多项式拟合到一组数据中。有时可能会发生 numpy.ployfit 返回的协方差矩阵只由 inf 组成，虽然 fit 看起来很有用。数据中没有 numpy.inf 或 'numpy.nan'!

例子:

import numpy as np
# sample data, does not contain really x**2-like behaviour, 
# but that should be visible in the fit results
x = [-449., -454., -459., -464., -469.]
y = [ 0.9677024,   0.97341953,  0.97724978,  0.98215678,  0.9876293]

fit, cov = np.polyfit(x, y, 2, cov=True)

print 'fit: ', fit
print 'cov: ', cov

结果:

fit: [  1.67867158e-06   5.69199547e-04   8.85146009e-01]
cov: [[ inf  inf  inf]
      [ inf  inf  inf]
      [ inf  inf  inf]]

np.cov(x,y) 给出

[[  6.25000000e+01  -6.07388099e-02]
 [ -6.07388099e-02   5.92268942e-05]]

因此 np.cov 与 np.polyfit 返回的协方差不同。有人知道发生了什么事吗？

编辑: 我现在明白了 numpy.cov 不是我想要的。我需要多项式系数的方差，但如果 (len(x) - order - 2.0) == 0，我就得不到它们。还有其他方法可以得到拟合多项式系数的方差吗？

Answer 1

作为rustil的回答说，这是由应用于协方差方程的分母的偏差校正引起的，这导致该输入的除法为零。此更正背后的原因与 Bessel's Correction 背后的原因类似。 .这确实表明数据点太少，无法以明确定义的方式估计协方差。

如何解决这个问题？好吧，这个版本的 polyfit 接受权重。您可以添加另一个数据点，但以 epsilon 对其进行加权。这相当于减少了this formula中的2.0到 1.0。

x = [-449., -454., -459., -464., -469.]
y = [ 0.9677024,   0.97341953,  0.97724978,  0.98215678,  0.9876293]

x_extra = x + x[-1:]
y_extra = y + y[-1:]
weights = [1.0, 1.0, 1.0, 1.0, 1.0, sys.float_info.epsilon]

fit, cov = np.polyfit(x, y, 2, cov=True)
fit_extra, cov_extra = np.polyfit(x_extra, y_extra, 2, w=weights, cov=True)

print fit == fit_extra
print cov_extra

输出。请注意，拟合值是相同的:

>>> print fit == fit_extra
[ True  True  True]
>>> print cov_extra
[[  8.84481850e-11   8.11954338e-08   1.86299297e-05]
 [  8.11954338e-08   7.45405039e-05   1.71036963e-02]
 [  1.86299297e-05   1.71036963e-02   3.92469307e+00]]

我非常不确定这是否特别有意义，但这是解决问题的一种方法。不过，这有点麻烦。对于更健壮的东西，您可以修改 polyfit 以接受它自己的 ddof 参数，也许可以代替 cov 当前接受的 bool 值。 (我只是 opened an issue 提出了同样多的建议。)

关于 cov 计算的最后一点快速说明:如果您查看 least squares regression 上的维基百科页面，您会看到系数协方差的简化公式是 inv(dot(dot(X, W), X))，它有一个 corresponding line在 numpy 代码中——至少粗略地说。在这种情况下，X 是 Vandermonde matrix , 权重已经是 multiplied in . numpy 代码还进行了一些缩放(我理解；这是最小化数值误差的策略的一部分)并将结果乘以残差的范数(我不理解；我只能猜测它是另一个版本的一部分协方差公式)。