需要返回协方差的 Python 多项式拟合函数

Question

我想对数据集 (X,Y,Yerr) 进行最小二乘多项式拟合并获得拟合参数的协方差矩阵。此外，由于我有很多数据集，CPU 时间是一个问题，所以我正在寻找一个分析（=快速）解决方案。我发现了以下（非理想）选项：

numpy.polyfit 进行拟合，但不考虑错误 Yerr，也不返回协方差；

numpy.polynomial.polynomial.polyfit 接受 Yerr 作为输入（以权重的形式），但也不返回协方差；

scipy.optimize.curve_fit 和 scipy.optimize.leastsq 可以定制以拟合多项式并返回协方差矩阵，但是 - 作为迭代方法 - 这些比 polyfit 例程（产生解析解）慢得多；

Python 是否提供了一个分析多项式拟合例程来返回拟合参数的协方差（或者我是否必须自己编写一个:-)？

更新： 似乎在 Numpy 1.7.0 中，numpy.polyfit 现在不仅确实接受权重，而且还返回系数的协方差矩阵......所以，问题解决了！ :-)

Answer 1

您想要一个快速加权最小二乘模型，该模型无需额外开销即可返回协方差矩阵？通常，正确的协方差矩阵将取决于数据生成过程 (DGP)，因为不同的 DGP（例如误差的异方差）意味着参数估计的不同分布（想想 White 与 OLS 标准误差）。但是，如果您可以假设 WLS 是正确的方法，并且我相信您会使用 WLS 的 beta 的渐近方差估计 (1/n X'V^-1X)^-1，其中 V 是加权矩阵由 Yerrs 创建。如果 numpy.polynomial.polynomial.polyfit 为您工作，这是一个非常简单的公式。

我查找了在线参考，但找不到。但请参阅 Fumio Hayashi 的 Ecomometrics，2000 年，普林斯顿大学出版社，p。 133 - 137 进行推导和讨论。

12/4/12 更新： 还有另一个接近的堆栈溢出问题： numpy.polyfit has no keyword 'cov' 很好地解释了如何使用 scikits.statsmodels 来做你想做的事。我相信您会想要替换该行：

result = sm.OLS(Y,reg_x_data).fit()

到

result = sm.WLS(Y,reg_x_data, weights).fit()

您将权重定义为 Yerr 的函数，就像之前使用 numpy.polynomial.polynomial.polyfit 一样。有关将 statsmodels 与 WLS 结合使用的更多详细信息，请参见 statsmodels website。

Answer 2

这里使用的是 scipy.linalg.lstsq

import numpy as np,numpy.random, scipy.linalg
#generate the test data
N = 100
xs = np.random.uniform(size=N)
errs = np.random.uniform(0, 0.1, size=N) # errors
ys = 1 + 2 * xs + 3 * xs ** 2 + errs * np.random.normal(size=N)

# do the fit
polydeg = 2
A = np.vstack([1 / errs] + [xs ** _ / errs for _ in range(1, polydeg + 1)]).T
result = scipy.linalg.lstsq(A, (ys / errs))[0]
covar = np.matrix(np.dot(A.T, A)).I
print result, '\n', covar

>> [ 0.99991811  2.00009834  3.00195187]
[[  4.82718910e-07  -2.82097554e-06   3.80331414e-06]
 [ -2.82097554e-06   1.77361434e-05  -2.60150367e-05]
 [  3.80331414e-06  -2.60150367e-05   4.22541049e-05]]