3

我想对数据集（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不仅接受权重，还返回系数的协方差矩阵......所以，问题解决了！:-)

4

2 回答 2

0

您想要一个快速加权最小二乘模型，无需额外开销即可返回协方差矩阵？通常，正确的协方差矩阵将取决于数据生成过程 (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 用于推导和讨论。

2012 年 12 月 4 日更新：还有另一个堆栈溢出问题即将到来： numpy.polyfit 没有关键字“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 网站。

于 2012-12-03T22:51:37.300 回答
0

这里使用 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]]
于 2012-12-04T19:45:01.610 回答