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?
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)