Fitting Voigt Function To Data In Python

I recently got a script running to fit a gaussian to my absorption profile with help of SO. My hope was that things would work fine if I simply replace the Gauss function by a Voig

Solution 1:

I may be misunderstanding the model you're using, but I think you would need to include some sort of constant or linear background.

To do that with lmfit (which has Voigt, Gaussian, and many other models built in, and tries very hard to make these interchangeable), I would suggest starting with something like this:

import numpy as np
import matplotlib.pyplot as plt
from lmfit.models import GaussianModel, VoigtModel, LinearModel, ConstantModel

x = np.arange(13)
xx = np.linspace(0, 13, 100)
y = np.array([19699.959 , 21679.445 , 21143.195 , 20602.875 , 16246.769 ,
              11635.25  ,  8602.465 ,  7035.493 ,  6697.0337,  6510.092 ,
              7717.772 , 12270.446 , 16807.81  ])

# build model as Voigt + Constant## model = GaussianModel() + ConstantModel()
model = VoigtModel() + ConstantModel()

# create parameters with initial values
params = model.make_params(amplitude=-1e5, center=8, 
                           sigma=2, gamma=2, c=25000)

# maybe place bounds on some parameters
params['center'].min = 2
params['center'].max = 12
params['amplitude'].max = 0.# do the fit, print out report with results 
result = model.fit(y, params, x=x)
print(result.fit_report())

# plot data, best fit, fit interpolated to `xx`
plt.plot(x, y, 'b+:', label='data')
plt.plot(x, result.best_fit, 'ko', label='fitted points')
plt.plot(xx, result.eval(x=xx), 'r-', label='interpolated fit')
plt.legend()
plt.show()

And, yes, you can simply replace VoigtModel() with GaussianModel() or LorentzianModel() and redo the fit and compare the fit statistics to see which model is better.

For the Voigt model fit, the printed report would be

[[Model]]
    (Model(voigt) + Model(constant))
[[Fit Statistics]]
    # fitting method   = leastsq
    # functionevals   = 41
    # datapoints      = 13
    # variables        = 4
    chi-square         = 17548672.8
    reducedchi-square = 1949852.54
    Akaikeinfocrit   = 191.502014
    Bayesianinfocrit = 193.761811
[[Variables]]
    amplitude: -173004.338 +/- 30031.4068 (17.36%) (init = -100000)
    center:     8.06574198 +/- 0.16209266 (2.01%) (init = 8)
    sigma:      1.96247322 +/- 0.23522096 (11.99%) (init = 2)
    c:          23800.6655 +/- 1474.58991 (6.20%) (init = 25000)
    gamma:      1.96247322 +/- 0.23522096 (11.99%) == 'sigma'
    fwhm:       7.06743644 +/- 0.51511574 (7.29%) == '1.0692*gamma+sqrt(0.8664*gamma**2+5.545083*sigma**2)'
    height:    -18399.0337 +/- 2273.61672 (12.36%) == '(amplitude/(max(2.220446049250313e-16, sigma*sqrt(2*pi))))*wofz((1j*gamma)/(max(2.220446049250313e-16, sigma*sqrt(2)))).real'[[Correlations]] (unreported correlations are < 0.100)
    C(amplitude, c)     = -0.957
    C(amplitude, sigma) = -0.916
    C(sigma, c)         =  0.831
    C(center, c)        = -0.151

Note that by default gamma is constrained to be the same value as sigma. This constraint can be lifted and gamma made to vary independently with params['gamma'].set(expr=None, vary=True, min=1.e-9). I think that you may not have enough data points in this data set to robustly and independently determine gamma.

The plot for that fit would look like this:

Solution 2:

I managed to get something, but not very satisfying. If you remove the offset as a parameter and add 20000 directly in the Voigt function, with starting values [8, 126000, 0.71, 2] (the particular values don't' matter much) you'll get something like

Now, the fit produces a value for gamma which is negative which I cannot really justify. I would expect gamma to always be positive, but maybe I'm wrong and it's completely fine.

One thing you could try is to mirror your data so that its a "positive" peak (and while at it removing the background) and/or normalize the values. That might help you in the convergence.

I have no idea why when using the offset as a parameter the solver has problems finding an optimum. Maybe you need a different optimizer routine.

Maybe it'll be a better option to use the lmfit package that it's a wrapper over scipy to fit nonlinear functions with many prebuilt lineshapes. There is even an example of fitting a Voigt profile.