python curve_fit 没有给出合理的拟合结果答案

【问题标题】：python curve_fit does not give reasonable fitting resultpython curve_fit 没有给出合理的拟合结果
【发布时间】：2015-07-22 03:18:53
【问题描述】：

我正在尝试将高斯拟合到光谱中，并且 y 值大约为 10^(-19)。在我将整个数据乘以 10^(-19) 之前和之后，Curve_fit 给了我很差的拟合结果。附件是我的代码，它是一组相当简单的数据，只是值非常小。如果我想保留我的原始值，我将如何获得一个合理的高斯拟合来给我正确的参数？

#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 x 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)

【问题讨论】：

您是否尝试使用一些合理大小的 y 值来测试您的代码？
如果您打印出aaa 的一部分，您的帖子也会受益匪浅..
刚才我尝试将所有内容乘以 10^(-19)，但结果仍然相同。所以这不是一个精确的事情。我会修改我的问题。谢谢！
注意：不需要 yn，它可能来自添加噪声并将噪声数据拟合到原始数据 (y) 的教程。

标签： python gaussian astronomy

【解决方案1】：

你打错了curve_fit，这里是用法

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

f 是你的函数，它的第一个 arg 是一个自变量数组，后面的 args 是函数参数（例如振幅、中心等）
xdata 是自变量
ydata 是因变量
p0 是对函数参数的初始猜测（对于 Guassian，这是幅度、宽度、中心）

默认情况下，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

猜测宽度很棘手，我会先找到半最大值所在的位置（有两个，但你只需要找到一个，因为高斯是对称的）：

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。希望对您有所帮助。

【讨论】：

【解决方案2】：

忘记重新缩放、进行线性更改或使用 p0 参数，这些通常不起作用！尝试在 curve_fit 中为 n 个参数使用 bounds 参数，如下所示：

a0=np.array([a01,...,a0n])
af=np.array([af1,...,afn])
method="trf",bounds=(a0,af)

希望它有效！ ;)

【讨论】：