Different constraints for fit parameter in lmfit model

Question

I am trying to create a multible voigt/Gaussian/Lorentizan-peak fit function with lmfit. Therefore, I wrote the following Function:

def apply_fit_mix_multy(data,modelPeak,peakPos,amplitud,**kwargs):
peakPos=np.array(peakPos)
Start=kwargs.get('Start',data[0,0])
length_data=len(data)-1
End=kwargs.get('End',data[length_data,0])
StartPeak=kwargs.get('StartPeak',data[0,0])
EndPeak=kwargs.get('EndPeak',data[length_data,0])
BackFunc=kwargs.get('BackFunc',False)
BackCut=kwargs.get('BackCut',False)
dataN=data_intervall(data,Start,End)
y=dataN[:, 1]
x=dataN[:, 0]
amplitud=amplitud
center=peakPos

mod = None
for i in range(len(peakPos)):
    this_mod = make_model(i,amplitud,center,modelPeak)
    if mod is None:
        mod = this_mod
    else:
        mod = mod + this_mod

bgy=[list() for f in range(len(x))]
if(BackFunc==True):
    bg,bgx=BackFunc
    for i in range(len(x)):
        bgy[i]=bg.best_values.get('c')        

elif(BackCut!=False):
    slope,intercept=back_ground_cut(data,BackCut[0],BackCut[1])  
    for i in range(len(x)):
        bgy[i]=slope*x[i]+intercept       

if(BackCut!=False):
    print('Background substraction model is used! (Sign=Sign-backgr(linear between two points))')
    y=y-bgy
    out = mod.fit(y, x=x)
else:  
    print('Combination model is used! (offset+Gauss/Lor/Voig)')
    offset=ConstantModel()
    mod=mod+offset

out = mod.fit(y, x=x)#out is the fitted function

area=[list() for f in range(len(peakPos))]
comps=out.eval_components(x=x)
if(BackCut!=False):
    for i in range(len(peakPos)):
        area[i]=simps(comps['peak'+str(i)+'_'],x=x,even='avg')-simps(bgy,x=x,even='avg')
    fit_dict={'signal':y, 'convol':out.best_fit,'x':x,'peak_area':area,'backgr':bgy,'comps':comps}
else:
    for i in range(len(peakPos)):
        area[i]=simps(comps['peak'+str(i)+'_'],x=x,even='avg')
    fit_dict={'convol':out.best_fit,'x':x,'peak_area':area,'comps':comps} #comps is inf. of  sperate peaks

return fit_dict

The function reads in a data set, the modelPeak (e.g. GaussianModel) an initial guess of peak positions and amplitudes (peakPos, amplitude) .

In the first Part I initializing the model of the peaks (how many peaks...)

    for i in range(len(peakPos)):
    this_mod = make_model(i,amplitud,center,modelPeak)
    if mod is None:
        mod = this_mod
    else:
        mod = mod + this_mod

With the make_model funktion:

def make_model(num,amplitud,center,mod):
pref = "peak{0}_".format(num)
model = mod(prefix = pref)
model.set_param_hint(pref+'amplitud', value=amplitud[num], min=0, max=5*amplitud[num])
model.set_param_hint(pref+'center', value=center[num], min=center[num]-0.5, max=center[num]+0.5)
if(num==0):
    model.set_param_hint(pref+'sigma', value=0.3, min=0.01, max=1)
else:
    model.set_param_hint(pref+'sigma', value=0.3, min=0.01, max=1)
#print('Jetzt',center[num],amplitud[num])
return model

here is now my Problem: I I whant to fit e.g. 3 Peaks I whant that e.g. the sigma of the first peak is varies during the fit while the sigmas of the other peaks depend on the sigma of the first peak! any idea? thx maths

FYI this is how a fit looks like: enter image description here

Hi Matthias, as M Newville pointed out in his answer there is plenty of code unrelated to your question. You can easily remove stuff like back ground subtraction etc. Try to give a [minimum working example](https://stackoverflow.com/help/mcve) (e.g. what is `data_intervall()`?). This might include a small code part that generates some generic data. ...and please fix the indentation. — mikuszefski, May 16 '18 at 07:30

score 2 · Accepted Answer · answered May 08 '18 at 18:44

If I understand your long question (it would be helpful to remove the extraneous stuff - and there is quite a lot of it), you want to create a Model with multiple peaks, allowing sigma from the 1st peak to vary freely, and constraining sigma for the other peaks to depend on this.

To do that, you can either use parameter hints (as you use in your make_model() function) or set expressions for the parameters after the Parameters object is created. For the first approach, something like this

def make_model(num,amplitud,center,mod):
    pref = "peak{0}_".format(num)
    model = mod(prefix = pref)
    model.set_param_hint(pref+'amplitud', value=amplitud[num], min=0, max=5*amplitud[num])
    model.set_param_hint(pref+'center', value=center[num], min=center[num]-0.5, max=center[num]+0.5)
    if(num==0):
        model.set_param_hint(pref+'sigma', value=0.3, min=0.01, max=1)
    else:
        ## instead of 
        # model.set_param_hint(pref+'sigma', value=0.3, min=0.01, max=1)
        ## set peakN_sigma == peak0_sigma
        model.set_param_hint(pref+'sigma', expr='peak0_sigma') 
        ## or maybe set peakN_sigma == N * peak0_sigma 
        model.set_param_hint(pref+'sigma', expr='%d*peak0_sigma' % num)
    return model

You could also make the full model (simplified somewhat from your code, but the same idea):

model = (VoigtModel(prefix='peak0_') + VoigtModel(prefix='peak1_') +
         VoigtModel(prefix='peak2_') + LinearModel(prefix='const_'))

# create parameters with default values
params = model.make_params(peak0_amplitude=10, peak0_sigma=2, ....)

# set constraints for `sigma` params:
params['peak1_sigma'].expr = 'peak0_sigma'
params['peak2_sigma'].expr = 'peak0_sigma'

# similarly, set bounds as needed:
params['peak1_sigma'].min = 0
params['peak1_amplitude'].min = 0

Hope that helps...

Different constraints for fit parameter in lmfit model

1 Answers1

Linked