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

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

numpy.polyfit是否合适，但没有考虑错误Yerr，也不返回协方差;

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

scipy.optimize.curve_fit和scipy.optimize.leastsq可以被定制，以适应多项式和返回的协方差矩阵，但 - 是迭代方法 - 这些是比polyfit例程（其产生的解析解）慢得多;

Python提供了一个解析多项式拟合程序，它返回拟合参数的协方差（或者我必须自己写一个:-)？

更新： 看来，与NumPy 1.7.0，numpy.polyfit现在不仅不接受权，而且还返回系数的协方差矩阵。所以，问题解决了！ :-)

来源

2012-10-08 Rolf Bartstra

查找到mpfit或kmpfit。 http://www.astro.rug.nl/software/kapteyn/kmpfit.html – reptilicus

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

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

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

我查找了一个在线参考，但找不到一个。但请参阅林雄二的“计量经济学”，2000年，普林斯顿大学出版社， 133 - 137进行推导和讨论。

更新12年12月4日：有一种接近另一个堆栈溢出问题： 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之前。更多关于在WLS结束时使用statsmodels的更多详情，请致电 statsmodels website。

来源

2012-12-03 22:51:37 BKay

Thnx，我知道做计算的公式，我只是希望相应的代码已经在Python/Numpy中实现 - 这似乎不是这种情况:-( –

这是使用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-04 19:45:01

谢谢，这将工作正常一个数据集，甚至多个集合，只要每个集合中的错误都是相同的，然而，一般来说，对于不同的集合，错误可能是不同的，在每种情况下，产生不同的矩阵A'linalg.lstsq'算法然后需要放置在一个循环中 - 这正是我不想要的（因为计算速度）。在这种一般情况下，解决方案可能在一个巨大的数组操作中计算，这将大大加快速度。因为我知道这样的函数不存在（即：但是 - 因为我要自己构建它） –

如果您将有不同的数据集，您必须再次构建矩阵（无论如何这是非常轻量级的操作），并解决sy再干。没有其他办法。我认为将不同的卡方问题组合成一个大问题是没有任何好处的，因为矩阵计算的性能将会是〜N^2，所以你可以更好地解决多个小问题，而不是一个大问题很多参数。 –

你是对的，但将不同的卡方问题合并成一个大问题并不是我的意思。我致力于在单个3D阵列操作中并行解决各个问题，并在第三维上使用不同的数据集。我已经尝试过这种'快速和肮脏的'，在我的情况下（200万数据集），它比循环遍历单个数据集要快500倍（！）。 –

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

回答

相关问题