1

I have a function Imaginary which describes a physics process and I want to fit this to a dataset x_interpolate, y_interpolate. The function is a form of a Lorentzian peak function and I have some initial values that are user given, except for f_peak (the peak location) which I find using a peak finding algorithm. All of the fit parameters, except for the offset, are expected to be positive and thus I have set bounds_I accordingly. 

def Imaginary(freq, alpha, res, Ms, off):
    numerator = (2*alpha*freq*res**2)
    denominator = (4*(alpha*res*freq)**2) + (res**2 - freq**2)**2
    Im = Ms*(numerator/denominator) + off
    return Im

pI = np.array([alpha_init, f_peak, Ms_init, 0])

bounds_I = ([0,0,0,0, -np.inf], [np.inf,np.inf,np.inf, np.inf])

poptI, pcovI = curve_fit(Imaginary, x_interpolate, y_interpolate, pI, bounds=bounds_I)

In some situations I want to keep the parameter f_peak fixed during the fitting process. I tried an easy solution by changing bounds_I to:

bounds_I = ([0,f_peak+0.001,0,0, -np.inf], [np.inf,f_peak-0.001,np.inf, np.inf])

This is for many reasons not an optimal way of doing this so I was wondering if there is a more Pythonic way of doing this? Thank you for your help

STJ
  • 113
  • 1
  • 2
  • 12

3 Answers3

4

If a parameter is fixed, it is not really a parameter, so it should be removed from the list of parameters. Define a model that has that parameter replaced by a fixed value, and fit that. Example below, simplified for brevity and to be self-contained: 

x = np.arange(10)
y = np.sqrt(x)    
def parabola(x, a, b, c):
  return a*x**2 + b*x + c

fit1 = curve_fit(parabola, x, y)  #  [-0.02989396,  0.56204598,  0.25337086]
b_fixed = 0.5
fit2 = curve_fit(lambda x, a, c: parabola(x, a, b_fixed, c), x, y)

The second call to fit returns [-0.02350478, 0.35048631], which are the optimal values of a and c. The value of b was fixed at 0.5.

Of course, the parameter should be removed from the initial vector pI and the bounds as well.

0

You might find lmfit (https://lmfit.github.io/lmfit-py/) helpful. This library adds a higher-level interface to the scipy optimization routines, aiming for a more Pythonic approach to optimization and curve fitting. For example, it uses Parameter objects to allow setting bounds and fixing parameters without having to modify the objective or model function. For curve-fitting, it defines high level Model functions that can be used.

For you example, you could use your Imaginary function as you've written it with

from lmfit import Model
lmodel = Model(Imaginary)

and then create Parameters (lmfit will name the Parameter objects according to your function signature), providing initial values:

params = lmodel.make_params(alpha=alpha_init, res=f_peak, Ms=Ms_init, off=0)

By default all Parameters are unbound and will vary in the fit, but you can modify these attributes (without rewriting the model function):

params['alpha'].min = 0
params['res'].min = 0
params['Ms'].min = 0

You can set one (or more) of the parameters to not vary in the fit as with:

params['res'].vary = False

To be clear: this does not require altering the model function, making it much easier to change with is fixed, what bounds might be imposed, and so forth.

You would then perform the fit with the model and these parameters:

result = lmodel.fit(y_interpolate, params, freq=x_interpolate)

you can get a report of fit statistics, best-fit values and uncertainties for parameters with

print(result.fit_report())

The best fit Parameters will be held in result.params.

FWIW, lmfit also has builtin Models for many common forms, including Lorentzian and a Constant offset. So, you could construct this model as

from lmfit.models import LorentzianModel, ConstantModel

mymodel = LorentzianModel(prefix='l_') + ConstantModel()

params = mymodel.make_params()

which will have Parameters named l_amplitude, l_center, l_sigma, and c (where c is the constant) and the model will use the name x for the independent variable (your freq). This approach can become very convenient when you may want to change the functional form of the peaks or background, or when fitting multiple peaks to a spectrum.

M Newville
  • 7,486
  • 2
  • 16
  • 29
0

I was able to solve this issue regarding arbitrary number of parameters and arbitrary positioning of the fixed parameters:

def d_fit(x, y, param, boundMi, boundMx, listparam):
    Sparam, SboundMi, SboundMx = asarray([]), asarray([]), asarray([])
    Nparam, NboundMi, NboundMx = asarray([]), asarray([]), asarray([])
    for i in range(len(param)):
        if(listparam[i] == 1):
            Sparam = append(Sparam,asarray(param[i]))
            SboundMi = append(SboundMi,asarray(boundMi[i]))
            SboundMx = append(SboundMx,asarray(boundMx[i]))
        else:
            Nparam = append(Nparam,asarray(param[i]))
    def funF(x, Sparam):
        j = 0
        for i in range(len(param)):
            if(listparam[i] == 1):
                param[i] = Sparam[i-j]
            else:
                param[i] = Nparam[j]
                j = j + 1
        return fun(x, param)
    return curve_fit(lambda x, *Sparam: funF(x, Sparam), x, y, p0 = Sparam, bounds = (SboundMi,SboundMx))

In this case:

param = [a,b,c,...] # parameters array (any size)
boundMi = [min_a, min_b, min_c,...] # minimum allowable value of each parameter
boundMx = [max_a, max_b, max_c,...] # maximum allowable value of each parameter
listparam = [0,1,1,0,...] # 1 = fit and 0 = fix the corresponding parameter in the fit routine

and the root function is define as

def fun(x, param):
    a,b,c,d.... = param
    return a*b/c... # any function of the params a,b,c,d...

This way, you can change the root function and the number of parameters without changing the fit routine. And, at any time, you can fix or let fit any parameter by changing "listparam".

Use like this:

popt, pcov = d_fit(x, y, param, boundMi, boundMx, listparam)

"popt" and "pcov" are 1D arrays of the size of the number of "1" in "listparam" bringing the results of the fitted parameters (best value and err matrix)

"param" will ramain an 1D array of the same size of the original (input) "param", HOWEVER IT WILL BE UPDATED AUTOMATICALLY TO THE FITTED VALUES (same as "popt") for the fitted values, keeping the fixed values according to "listparam"

Hope can be usefull!

Obs1: x = 1D-array independent values and y = 1D-array dependent values

Obs2: This is my first post. Please let me know if I can improove it!

Douglas Leite
  • 1
  • 1