蟒蛇curve_fit沒有給出合理的擬合結果

Question

我想要適合高斯光譜和y值在10 ^（ - 19）的秩序。 Curve_fit在我將整個數據乘以10 ^（ - 19）之前和之後都給了我一個很差的擬合結果。附件是我的代碼，它是相當簡單的一組數據，除了值非常小。如果我想保留我的原始值，我將如何得到一個合理的高斯擬合，這會給我正確的參數？

#get fits data 
aaa=pyfits.getdata('p1.cal.fits') 

aaa=np.matrix(aaa) 
nrow=np.shape(aaa)[0] 
ncol=np.shape(aaa)[1] 

ylo=79 
yhi=90 
xlo=0 
xhi=1023 
glo=430 
ghi=470 

#sum all the rows to get spectrum 
ysum=[] 
for x in range(xlo,xhi): 
sum=np.sum(aaa[ylo:yhi,x]) 
ysum.append(sum) 

wavelen_pix=range(xhi-xlo) 
max=np.max(ysum) 
print "maximum is at x=", np.where(ysum==max) 

##fit gaussian 
#fit only part of my data in the chosen range [glo:ghi] 
x=wavelen_pix[glo:ghi] 
y=ysum[glo:ghi] 
def func(x, a, x0, sigma): 
    return a*np.exp(-(x-x0)**2/float((2*sigma**2))) 

sig=np.std(ysum[500:1000]) #std of background noise 

popt, pcov = curve_fit(func, x, sig) 
print popt 
#this gives me [1.,1.,1.], which is obviously wrong 
gaus=func(x,popt[0],popt[1],popt[2]) 

AAA是153×1024的圖像矩陣，部分看起來是這樣的：

matrix([[ -8.99793629e-20, 8.57133275e-21, 4.83523386e-20, ..., 
-1.54811004e-20, 5.22941515e-20, 1.71179195e-20], 
[ 2.75769318e-20, 1.03177243e-20, -3.19634928e-21, ..., 
1.66583803e-20, -9.88712568e-22, -2.56897725e-20], 
[ 2.88121935e-20, 8.57964252e-21, -2.60784327e-20, ..., 
1.72335180e-20, -7.61189937e-21, -3.45333075e-20], 
..., 
[ 1.04006903e-20, 1.61200683e-20, 7.04195205e-20, ..., 
1.72459645e-20, 4.29404029e-20, 1.99889374e-20], 
[ 3.22315752e-21, -5.61394194e-21, 3.28763096e-20, ..., 
1.99063583e-20, 2.12989880e-20, -1.23250648e-21], 
[ 3.66591810e-20, -8.08647455e-22, -6.22773168e-20, ..., 
-4.06145681e-21, 4.92453132e-21, 4.23689309e-20]], dtype=float32) 

Answer 1

您錯誤地調用curve_fit，這裏是使用

curve_fit(f, xdata, ydata, p0=None, sigma=None, absolute_sigma=False, check_finite=True, **kw) 

• ˚F是你的函數，它的第一個參數是一個獨立變量的數組，其後的參數是a重函數的參數（諸如幅度，中心等）
• XDATA是獨立變量
• YDATA是dependedent可變
• P0是在函數參數的初始猜測（爲高斯這是振幅，寬度，中心）

默認情況下，p0設置爲一個[1,1，...]的列表，這可能是您爲什麼得到這個結果，因爲你錯誤地稱呼它。

嘗試從數據估計振幅，中心和寬度，然後進行P0對象（詳見下文）

init_guess = (a_i, x0_i, sig_i) # same order as they are supplied to your function 
popt, pcov = curve_fit(func, xdata=x,ydata=y,p0=init_guess) 

下面是一個簡單的例子

xdata = np.linspace(0, 4, 50) 
mygauss = (10,2,0.5) #(amp, center, width) 
y  = func(xdata, *mygauss ) # using your func defined above  
ydata = y + 2*(np.random.random(50)- 0.5) # add some noise to create fake data 

現在我可以猜測擬合參數

ai = np.max(ydata) # guess the amplitude 
xi = xdata[ np.argmax(ydata)] # guess the position of center 

猜測寬度是棘手的，我會先找到半最大值是l ocated（有兩個，但你只需要找到一個，因爲高斯是對稱的）：

pos_half = argmin(np.abs(ydata-ao/2)) # subtract half the amplitude and find the minimum 

現在評估，這到底是高斯（XI）的中心：

sig_i = np.abs(xi - xdata[ pos_half]) # estimate the width 

現在您可以做出初步猜測

init_guess = (ai, xi sig_i) 

，適合

params, variance = curve_fit(func, xdata=xdata, ydata=ydata, p0=init_guess) 
print params 
#array([ 9.99457443, 2.01992858, 0.49599629]) 

這非常接近mygauss。希望能幫助到你。