Я пытаюсь подогнать данные к двунаправленной экспоненциальной функции Гаусса (уравнение ниже) с помощью lmfit. Уравнение имеет четыре переменные: центр, амплитуду, сигму и гамма. Первые три всегда положительны, и я могу справиться с минимальными и максимальными ограничениями. Однако гамма может быть отрицательной (отражая фронтальный, а не хвостатый пик). Если я устанавливаю границы параметра всегда отрицательными или всегда положительными, у меня не возникает проблем с подгонкой, но возникают проблемы, когда мне нужно «пересечь» ноль, когда функция недействительна.
Функция
def model_n_expgaus(x,y,number_peaks,peak_locations,peak_heights,gamma_min=-2,gamma_max=2,gamma_center=0.1,peak_variation=0.1, sigma_min=0.005,sigma_max=1,sigma_center=0.1,height_scale=0.01,height_scale_high=0.8):
"""
Fits n expoential gaussians to chromatographic data and returns list of parameters and the area of each peak
Inputs:
x: Time axis, list
y: CAD, FLD axis, list
number_peaks: Number of peaks, int
peak_locations: locations of peaks along x axis, float
peak_heights: Value on y axis of each peak, float
gamma_min: float:
mimimum value of gamma for exponential gaussian fit (see lmfit documentation)
gamma_max: float:
maximum value of gamma for exponential gaussian fit (see lmfit documentation)
gamma_center: float:
initial guess of value of gamma for exponential gaussian fit (see lmfit documentation)
peak_variation: float
max variation of center for exponential gaussian fits
sigma_min: float
min sigma for exponential gaussian fits
sigma_max: float
max sigma for exponential gaussian fits
sigma_center: float
guess for sigma for exponential gaussian fits
height_scale: float
minimum amplitude for exponential gaussian fits
height_scale_high: float
maximum amplitude for exponential gaussian fits
Outputs: Array of tuples (parameters, area). Tuples are define each peak, and each row in the array is an individual peak
Prints(optional):
Summary graphs
1. Peak identification
2. Intial Fit
3. Final Fit+ individual chromatographs
"""
#goes from max to min number of peaks, finds bic, removes smallest peak, repeat
model=''
for n in range(0,number_peaks):
gaus=lmfit.Model(bemg,prefix='p'+str(n)+"_")
#addition of first model for first peak
if model=='':
model=gaus
pars=gaus.make_params(prefix='p'+str(n)+"_")
#peak locations can vary by 1 min
pars['p'+str(n)+'_center'].set(value=peak_locations[n],min=peak_locations[n]-peak_variation,max=peak_locations[n]+peak_variation)
#sigma can only be so large and small
pars['p'+str(n)+'_sigma'].set(value=sigma_center, min=sigma_min,max=sigma_max)
#amplitude can only be so small to be a real peak. Height correction from Amplitude to height is roughly /gamma~2
pars['p'+str(n)+'_amplitude'].set(value=peak_heights[n]*0.5, min=height_scale*peak_heights[n],max=peak_heights[n]*height_scale_high)
pars['p'+str(n)+'_gamma'].set(value=gamma_center,max=gamma_max,min=gamma_min)
#all subsequent peaks
else:
model+=gaus
pars.update(gaus.make_params(prefix='p'+str(n)+"_"))
pars['p'+str(n)+'_center'].set(value=peak_locations[n],min=peak_locations[n]-peak_variation,max=peak_locations[n]+peak_variation)
#sigma can only be so large and small
pars['p'+str(n)+'_sigma'].set(value=sigma_center, min=sigma_min,max=sigma_max)
#amplitude can only be so small to be a real peak Height correction from Amplitude to height is roughly /gamma~2
pars['p'+str(n)+'_amplitude'].set(value=peak_heights[n]*0.5, min=height_scale*peak_heights[n],max=peak_heights[n]*height_scale_high)
pars['p'+str(n)+'_gamma'].set(value=gamma_center,max=gamma_max,min=gamma_min)
#Uses gradient descent to best fit model
try:
out = model.fit(y, pars, x=x,nan_policy='omit')
print(out.fit_report())
except AttributeError:
return None
comps=out.eval_components(x=x)
#plots each peak
areas=[]
for n in range(0,number_peaks):
area=sum(comps['p'+str(n)+"_"])*0.01
center=out.params['p'+str(n)+"_center"].value
sigma=out.params['p'+str(n)+"_sigma"].value
amplitude=out.params['p'+str(n)+"_amplitude"].value
gamma=out.params['p'+str(n)+"_gamma"].value
areas.append([center,sigma,amplitude, gamma,area/6])
return(areas)
Я попробовал пару вещей.
Первым было попробовать установить lmfit nan_policy='omit'. Кажется, это не всегда работает, и я получаю следующую ошибку: ValueError: Массив, возвращаемый функцией, изменил размер между вызовами.
Я также пытался установить шаг грубого перебора на 1, чтобы попытаться избежать нулевой области, но это, похоже, тоже не работает.
Я пытаюсь подогнать данные к двунаправленной экспоненциальной функции Гаусса (уравнение ниже) с помощью lmfit. Уравнение имеет четыре переменные: центр, амплитуду, сигму и гамма. Первые три всегда положительны, и я могу справиться с минимальными и максимальными ограничениями. Однако гамма может быть отрицательной (отражая фронтальный, а не хвостатый пик). Если я устанавливаю границы параметра всегда отрицательными или всегда положительными, у меня не возникает проблем с подгонкой, но возникают проблемы, когда мне нужно «пересечь» ноль, когда функция недействительна. Функция [code]def bemg(x,amplitude,center,sigma,gamma): return np.exp(sigma*sigma/(2*gamma**2)+(center-x)/(gamma))*special.erfc((-1*math.copysign(1, gamma)*(x-center)/sigma-sigma/(gamma))/math.sqrt(2))*math.copysign(1, gamma)*amplitude/(2*gamma) [/code] Функция lmfit [code]def model_n_expgaus(x,y,number_peaks,peak_locations,peak_heights,gamma_min=-2,gamma_max=2,gamma_center=0.1,peak_variation=0.1, sigma_min=0.005,sigma_max=1,sigma_center=0.1,height_scale=0.01,height_scale_high=0.8): """ Fits n expoential gaussians to chromatographic data and returns list of parameters and the area of each peak Inputs: x: Time axis, list y: CAD, FLD axis, list
number_peaks: Number of peaks, int peak_locations: locations of peaks along x axis, float peak_heights: Value on y axis of each peak, float
gamma_min: float: mimimum value of gamma for exponential gaussian fit (see lmfit documentation) gamma_max: float: maximum value of gamma for exponential gaussian fit (see lmfit documentation) gamma_center: float: initial guess of value of gamma for exponential gaussian fit (see lmfit documentation) peak_variation: float max variation of center for exponential gaussian fits sigma_min: float min sigma for exponential gaussian fits sigma_max: float max sigma for exponential gaussian fits sigma_center: float guess for sigma for exponential gaussian fits height_scale: float minimum amplitude for exponential gaussian fits height_scale_high: float maximum amplitude for exponential gaussian fits Outputs: Array of tuples (parameters, area). Tuples are define each peak, and each row in the array is an individual peak
Prints(optional): Summary graphs 1. Peak identification 2. Intial Fit 3. Final Fit+ individual chromatographs """
#goes from max to min number of peaks, finds bic, removes smallest peak, repeat model='' for n in range(0,number_peaks): gaus=lmfit.Model(bemg,prefix='p'+str(n)+"_") #addition of first model for first peak if model=='': model=gaus pars=gaus.make_params(prefix='p'+str(n)+"_") #peak locations can vary by 1 min pars['p'+str(n)+'_center'].set(value=peak_locations[n],min=peak_locations[n]-peak_variation,max=peak_locations[n]+peak_variation) #sigma can only be so large and small pars['p'+str(n)+'_sigma'].set(value=sigma_center, min=sigma_min,max=sigma_max) #amplitude can only be so small to be a real peak. Height correction from Amplitude to height is roughly /gamma~2 pars['p'+str(n)+'_amplitude'].set(value=peak_heights[n]*0.5, min=height_scale*peak_heights[n],max=peak_heights[n]*height_scale_high) pars['p'+str(n)+'_gamma'].set(value=gamma_center,max=gamma_max,min=gamma_min) #all subsequent peaks else: model+=gaus pars.update(gaus.make_params(prefix='p'+str(n)+"_")) pars['p'+str(n)+'_center'].set(value=peak_locations[n],min=peak_locations[n]-peak_variation,max=peak_locations[n]+peak_variation) #sigma can only be so large and small pars['p'+str(n)+'_sigma'].set(value=sigma_center, min=sigma_min,max=sigma_max) #amplitude can only be so small to be a real peak Height correction from Amplitude to height is roughly /gamma~2 pars['p'+str(n)+'_amplitude'].set(value=peak_heights[n]*0.5, min=height_scale*peak_heights[n],max=peak_heights[n]*height_scale_high) pars['p'+str(n)+'_gamma'].set(value=gamma_center,max=gamma_max,min=gamma_min) #Uses gradient descent to best fit model try: out = model.fit(y, pars, x=x,nan_policy='omit') print(out.fit_report()) except AttributeError: return None
comps=out.eval_components(x=x)
#plots each peak areas=[] for n in range(0,number_peaks): area=sum(comps['p'+str(n)+"_"])*0.01 center=out.params['p'+str(n)+"_center"].value sigma=out.params['p'+str(n)+"_sigma"].value amplitude=out.params['p'+str(n)+"_amplitude"].value gamma=out.params['p'+str(n)+"_gamma"].value
return(areas) [/code] Я попробовал пару вещей. Первым было попробовать установить lmfit nan_policy='omit'. Кажется, это не всегда работает, и я получаю следующую ошибку: ValueError: Массив, возвращаемый функцией, изменил размер между вызовами. Я также пытался установить шаг грубого перебора на 1, чтобы попытаться избежать нулевой области, но это, похоже, тоже не работает.
Мобильная версия