numpy.polyfit与scipy.odr的区别

Question

我有一组数据集，理论上应该由二次多项式描述。我想要拟合这些数据，使用了numpy.polyfit方法。但是，不幸的是，返回系数的误差不可用。因此，我决定也使用scipy.odr方法来拟合数据。奇怪的是，多项式的系数相互偏离。

我不明白这一点，因此决定在我自己生成的一组数据上测试两种拟合程序：

import numpy
import scipy.odr
import matplotlib.pyplot as plt

x = numpy.arange(-20, 20, 0.1)
y = 1.8 * x**2 -2.1 * x + 0.6 + numpy.random.normal(scale = 100, size = len(x))

#Define function for scipy.odr
def fit_func(p, t):
  return p[0] * t**2 + p[1] * t + p[2]

#Fit the data using numpy.polyfit
fit_np = numpy.polyfit(x, y, 2)

#Fit the data using scipy.odr
Model = scipy.odr.Model(fit_func)
Data = scipy.odr.RealData(x, y)
Odr = scipy.odr.ODR(Data, Model, [1.5, -2, 1], maxit = 10000)
output = Odr.run()
#output.pprint()
beta = output.beta
betastd = output.sd_beta

print "poly", fit_np
print "ODR", beta

plt.plot(x, y, "bo")
plt.plot(x, numpy.polyval(fit_np, x), "r--", lw = 2)
plt.plot(x, fit_func(beta, x), "g--", lw = 2)

plt.tight_layout()

plt.show()

一个示例的结果如下：

poly [ 1.77992643 -2.42753714  3.86331152]
ODR [   3.8161735   -23.08952492 -146.76214989]

在这张图片中，numpy.polyfit的解决方案（红色虚线）相当好。而scipy.odr的解决方案（绿色虚线）完全偏离了。我必须指出，在我要拟合的实际数据集中，numpy.polyfit和scipy.odr之间的差异较小。然而，我不明白两者之间的差异来自何处，为什么在我的测试示例中差异非常大，以及哪种拟合程序更好？我希望你能提供答案，帮助我更好地理解这两个拟合程序，并在此过程中回答我所提出的问题。

- The Dude

Answer 1

在您使用ODR时，它会进行完整的正交距离回归。要让它执行普通的非线性最小二乘拟合，请添加：

Odr.set_job(fit_type=2)

在开始优化之前，您将获得您期望的结果。

完整ODR失败的原因是由于没有指定权重/标准偏差。很明显，它难以解释点云，并假定x和y具有相等的权重。如果您提供估计的标准偏差，ODR也会产生一个好的（尽管不同）结果。

Data = scipy.odr.RealData(x, y, sx=0.1, sy=10)

- Christian K.

---

当我这样做时，我会得到一个分段错误。此外，如果我像http://docs.scipy.org/doc/scipy/reference/generated/scipy.odr.ODR.html建议的那样将“job”作为Odr选项包含。 - The Dude

@The Dude：我犯了一个错误，应该是fit_type=2，我刚刚编辑了我的帖子。请再试一次。fit_type=1用于隐式拟合。 - Christian K.

可以的！你有关于fit_type可能取值的不同整数的文档链接吗？ - The Dude

http://docs.scipy.org/doc/scipy/reference/generated/scipy.odr.ODR.set_job.html#scipy.odr.ODR.set_job - Christian K.

Answer 2

实际问题是odr输出的beta系数顺序与numpy.polyfit相反。因此，绿色曲线计算不正确。要绘制它，请改用：

plt.plot(x, fit_func(beta[::-1], x), "g--", lw = 2)

- kirkt

---

不，polyfit 首先返回最高次幂系数，就像在他的 fit_func() 中定义的一样。你说得对，如果有人弄错了，这可能会毁掉一切。 - EL_DON