如何使用 lmfit 的优化参数找到函数的面积（伪 Voigt）？答案

【问题标题】：How to find the area of a function (Pseudo Voigt) using optimized parameters from lmfit?如何使用 lmfit 的优化参数找到函数的面积（伪 Voigt）？
【发布时间】：2019-10-17 18:53:07
【问题描述】：

我正在尝试确定曲线的面积（峰值）。我能够使用 Pseudo Voigt 轮廓和指数背景成功拟合峰值（数据），并获得与使用商业软件获得的参数一致的拟合参数。现在的问题是试图将这些拟合峰参数与峰面积联系起来。

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

拟合效果很好（通过绘图确认...），但计算区域并不接近。在尝试使用传递给它的最佳拟合值绘制我的 psuedo_voigt_func 函数后，我发现它是一个过于强烈的峰值。这样，集成可能是正确的，那么错误将在于我如何创建峰值，即通过将拟合参数传递给我的 psuedo_voigt_func 函数，其中该函数是从 lmfit 模型网站 (https://lmfit.github.io/lmfit-py/builtin_models.html) 转录的。我相信我正确编写了伪 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.

【问题讨论】：

嗨，我的印象是恒定的背景效果很好。使用指数模型有什么特别的原因吗？ ...( -3.8e7 +/- 5.6e12 ???)
实际上，您的 Pseud-Voigt 中有两个错误：term2 应为：np.exp( -(x - cen)**2 / (2 * sig_gauss**2) ) ter3 应为 ((amp*alpha) / np.pi) * ( sig / ((x - cen)**2 + sig**2) )。因此，出现了两次括号问题。
感谢您的 cmets。我决定不使用 lmfit 模块，而只是写出我自己的方程式，我现在有一个工作脚本。因为我只是在尝试不同的数据，所以没有特别的理由说明对峰值数据的指数基线贡献。我确定了基线项的二次方和峰值的 Voigt 函数。

标签： python curve-fitting integrate lmfit

【解决方案1】：

OP 的pseudo_voigt 格式不正确，但似乎也没有错，不过pseudo_voigt 有不同的定义。下面我从 Wikipedia 实现了一个（代码中的链接），通常会产生很好的结果。然而，从对数尺度来看，这个数据并不是那么好。我还使用复杂的定义来获得真正的 Voigtusing Fedeeva 函数，如 lmfit。

代码如下：

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()

提供：

[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 的事实，这并不是什么大惊喜。

【讨论】：

注意拟合结果显示，sigma 和gamma 非常相似。我的印象是，如果一个占主导地位的 Pseudo-Voigt 函数可能会更好地工作，这里解释了“失败”。
非常感谢您的所有工作！这看起来很棒，而且比我能做的要好得多。我绝对可以将它用于我的代码。我看到的唯一问题是它返回的区域比预期的要大，我想我已经弄清楚了原因。这是在基线（y = 0）处与背景进行整合，但是我需要相对于数据基线的峰面积，这就是我包含 exp 背景模型的原因。我用商业软件证实了这一点，前者和后者的情况分别为 34239 和 18704（正确）。背景如何融入？
@KillianTallman 检查更新。我使用恒定的背景。减法很容易，积分约为 19000。您可以轻松修改代码以使用指数背景。干杯。
哇，太棒了！这就是我一直在寻找的。为了让我和其他观众更容易理解，您能否添加一些 cmets 或使变量名称更具体（例如 a 可以是 alpha 或amplitude？）。主要只需要辅助函数中的变量（我已经掌握得很好）。也许还有一点关于你用 solvNoBack 做了什么。看起来您正在将 y0（不完全确定该变量代表什么）设置为零，因此它不会传递给 Voigt 函数进行积分。任何添加的 cmets 将不胜感激！！！！！！
意识到原始问题在打印结果中有正确答案可能是有益的：该区域确实是 ~19000，与 lmfit 报告的完全一样，由于细节的原因存在一些不确定性使用的背景。这里的答案并没有错，但它比使用 lmfit 的原始帖子要工作得多（而且不太清楚）。

【解决方案2】：

与lmfit 中的其他类似峰的线形和模型一样，amplitude 参数应该给出该组件的面积。

就其价值而言，使用真正的 Voigt 函数可能比使用伪 Voigt 函数更好。

【讨论】：