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

Question

我正在尝试确定曲线的面积（峰值）。我能够使用 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])

Answer 1

OP 的pseudo_voigt格式不正确，但似乎也没有错，但是 有不同的定义pseudo_voigt。下面我从 Wikipedia 实现了一个（代码中的链接），通常会产生很好的结果。然而，从对数尺度来看，这个数据并不是那么好。我还使用复杂的定义来获得真正的 VoigtusingFedeeva函数，例如lmfit.代码如下所示：import matplotlib.pyplot as pltimport numpy as npfrom scipy.optimize import curve_fitfrom scipy.special import wofzfrom scipy.integrate import quaddef 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_approximationsdef 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 outdef 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这一事实并不是一个大惊喜。

Answer 2

与 中的其他类似峰的线形和模型一样lmfit，该amplitude参数应给出该组件的面积。对于它的价值，使用真正的 Voigt 函数而不是伪 Voigt 函数可能会更好。