2015-07-22 8 views
2

Я пытаюсь свести гауссовский к спектру и значения y составляют порядка 10^(- 19). Curve_fit дает мне плохой результат, как до, так и после того, как я умножу все свои данные на 10^(- 19). Прикрепленный мой код, это довольно простой набор данных, за исключением того, что значения очень малы. Если я хочу сохранить свои первоначальные значения, как мне получить разумную гауссовскую посадку, которая даст мне правильные параметры?python 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]) 

ааа является 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) 
+0

Вы пытались протестировать свой код с некоторыми разумными размерами y? – dermen

+0

также ваше сообщение принесет большую пользу, если вы напечатаете часть того, что '' 'aaa''' выглядит .. – dermen

+0

Сейчас я попытался размножить все на 10^(- 19), и он все равно дает мне такой же результат. Так что это не точная вещь. Я пересмотрю свой вопрос. Спасибо! – RaynDrop

ответ

1

Вы звоните curve_fit неправильно, здесь является использование

curve_fit(f, xdata, ydata, p0=None, sigma=None, absolute_sigma=False, check_finite=True, **kw) 
  • е является вашей функцией, первый аргумент которой представляет собой массив независимых переменных и последующие аргументы повторно параметры функции (например, амплитуды, в центре, и т.д.)
  • XData являются независимым переменными
  • ydata является dependedent переменного
  • р0 является начальным приближением при параметрах функции (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 

Теперь я могу Угадай подходят PARAMS

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

Гадание ширина сложно, я бы сначала найти, где половина МАХ л ocated (есть два, но вам нужно только найти, как Gaussian симметрична):

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. Надеюсь, поможет.

+0

Спасибо большое! Он отлично работает! @dermen – RaynDrop