如何使用lmfit中的优化参数来查找函数(伪Voigt)的区域?

问题描述 投票:0回答:1

我正在尝试确定曲线的面积(峰值)。我能够使用伪Voigt配置文件和指数背景成功拟合峰(数据),并获得与使用商业软件获得的参数相符的拟合参数。现在的问题是试图将那些拟合的峰参数与峰面积相关联。

与高斯线形的情况不同,我找不到使用拟合参数计算峰面积的简单方法。所以我正在尝试使用scipy quad函数来集成我的拟合函数。我知道该区域应该由商业软件确定为19,000左右,但是我得到的数值非常大。

拟合效果很好(通过绘图确认...),但是计算区域不紧密。尝试绘制传递有最佳拟合值的psuedo_voigt_func函数后,我发现它的峰过于密集。这样,积分可能是正确的,那么错误将出在我如何创建峰的过程中,即如何将拟合的参数传递给我的psuedo_voigt_func函数,该函数是从lmfit模型网站(https://lmfit.github.io/lmfit-py/builtin_models.html)转录而来的。我相信我已正确编写psuedo voigt函数的脚本,但无法正常工作。

#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.
python curve-fitting integrate lmfit
1个回答
0
投票

lmfit中的其他峰状线型和模型一样,amplitude参数应给出该分量的面积。

© www.soinside.com 2019 - 2024. All rights reserved.