How to do a fit with a Split Pearson 7 function in Python?

I have an asymmetric peak that I want to fit using a split Pearson 7 (SP7) function with lmfit. First, the SP7 function has been defined, then I re-called it with lmfit, but at the end I have this error: ValueError: varargs '*p' is not supported. Why?

def polynomial(x, *p):
    poly = Polynomial(p)
    y = poly(x)
    return y

def pearson7(x, *p):
    x0, IM, H, m = p
    y = IM / (1+4*(2**(1/m)-1)*(x-x0)**2/H**2)**m
    return y

def splitpearson7(x, *p):
    x0, IM, Hl, ml, Hr, mr = p
    pl = [x0, IM, Hl, ml]
    pr = [x0, IM, Hr, mr]
    y = pearson7(x, *pl)*np.heaviside(x-x0, 0.5) + pearson7(x, *pr)*np.heaviside(x0-x, 0.5)
    return y

def splitpearson7_backg(x, *p):
    """
    Split PearsonVII distribution with linear background (8 parameters)
    * background y = A + Bx
    """
    y = splitpearson7(x, *p[0:6]) + polynomial(x, *p[6:8])
    return y ```

first_params_fit=[0, np.max(y), 0.05, 0.5, 0.05, 0.5, 10, 0.3]
mod = Model(splitpearson7_backg)
pars = mod.make_params(first_params_fit)
result = mod.fit(y, pars, x=x)
print(result.fit_report()) ```

Solution

ValueError: varargs '*p' is not supported.

Indeed.

Why?

Well, partly because it is a horrible way to define a function, and partly because it breaks the concept of lmfit.Model to be able to determine and name parameters for your model function. In short: readability counts. Do the person reading your code (including your future self) a favor and write readable code.

With your splitpearson7 the meaning and order of the parameters for the function are lost. The order is still critical of course, but it is not described in any way. It does not even say how many parameters there are. If, instead, you used

def splitpearson7(x, x0, IM, Hl, ml, Hr, mr):
    pl = [x0, IM, Hl, ml]
    pr = [x0, IM, Hr, mr]
    y = pearson7(x, *pl)*np.heaviside(x-x0, 0.5) + pearson7(x, *pr)*np.heaviside(x0-x, 0.5)
    return y

then you could create a lmfit.Model with

 sp7_model = lmfit.Model(splitpearson7)

and have named parameters:

 params = sp7_model.make_params(x0=0, IM=1, Hl=2, ml=3.3, Hr=2.1, mr=5)

Indeed, your pearson7 function could also use named parameters. Or, you could use the version builtin to lmfit:

from lmfit.lineshapes import pearson7
def splitpearson7(x, center, amplitude, sigma_left, m_left, sigma_right, m_right):
      left = heaviside(x-center, 0.5)*pearson7(x, amplitude, center, sigma_left, m_left) 
      right = heaviside(center-x, 0.5)*pearson7(x, amplitude, center, sigma_right, m_right) 
      return left+right

In addition, it is generally "the lmfit Model way" to not build a model that adds a function with a background, but to make a composite model, as with

from numpy import heaviside
from lmfit.lineshapes import pearson7
from lmfit.model import Model
from lmfit.models import PolynomialModel

def splitpearson7(x, center, amplitude, sigma_left, m_left, sigma_right, m_right):
    left = heaviside(x-center, 0.5)*pearson7(x, amplitude, center, sigma_left, m_left) 
    right = heaviside(center-x, 0.5)*pearson7(x, amplitude, center, sigma_right, m_right) 
    return left+right

# Note that PolynomialModel(2) will be quadratic with parameters 
# c0, c1, c2, and a form = c0 + c1*x + c2*x*x
mod = Model(splitpearson7) + PolynomialModel(2)
params = mod.make_params(center=0, amplitude=200, 
                         sigma_left=2, m_left=1, 
                         sigma_right=2.2, m_right=1.1, 
                         c0=10, c1=0, c2=0)

That will make it much easier to alter the model for the background.