蟒蛇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。希望能帮助到你。