Search code examples
pythonmodel-fittinglmfit

Creating more elaborate composite models with lmfit


I want to create more elaborate / complex models from simple ones with lmfit.

I have two functions like e.g. a gaussian (norm. to unity at peak) and a lorentzian (norm. to unity at peak) and want to fit e.g. a linear combination of them. Still adding to unity. So I could write a new function like

def voigt(*all the parameters*, alpha)
   return alpha*gaussian(...) + (1-alpha)*lorentzian(...)

but this is not very adaptive. So I now do instead:

mod = ConstantModel(prefix = 'a1_') + Model(gauss) + ConstantModel(prefix='a2_') * Model(lorentz)
pars = Parameters()
pars.add('a1_c', value = 1, min = 0, max = 1)
pars.add('a2_c', expr = '1-a1_c')

Still feels a bit clumsy. Is there a more elegant way?


Solution

  • I think it mostly depends on what sort of complexity and adaptability you're looking for. If I understand the question (perhaps an example or simplification of your real goal), you want to constrain two peak functions to have maximal heights that sum to a fixed value (say, 1). Why not define basis functions that scale as you want, perhaps

     import numpy as np
     def gauss(x, amp, cen, wid):
          return amp * np.exp(-(x-cen)**2/wid)  # deliberately Gaussian-like
    
     def loren(x, amp, cen, wid): 
          return amp * wid / ( (x-cen)**2 + wid) # deliberately Lorentzian-like
    

    and then sum those and constrain the two amp parameters to sum to 1:

    from lmfit import Model
    mod = Model(gauss, prefix='g_') + Model(loren, prefix='l_')
    
    pars = mod.make_params(g_cen=1, g_wid=1, l_cen=1, l_wid=1)
    pars.add('g_amp', value=0.5, min=0, max=1)
    pars.add('l_amp', expr='1 - g_amp', min=0, max=1)