纽比.波利菲特给出有用的拟合，但协方差矩阵是无穷大的

Question

我试图用多项式拟合一组数据。有时，numpy.ployfit返回的协方差矩阵只由inf组成，尽管拟合似乎很有用。没有numpy.inf或'努比·楠“在数据里！在

示例：

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

结果：

^{pr2}$

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

---

区别应该在自由度上。在polyfit方法中，它已经考虑到你的学位是2，因此导致：

RuntimeWarning: divide by zero encountered in true_divide
  fac = resids / (len(x) - order - 2.0)

你可以传递你的np.cov一个ddof=关键字（ddof=delta自由度），你会遇到同样的问题

Answer 2

正如rustil的回答所说，这是由应用于协方差方程分母的偏差校正引起的，这将导致此输入的零除。这一更正背后的推理与Bessel's Correction背后的推理相似。这实际上是一个信号，表明有太少的数据点，以一个明确的方式估计协方差。在

如何避开这个问题？这个版本的权重。您可以添加另一个数据点，但要在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

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

^{pr2}$

我不确定这是否特别有意义，但这是解决问题的一种方法。不过，这有点糊涂。对于更健壮的方法，可以修改polyfit以接受它自己的ddof参数，也许可以代替cov当前接受的布尔值。（我只是opened an issue建议这么多。）

关于cov计算的最后一点提示：如果你看一下least squares regression上的维基百科页面，你会发现系数协方差的简化公式是inv(dot(dot(X, W), X))，至少粗略地说，numpy代码中有一个corresponding line。在本例中，X是Vandermonde matrix，权重已经是multiplied in。numpy代码也做一些缩放（我理解；这是最小化数值误差的策略的一部分）并将结果乘以残差的范数（我不明白；我只能猜测它是另一个版本的协方差公式的一部分）。在