如何使用lmfit中的优化参数查找函数(伪Voigt)的区域? [英] How to find the area of a function (Pseudo Voigt) using optimized parameters from lmfit?
问题描述
我正在尝试确定曲线的面积(峰值).我能够使用伪Voigt配置文件和指数背景成功拟合峰(数据)，并获得与使用商业软件获得的参数相符的拟合参数.现在的问题是试图将那些拟合的峰参数与峰面积相关联.
与高斯线形的情况不同，我找不到使用拟合参数计算峰面积的简单方法.因此，我尝试使用scipy quad函数来集成我的拟合函数.我知道该区域应该由商业软件确定为19,000左右，但是我得到的数值非常大.
拟合效果很好(通过绘图确认...)，但是计算区域不紧密.尝试绘制具有最佳拟合值的psuedo_voigt_func函数后，我发现它的峰值太强了.这样，积分可能是正确的，那么错误将出在我如何创建峰的方法中，即如何将拟合的参数传递给我的psuedo_voigt_func函数，在该函数中，该函数是从lmfit模型网站进行转录的(
[3.52039054e + 01 8.13244777e-02 7.80206967e-02 1.96178358e + 04 4.48314849e + 03](34264.98814344757，0.00017531957481189617)(34241.971907301166，0.0002019796740206914)(18999.266974139795，0.0002019796990069267) 从图中可以明显看出 pseudo_voigt 不是很好.但是，积分相差不大.考虑到拟合优化 chi ** 2 的事实，这并不是什么大惊喜.
I am trying to determine the area of a curve (peak). I was able to successfully fit the peak (data) using a Pseudo Voigt profile and a exponential background and get the fitting parameters out which agree with parameters obtained using commercial software. The issue is now trying to relate those fitted peak parameters to the area of the peak.
I could not find an simple method of using the fitted parameters to calculate the area of the peak unlike in the case with a Gaussian line shape. So I am trying to use the scipy quad function to integrate my fitted function. I know the area should be around 19,000 determined by commercial software but I am getting very large incorrect values.
The fitting works very well (confirmed by plotting...) but the calculate area is not close. After trying to plot the my psuedo_voigt_func function with the best fit values passed to it, I found that it was a much too intense peak. With that, the integration might have been correct then the error would be in how I create my peak which was by passing the fitted parameters to my psuedo_voigt_func function, in which the function was transcribed from the lmfit model website (https://lmfit.github.io/lmfit-py/builtin_models.html). I believe I scripted the psuedo voigt function correctly but it isn't working.
#modules
import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from lmfit.models import GaussianModel, LinearModel, VoigtModel, Pearson7Model, ExponentialModel, PseudoVoigtModel
from scipy.integrate import quad
#data
x = np.array([33.05, 33.1 , 33.15, 33.2 , 33.25, 33.3 , 33.35, 33.4 , 33.45, 33.5 , 33.55, 33.6 , 33.65, 33.7 , 33.75, 33.8 , 33.85, 33.9 , 33.95, 34. , 34.05, 34.1 , 34.15, 34.2 , 34.25, 34.3 , 34.35, 34.4 , 34.45, 34.5 , 34.55, 34.6 , 34.65, 34.7 , 34.75, 34.8 , 34.85, 34.9 , 34.95, 35. , 35.05, 35.1 , 35.15, 35.2 , 35.25, 35.3 , 35.35, 35.4 , 35.45, 35.5 , 35.55, 35.6 , 35.65, 35.7 , 35.75, 35.8 , 35.85, 35.9 , 35.95, 36. , 36.05, 36.1 , 36.15, 36.2 , 36.25, 36.3 , 36.35, 36.4 , 36.45])
y = np.array([4569, 4736, 4610, 4563, 4639, 4574, 4619, 4473, 4488, 4512, 4474, 4640, 4691, 4621, 4671, 4749, 4657, 4751, 4921, 5003, 5071, 5041, 5121, 5165, 5352, 5304, 5408, 5393, 5544, 5625, 5859, 5851, 6155, 6647, 7150, 7809, 9017, 10967, 14122, 19529, 28029, 39535, 50684, 55730, 52525, 41356, 30015, 20345, 14368, 10736, 9012, 7765, 7064, 6336, 6011, 5806, 5461, 5283, 5224, 5221, 4895, 4980, 4895, 4852, 4889, 4821, 4872, 4802, 4928])
#model
bkg_model = ExponentialModel(prefix='bkg_') #BACKGROUND model
peak_model = PseudoVoigtModel(prefix='peak_') #PEAK model
model = peak_model + bkg_model
#parameters
pars = bkg_model.guess(y, x=x) #BACKGROUND parameters
pars.update(peak_model.make_params()) #PEAK parameters
pars['peak_amplitude'].set(value=17791.293, min=0)
pars['peak_center'].set(value=35.2, min=0, max=91)
pars['peak_sigma'].set(value=0.05, min=0)
#fitting
init = model.eval(pars, x=x) #initial parameters
out = model.fit(y, pars, x=x) #fitting
#integration part
def psuedo_voigt_func(x, amp, cen, sig, alpha):
sig_gauss = sig / np.sqrt(2*np.log(2))
term1 = (amp * (1-alpha)) / (sig_gauss * np.sqrt(2*np.pi))
term2 = np.exp(-(x-cen)**2) / (2 * sig_gauss**2)
term3 = ((amp*alpha) / np.pi) * ( sig / ((x-cen)**2) + sig**2)
psuedo_voigt = (term1 * term2) + term3
return psuedo_voigt
fitted_amp = out.best_values['peak_amplitude']
fitted_cen = out.best_values['peak_center']
fitted_sig = out.best_values['peak_sigma']
fitted_alpha = out.best_values['peak_fraction']
print(quad(psuedo_voigt_func, min(x), max(x), args=(fitted_amp, fitted_cen, fitted_sig, fitted_alpha)))
#output result of fit:
[[Model]]
(Model(pvoigt, prefix='peak_') + Model(exponential, prefix='bkg_'))
[[Fit Statistics]]
# fitting method = leastsq
# function evals = 542
# data points = 69
# variables = 6
chi-square = 2334800.79
reduced chi-square = 37060.3300
Akaike info crit = 731.623813
Bayesian info crit = 745.028452
[[Variables]]
bkg_amplitude: 4439.34760 +/- 819.477320 (18.46%) (init = 10.30444)
bkg_decay: -38229822.8 +/- 5.6275e+12 (14720258.94%) (init = -5.314193)
peak_amplitude: 19868.0711 +/- 106.363477 (0.54%) (init = 17791.29)
peak_center: 35.2039076 +/- 3.3971e-04 (0.00%) (init = 35.2)
peak_sigma: 0.14358871 +/- 5.6049e-04 (0.39%) (init = 0.05)
peak_fraction: 0.62733180 +/- 0.01233108 (1.97%) (init = 0.5)
peak_fwhm: 0.28717742 +/- 0.00112155 (0.39%) == '2.0000000*peak_sigma'
peak_height: 51851.3174 +/- 141.903066 (0.27%) == '(((1 peak_fraction)*peak_amplitude)/max(2.220446049250313e-16, (peak_sigma*sqrt(pi/log(2))))+(peak_fraction*peak_amplitude /max(2.220446049250313e-16, (pi*peak_sigma)))'
[[Correlations]] (unreported correlations are < 0.100)
C(bkg_amplitude, bkg_decay) = -0.999
C(peak_amplitude, peak_fraction) = 0.838
C(peak_sigma, peak_fraction) = -0.481
C(bkg_decay, peak_amplitude) = -0.338
C(bkg_amplitude, peak_amplitude) = 0.310
C(bkg_decay, peak_fraction) = -0.215
C(bkg_amplitude, peak_fraction) = 0.191
C(bkg_decay, peak_center) = -0.183
C(bkg_amplitude, peak_center) = 0.183
C(bkg_amplitude, peak_sigma) = 0.139
C(bkg_decay, peak_sigma) = -0.137
#output of integration:
(4015474.293103768, 3509959.3601876567)
C:/Users/script.py:126: IntegrationWarning: The integral is probably divergent, or slowly convergent.
THe OP's pseudo_voigt is not well formatted, but does not seem to be wrong, either, There are different definitions of the pseudo_voigt, though. Below I implemented one from Wikipedia (link in code), which usually gives good results. Looking on the logarithmic scale, however, it is not so good with this data. I also used the complex definition to get the true Voigtusing the Fedeeva function like lmfit.
Code looks as follows:
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit
from scipy.special import wofz
from scipy.integrate import quad
def cauchy(x, x0, g):
return 1. / ( np.pi * g * ( 1 + ( ( x - x0 ) / g )**2 ) )
def gauss( x, x0, s):
return 1./ np.sqrt(2 * np.pi * s**2 ) * np.exp( - (x-x0)**2 / ( 2 * s**2 ) )
# https://en.wikipedia.org/wiki/Voigt_profile#Numeric_approximations
def pseudo_voigt( x, x0, s, g, a, y0 ):
fg = 2 * s * np.sqrt( 2 * np.log(2) )
fl = 2 * g
f = ( fg**5 + 2.69269 * fg**4 * fl + 2.42843 * fg**3 * fl**2 + 4.47163 * fg**2 * fl**3 + 0.07842 * fg * fl**4+ fl**5)**(1./5.)
eta = 1.36603 * ( fl / f ) - 0.47719 * ( fl / f )**2 + 0.11116 * ( f / fl )**3
return y0 + a * ( eta * cauchy( x, x0, f) + ( 1 - eta ) * gauss( x, x0, f ) )
def voigt( x, s, g):
z = ( x + 1j * g ) / ( s * np.sqrt( 2. ) )
v = wofz( z ) #Feddeeva
out = np.real( v ) / s / np.sqrt( 2 * np.pi )
return out
def fitfun( x, x0, s, g, a, y0 ):
return y0 + a * voigt( x - x0, s, g )
if __name__ == '__main__':
xlist = np.array( [ 33.05, 33.1 , 33.15, 33.2 , 33.25, 33.3 , 33.35, 33.4 , 33.45, 33.5 , 33.55, 33.6 , 33.65, 33.7 , 33.75, 33.8 , 33.85, 33.9 , 33.95, 34. , 34.05, 34.1 , 34.15, 34.2 , 34.25, 34.3 , 34.35, 34.4 , 34.45, 34.5 , 34.55, 34.6 , 34.65, 34.7 , 34.75, 34.8 , 34.85, 34.9 , 34.95, 35. , 35.05, 35.1 , 35.15, 35.2 , 35.25, 35.3 , 35.35, 35.4 , 35.45, 35.5 , 35.55, 35.6 , 35.65, 35.7 , 35.75, 35.8 , 35.85, 35.9 , 35.95, 36. , 36.05, 36.1 , 36.15, 36.2 , 36.25, 36.3 , 36.35, 36.4 , 36.45])
ylist = np.array( [ 4569, 4736, 4610, 4563, 4639, 4574, 4619, 4473, 4488, 4512, 4474, 4640, 4691, 4621, 4671, 4749, 4657, 4751, 4921, 5003, 5071, 5041, 5121, 5165, 5352, 5304, 5408, 5393, 5544, 5625, 5859, 5851, 6155, 6647, 7150, 7809, 9017, 10967, 14122, 19529, 28029, 39535, 50684, 55730, 52525, 41356, 30015, 20345, 14368, 10736, 9012, 7765, 7064, 6336, 6011, 5806, 5461, 5283, 5224, 5221, 4895, 4980, 4895, 4852, 4889, 4821, 4872, 4802, 4928])
sol, err = curve_fit( pseudo_voigt, xlist, ylist, p0=[ 35.25,.05,.05, 30000., 3000] )
solv, errv = curve_fit( fitfun, xlist, ylist, p0=[ 35.25,.05,.05, 20000., 3000] )
print solv
xth = np.linspace( xlist[0], xlist[-1], 500)
yth = np.fromiter( ( pseudo_voigt(x ,*sol) for x in xth ), np.float )
yv = np.fromiter( ( fitfun(x ,*solv) for x in xth ), np.float )
print( quad(pseudo_voigt, xlist[0], xlist[-1], args=tuple( sol ) ) )
print( quad(fitfun, xlist[0], xlist[-1], args=tuple( solv ) ) )
solvNoBack = solv
solvNoBack[-1] =0
print( quad(fitfun, xlist[0], xlist[-1], args=tuple( solvNoBack ) ) )
fig = plt.figure()
ax = fig.add_subplot( 1, 1, 1 )
ax.plot( xlist, ylist, marker='o', linestyle='', label='data' )
ax.plot( xth, yth, label='pseudo' )
ax.plot( xth, yv, label='voigt with hack' )
ax.set_yscale('log')
plt.legend( loc=0 )
plt.show()
providing:
[3.52039054e+01 8.13244777e-02 7.80206967e-02 1.96178358e+04 4.48314849e+03]
(34264.98814344757, 0.00017531957481189617)
(34241.971907301166, 0.0002019796740206914)
(18999.266974139795, 0.0002019796990069267)
From the graph it is obvious that the pseudo_voigt is not very good. The integral, however, does not differ very much. Considering the fact that the fit optimises chi**2 this is not such a big surprise, though.
这篇关于如何使用lmfit中的优化参数查找函数(伪Voigt)的区域?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
