需要一个Python多项式拟合函数，该函数返回协方差。

Question

需要一个Python多项式拟合函数，该函数返回协方差。

pythonnumpyscipycovariancepolynomial-math

4

我希望在数据集 (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 现在不仅接受权重，而且还返回系数的协方差矩阵... 因此，问题已解决！ :-)

- Rolf Bartstra

请查看mpfit或kmpfit。http://www.astro.rug.nl/software/kapteyn/kmpfit.html - reptilicus

根据链接，这是另一个（通用的）迭代求解器。由于速度的原因，我正在寻找分析（=非迭代）解决方案 - 这对于多项式来说是完全可能的。 - Rolf Bartstra

4

Statsmodels是什么？https://groups.google.com/forum/?fromgroups=#!topic/pystatsmodels/paCNa5sXbOo http://statsmodels.sourceforge.net/devel/generated/statsmodels.regression.linear_model.OLS.html - joris

@joris 这个链接可能确实有用，我会去看一下的。谢谢 - Rolf Bartstra

进行多项式拟合的代码非常简单…为什么不只是改编现有的numpy程序以同时返回协方差矩阵呢？ - Brian B

2个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- BKay · Answer 1

你想要一个快速的加权最小二乘模型，可以返回协方差矩阵而不需要额外的开销吗？一般来说，正确的协方差矩阵取决于数据生成过程（DGP），因为不同的DGP（比如误差异方差）意味着参数估计的不同分布（想想白噪声和OLS标准误之间的区别）。但是如果你可以假设WLS是正确的方法，我相信你会使用WLS的渐近方差估计值beta，(1/n X'V^-1X)^-1，其中V是从Yerrs创建的加权矩阵。如果numpy.polynomial.polynomial.polyfit对你有用，那么这是一个非常简单的公式。

我找了一个在线参考资料，但没有找到。但是请参见Fumio Hayashi的《经济计量学》，2000年，普林斯顿大学出版社，第133-137页的推导和讨论。

更新12/4/12：还有另一个与之相似的stackoverflow问题：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()

在使用numpy.polynomial.polynomial.polyfit时，您可以像以前一样将权重定义为Yerr的函数。有关使用WLS和statsmodels的更多详细信息，请访问statsmodels网站。

- sega_sai · 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]]