Absolute Wertgrenzen für lmfit? Ich versuche, lmfit zur nichtlinearen Anpassung der kleinsten Quadrate an eine Funktion ⇐ Python

[Anonymous]

Ich versuche, Daten mithilfe von lmfit an eine bidirektionale exponentielle Gaußsche Funktion (Gleichung unten) anzupassen. Die Gleichung hat vier Variablen: Zentrum, Amplitude, Sigma und Gamma. Die ersten drei sind immer positiv und ich komme mit Min/Max-Grenzwerten zurecht. Allerdings kann Gamma negativ sein (was eher einen vorderen als einen nachlaufenden Peak widerspiegelt). Wenn ich die Parametergrenzen auf immer negativ oder immer positiv einstelle, habe ich keine Probleme mit der Anpassung, aber es gibt Probleme, wenn ich Null „kreuzen“ muss, wenn die Funktion ungültig ist.
Funktion

Code: Select all

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-Funktion

Code: Select all

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)

Ich habe ein paar Dinge ausprobiert.
Das erste war, zu versuchen, lmfit nan_policy='omit' zu setzen. Dies scheint nicht immer zu funktionieren und ich erhalte diese Fehlermeldung: ValueError: Das von einer Funktion zurückgegebene Array hat seine Größe zwischen den Aufrufen geändert
Ich habe auch versucht, den Brute-Schritt auf 1 zu setzen, um den Nullbereich zu vermeiden, aber das scheint auch nicht zu funktionieren.