需要返回協方差的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多項式擬合函數

回答

相關問題