Границы абсолютного значения для lmfit? Попытка использовать lmfit для нелинейного метода наименьших квадратов, подходящ ⇐ Python

[Anonymous]

Я пытаюсь подогнать данные к двунаправленной экспоненциальной функции Гаусса (уравнение ниже) с помощью lmfit. Уравнение имеет четыре переменные: центр, амплитуду, сигму и гамма. Первые три всегда положительны, и я могу справиться с минимальными и максимальными ограничениями. Однако гамма может быть отрицательной (отражая фронтальный, а не хвостатый пик). Если я устанавливаю границы параметра всегда отрицательными или всегда положительными, у меня не возникает проблем с подгонкой, но возникают проблемы, когда мне нужно «пересечь» ноль, когда функция недействительна.
Функция

Код: Выделить всё

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)

Функция 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, чтобы попытаться избежать нулевой области, но это, похоже, тоже не работает.

Подробнее здесь: https://stackoverflow.com/questions/795 ... st-squares