Using lmfit model to best fit the multiple peaks data

Question

import matplotlib.pyplot as plt
import numpy as np
from lmfit import Model
x=np.array([4698.031, 4698.027, 4698.024, 4698.021, 4698.017, 4698.014,4698.011, 4698.007, 4698.004, 4698.001, 4697.997, 4697.994, 4697.991, 4697.987, 4697.984, 4697.981, 4697.977, 4697.974, 4697.971, 4697.967, 4697.964, 4697.961, 4697.957, 4697.954, 4697.951, 4697.947, 4697.944, 4697.941, 4697.937, 4697.934, 4697.931, 4697.927, 4697.924, 4697.921, 4697.917])
y=np.array([0.56565, 0.586575, 0.70335, 0.991245, 1.447545, 4.944375, 11.97281, 18.22095, 19.7613, 17.13792, 13.35083, 10.26506, 7.898505, 5.084775, 2.4192, 1.34358, 0.829905, 1.31322, 3.2049, 4.0095, 2.83263, 1.51605, 0.643275, 0.48972, 0.432675, 0.084375, 0.135345, 0.362145, 0.34425, 0.307125, 0.469125, 0.297, 0.183255, 0.528855, 0.523125])
gmodel = Model(gaussian, prefix='p1_') + Model(gaussian, prefix='p2_')
params = gmodel.make_params(p1_amp=0.1, p1_cen=4697.97, p1_wid=0.005, p2_amp=0.5, p2_cen=4698.00, p2_wid=0.005)


params['p1_cen'].min = x.min()
params['p1_cen'].max = 4697.98
params['p2_cen'].min = 4697.98
params['p2_cen'].max = x.max()

result = gmodel.fit(y, params, x=x)

updated program few questions

1. When I run the program it gives me an error. Traceback (most recent call last):

File "", line 7, in gmodel = Model(gaussian, prefix='p1_') + Model(gaussian, prefix='p2_')

NameError: name 'gaussian' is not defined

For which i include from lmfit.lineshapes import gaussianand it works.

Is is right or wrong??

2. After adding the library/package the program says Traceback (most recent call last):

File "", line 9, in params['p1_'].min = x.min()

KeyError: 'p1_cen'

help me with this.

3. After removing the params the program runs without any errors but what do i do in order to show plots for my given data n the best fit data . I just did plt.plot(x,y) for my given data but no idea what to do for my best fit data.....

Answer 1

You need to give better starting guesses. Your data x ranges from 4697.8 to 4698.0, but your initial guess for cen is 1.0, way outside your data range. Change the initial values to more sensible values for your data.

Update to address the question of how to fit 2 peaks:

To fit 2 peaks, you can create a composite model, as with

gmodel = Model(gaussian, prefix='p1_') + Model(gaussian, prefix='p2_')

result  = gmodel.fit(y, x=x, p1_amp=0.1, p1_cen=4697.97, p1_wid=0.005,
                     p2_amp=0.5, p2_cen=4698.00, p2_wid=0.005)

This will create a model that has 2 Gaussians added together, and giving the parameter names for one model a prefix p1_ and the parameters for the other a prefix of p2_. You can then give initial values for each of these when doing the fit.

It is generally recommended to create an instance of lmfit.Parameters for each model, as this allows you to further manipulate the parameter settings, for example, to set upper and lower bounds. For this kind of model, it might be helpful to make sure the peaks cannot overlap too closely or overlap. You could do that like this:

gmodel = Model(gaussian, prefix='p1_') + Model(gaussian, prefix='p2_')

params = gmodel.make_params(p1_amp=0.1, p1_cen=4697.97, p1_wid=0.005,
                            p2_amp=0.5, p2_cen=4698.00, p2_wid=0.005)
params['p1_cen'].min = x.min()
params['p1_cen'].max = 4697.98
params['p2_cen'].min = 4697.98
params['p2_cen'].max = x.max()

result = gmodel.fit(y, params, x=x)

which sets min/max values for the cen parameters for each model so that they cannot swap.