I'm writing a script that is going to be used to fit some photoluminescence spectra with custom models and after creating a script using SciPy I learned that setting bounds on the fitting parameters is much easier if lmfit (cars9.uchicago.edu/software/python/lmfit/intro.html) is used instead so I decided to use that instead. Below is a script I wrote to fit one spectra with one Gaussian function (using 2 Gaussian functions is better, but I wanted to start with a simpler case).
from numpy import loadtxt, vstack, average, exp
import matplotlib.pyplot as plt
import lmfit
def sub_BG(data_array):
naujas_X = data_array[:,1]-average(data_array[1:10,1])
return vstack((data_array[:,0], naujas_X)).T
def gauss(x, a1, b1, c1):
return a1 * exp(-(x- b1)**2/ (2*c1**2))
#amp * exp(-(x-cen)**2 /wid)
data="single_spectra.txt"
spectra = loadtxt(open(data, 'r'), skiprows=2)
spectra_beBG = sub_BG(spectra)
plt.plot(spectra_beBG[:,0],spectra_beBG[:,1],'g')
mod1 = lmfit.Model(gauss)
pars = lmfit.Parameters()
# (Name, Value, Vary, Min, Max, Expr)
pars.add_many(('a1', 590, True, None, None, None),
('b1', 500, True, None, None, None),
('c1', 20, True, None, None , None))
out = mod1.fit(spectra_beBG[:,0], pars, x=spectra_beBG[:,1])
y = gauss(spectra_beBG[:,0],
out.best_values["a1"],
out.best_values["b1"],
out.best_values["c1"])
plt.plot(spectra_beBG[:,0], out.best_fit, "r--")
plt.plot(spectra_beBG[:,0], y, "b--")
print(out.fit_report())
This returns:
[[Model]]
Model(gauss)
[[Fit Statistics]]
# function evals = 77
# data points = 1024
# variables = 3
chi-square = 28469283.530
reduced chi-square = 27883.725
[[Variables]]
a1: 561.593868 +/- 8.255604 (1.47%) (init= 590)
b1: 100.129107 +/- 85.34384 (85.23%) (init= 500)
c1: 1.3254e+06 +/- 7.23e+05 (54.52%) (init= 20)
[[Correlations]] (unreported correlations are < 0.100)
C(b1, c1) = -0.892
C(a1, c1) = -0.763
C(a1, b1) = 0.685
If I change pars.add_many() to something closer to nature, for example:
pars.add_many(('a1', 590, True, 550, 630, None),
('b1', 500, True, 450, 650, None),
('c1', 30, True, 20, 70 , None))
I get this:
[[Model]]
Model(gauss)
[[Fit Statistics]]
# function evals = 147
# data points = 1024
# variables = 3
chi-square = 304708538.428
reduced chi-square = 298441.272
[[Variables]]
a1: 629.999937 +/- 0 (0.00%) (init= 590)
b1: 475.821359 +/- 0 (0.00%) (init= 500)
c1: 70 +/- 0 (0.00%) (init= 30)
[[Correlations]] (unreported correlations are < 0.100)
Help?
Hm, should that be
out = mod1.fit(spectra_beBG[:,1], pars, x=spectra_beBG[:,0])
That is, you want to fit "y", and pass in the "pars" and "x" array to help calculate the model with those parameters and independent variables.