2011-08-19 77 views
7

我一直在用Python進行一些蒙特卡洛物理模擬,我無法確定非線性最小二乘擬合係數的標準誤差。非線性迴歸中的標準誤差

最初,我使用SciPy的scipy.stats.linregress作爲我的模型,因爲我認爲它將是一個線性模型,但注意到它實際上是某種功能函數。然後我使用NumPy的polyfit,自由度爲2,但我無法找到確定係數的標準誤差。

我知道gnuplot可以確定我的錯誤,但我需要做適合超過30種不同的情況。我想知道是否有人知道Python從gnuplot讀取標準錯誤還是有一些其他庫可以使用?

回答

2

它看起來像gnuplot的使用levenberg - marquardt並有一個Python implementation available - 你可以從mpfit.covar屬性(順便說一下,你應該擔心什麼錯誤估計「的意思是」錯誤的估計 - 是其他參數允許調整補償,例如?)

+0

感謝您的鏈接!最後我沒有使用mpfit,但是文檔確實幫助我理解了scipy的curve_fit! – syntaxing

10

終於找到了這個長問問題的答案!我希望這至少可以爲一個人節省幾個小時對這個話題的無望研究。在其優化部分,Scipy有一個稱爲curve_fit的特殊功能。它使用最小二乘法確定係數,最好的是,它給出了協方差矩陣。協方差矩陣包含每個係數的方差。更準確地說,矩陣的對角線是方差並且通過平方根值可以確定每個係數的標準誤差! SciPy的沒有太多的文檔這所以這裏是爲了更好地理解示例代碼:

import numpy as np 
from scipy.optimize import curve_fit 
import matplotlib.pyplot as plot 


def func(x,a,b,c): 
    return a*x**2 + b*x + C#Refer [1] 

x = np.linspace(0,4,50) 
y = func(x,2.6,2,3) + 4*np.random.normal(size=len(x)) #Refer [2] 


coeff, var_matrix = curve_fit(func,x,y) 
variance = np.diagonal(var_matrix) #Refer [3] 

SE = np.sqrt(variance) #Refer [4] 

#======Making a dictionary to print results======== 
results = {'a':[coeff[0],SE[0]],'b':[coeff[1],SE[1]],'c':[coeff[2],SE[2]]} 

print "Coeff\tValue\t\tError" 
for v,c in results.iteritems(): 
    print v,"\t",c[0],"\t",c[1] 
#========End Results Printing================= 

y2 = func(x,coeff[0],coeff[1],coeff[2]) #Saves the y values for the fitted model 

plot.plot(x,y) 
plot.plot(x,y2) 

plot.show() 
  1. 什麼這個函數返回是至關重要的,因爲它定義將用於擬合的模型
  2. 使用什麼函數來創建某些任意數據+一些噪聲
  3. 保存的協方差矩陣的對角線至1D矩陣這僅僅是一個普通陣列
  4. 廣場生根的方差得到的標準誤差(SE)