2013-05-06 4 views
1

Я нашёл и скопировал этот код, чтобы получить FWHM от Finding the full width half maximum of a peak (от 2 до последнего ответа). Мой код ниже использует мои собственные данные. Сгенерированный сюжет выглядит неправильно, так как мои данные появляются в одной стороне, а зеленая коробка - с другой стороны. Что нужно изменить, чтобы я мог видеть гауссовские данные?Ошибка для моего участка и расчет FWHM

import numpy as np, scipy.optimize as opt 
from pylab import * 

def gauss(x, p): 
    return 1.0/(p[1]*np.sqrt(2*np.pi))*np.exp(-(x-p[0])**2/(2*p[1]**2)) 

x = [6711.19873047, 6712.74267578, 6714.28710938, 6715.83544922, \ 
    6717.38037109, 6718.92919922, 6720.47509766, 6722.02490234, \ 
    6723.57128906, 6725.11767578, 6726.66845703, 6728.21630859, \ 
    6729.76757812, 6731.31591797, 6732.86816406, 6734.41699219, \ 
    6735.96630859, 6737.51953125, 6739.06933594, 6740.62353516, \ 
    6742.17431641, 6743.72900391] 

y = [20.86093712, 23.60984612, 23.079916, 18.17703056, 18.24843597, \ 
    16.70049095, 19.48906136, 16.7509613, 19.09896088, 32.03007889, \ 
    54.56513977, 58.76417542, 40.93075562, 24.77710915, 17.68757629, \ 
    17.60736847, 18.89552498, 17.84486008, 17.49455452, 18.29696465, \ 
    18.55847931, 19.26465797] 

# Fit a guassian 
p0 = [0,70] 
errfunc = lambda p, x, y: gauss(x, p) - y # Distance to the target function 
p1, success = opt.leastsq(errfunc, p0[:], args=(x, y)) 

fit_mu, fit_stdev = p1 

FWHM = 2*np.sqrt(2*np.log(2))*fit_stdev 
print "FWHM", FWHM 

plot(x,y) 
plot(x, gauss(x,p1),lw=3,alpha=.5, color='r') 
axvspan(fit_mu-FWHM/2, fit_mu+FWHM/2, facecolor='g', alpha=0.5) 
show() 

ответ

2

Похоже первоначальной догадке, p0, вы передаете scipy.optimize.leastsq() проблема здесь. Если он недостаточно близко (например, не в пределах диапазона данных), то возвращаемое решение не будет иметь большого смысла.

Чтобы получить разумный выход, вы можете инициализировать p0 = [6730,70], что приведет к следующему сюжету:

A better initial guess for the expected distribution.

Но это все еще не правильный ответ, потому что вы пропустили шаг нормализации PDF присутствует в original answer. Вот нормализация сниппет, поскольку это относится к коду (добавьте его, прежде чем подходит гауссовой):

# Renormalize to a proper PDF 
y /= ((max(x) - min(x))/len(x)) * np.sum(y) 

Добавление это приводит к следующему сюжету:

Properly renormalized PDF

гауссовой подходит не очень но по крайней мере он имеет разумное среднее и стандартное отклонение, как я изначально описывал ниже.


Проблема связана с выходом функции наименьших квадратов, p1, которая является такой же, как p0.

Другими словами, учитывая, что полученное среднее значение равно 0 вместо np.mean(x) = 6727.45, блок будет нанесен на график в совершенно другом месте. Кроме того, ширина блока строится как 329.67 на основе стандартного отклонения 70, а не 46.28 на основе np.std(x) = 9.83.

+0

Эти данные должны иметь постоянное смещение в подгонке ..... – tacaswell

+0

@tcaswell: Конечно, я просто следил за связанным ответом. Не стесняйтесь изменять мой ответ. – fgb

Смежные вопросы