scipy.curve_fit 与 numpy.polyfit 不同的协方差矩阵答案

【问题标题】：scipy.curve_fit vs. numpy.polyfit different covariance matricesscipy.curve_fit 与 numpy.polyfit 不同的协方差矩阵
【发布时间】：2018-08-23 17:40:51
【问题描述】：

我正在使用 Python 3.6 进行数据拟合。最近遇到以下问题，经验不足，不知道怎么处理。

如果我在同一组数据点上使用 numpy.polyfit(x, y, 1, cov=True) 和 scipy.curve_fit(lambda: x, a, b: a*x+b, x, y) ，我得到几乎相同的系数 a 和 b。但是 scipy.curve_fit 的协方差矩阵的值大约是 numpy.polyfit 的值的一半。

由于我想使用协方差矩阵的对角线来估计系数的不确定性（u = numpy.sqrt(numpy.diag(cov))），所以我有三个问题：

哪个协方差矩阵是正确的（我应该使用哪个）？
为什么会有差异？
需要什么才能使它们相等？

谢谢！

编辑：

import numpy as np
import scipy.optimize as sc

data = np.array([[1,2,3,4,5,6,7],[1.1,1.9,3.2,4.3,4.8,6.0,7.3]]).T

x=data[:,0]
y=data[:,1]

A=np.polyfit(x,y,1, cov=True)
print('Polyfit:', np.diag(A[1]))

B=sc.curve_fit(lambda x,a,b: a*x+b, x, y)
print('Curve_Fit:', np.diag(B[1]))

如果我使用statsmodels.api，结果对应于curve_fit。

【问题讨论】：

欢迎来到 Stack Overflow！您能为我们准备两个可见的最小示例吗？看看minimal reproducible example。（拟合 3 或 4 个点就足够了）

标签： python numpy scipy curve-fitting

【解决方案1】：

我想这与this有关

593          # Some literature ignores the extra -2.0 factor in the denominator, but 
594          #  it is included here because the covariance of Multivariate Student-T 
595          #  (which is implied by a Bayesian uncertainty analysis) includes it. 
596          #  Plus, it gives a slightly more conservative estimate of uncertainty. 
597          if len(x) <= order + 2: 
598              raise ValueError("the number of data points must exceed order + 2 " 
599                               "for Bayesian estimate the covariance matrix") 
600          fac = resids / (len(x) - order - 2.0) 
601          if y.ndim == 1: 
602              return c, Vbase * fac 
603          else: 
604              return c, Vbase[:,:, NX.newaxis] * fac

在这种情况下，len(x) - order 是 4，(len(x) - order - 2.0) 是 2，这可以解释为什么您的值相差 2 倍。

这解释了问题 2。问题 3 的答案可能是“获取更多数据。”，对于较大的len(x)，差异可能可以忽略不计。

哪个公式正确（问题 1）可能是 Cross Validated 的问题，但我认为它是 curve_fit，因为它明确旨在计算您所说的不确定性.来自documentation

pcov：二维数组

popt 的估计协方差。对角线提供参数估计的方差。要计算参数的一个标准差误差，请使用 perr = np.sqrt(np.diag(pcov))。

虽然上面polyfit 的代码中的注释表明它的意图更多是用于学生-T 分析。

【讨论】：

【解决方案2】：

这两种方法以不同的方式计算协方差。我不太确定polyfit 使用的方法，但curve_fit 通过反转 J.T.dot(J) 来估计协方差矩阵，其中 J 是模型的雅可比矩阵。通过查看polyfit 的代码，他们似乎反转了 lhs.T.dot(lhs)，其中 lhs 被定义为 Vandermonde 矩阵，尽管我不得不承认我不知道第二种方法的数学背景。

现在，至于你的问题是正确的，polyfit 的代码有以下注释：

# Some literature ignores the extra -2.0 factor in the denominator, but
#  it is included here because the covariance of Multivariate Student-T
#  (which is implied by a Bayesian uncertainty analysis) includes it.
#  Plus, it gives a slightly more conservative estimate of uncertainty.

基于此和您的观察，polyfit 似乎总是给出比curve_fit 更大的估计值。这是有道理的，因为J.T.dot(J) 是协方差矩阵的一阶近似值。因此，当有疑问时，高估错误总是更好。

但是，如果您知道数据中的测量误差，我建议您也提供它们并使用absolute_sigma=True 致电curve_fit。在我自己的tests 看来，这样做确实符合人们预期的分析结果，所以我很想知道在提供测量误差时两者中哪一个表现更好。

【讨论】：