在Python中对Voigt函数进行数据拟合

问题描述

我最近运行了一个脚本，使用help of SO将高斯函数与我的吸收配置文件相匹配。我希望如果我简单地用Voigt函数替换Gauss函数，事情就会运行得很好，但情况似乎并非如此。我认为这主要是因为它是一款变速的Voigt。

编辑：轮廓是光学厚度不同的吸收线。在实践中，它们将是光学厚薄特征的混合体。就像这张图的底部。当前的数据将更像顶部的图像，但可能底部已经变平了一点。(我们只能看到轮廓的左侧，略高于中心)

对于高斯来说，它看起来是这样的，正如预测的那样，底部似乎没有拟合所希望的那么深，但仍然相当接近。不过，个人资料本身应该仍然是Voigt。但现在我意识到，这些中心点可能会让人感觉不对劲。那么也许应该根据机翼的位置来增加重量呢？

我最想知道的是移位函数是否定义错误，或者它是否是我的起始值。

import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import numpy as np
from scipy.special import wofz

x = np.arange(13)
xx = xx = np.linspace(0, 13, 100)
y = np.array([19699.959 , 21679.445 , 21143.195 , 20602.875 , 16246.769 ,
        11635.25  ,  8602.465 ,  7035.493 ,  6697.0337,  6510.092 ,
              7717.772 , 12270.446 , 16807.81  ])
# weighted arithmetic mean (corrected - check the section below)
#mean = 2.4
sigma = 2.4
gamma = 2.4

def Gauss(x, y0, a, x0, sigma):
    return y0 + a * np.exp(-(x - x0)**2 / (2 * sigma**2))


def Voigt(x, x0, y0, a, sigma, gamma):
    #sigma = alpha / np.sqrt(2 * np.log(2))

    return y0 + a * np.real(wofz((x - x0 + 1j*gamma)/sigma/np.sqrt(2))) / sigma /np.sqrt(2*np.pi)


popt, pcov = curve_fit(Voigt, x, y, p0=[8, np.max(y), -(np.max(y)-np.min(y)), sigma, gamma])
#p0=[8, np.max(y), -(np.max(y)-np.min(y)), mean, sigma])

plt.plot(x, y, 'b+:', label='data')
plt.plot(xx, Voigt(xx, *popt), 'r-', label='fit')

plt.legend()
plt.show()

---

解决方案

我可能误解了您正在使用的模型，但我认为您需要包含某种恒定或线性背景。

要使用lmfit(它内置了Voigt、Gauss和许多其他模型，并非常努力地使这些模型可互换)，我建议从以下内容开始：

import numpy as np
import matplotlib.pyplot as plt
from lmfit.models import GaussianModel, VoigtModel, LinearModel, ConstantModel

x = np.arange(13)
xx = np.linspace(0, 13, 100)
y = np.array([19699.959 , 21679.445 , 21143.195 , 20602.875 , 16246.769 ,
              11635.25  ,  8602.465 ,  7035.493 ,  6697.0337,  6510.092 ,
              7717.772 , 12270.446 , 16807.81  ])

# build model as Voigt + Constant
## model = GaussianModel() + ConstantModel()
model = VoigtModel() + ConstantModel()

# create parameters with initial values
params = model.make_params(amplitude=-1e5, center=8, 
                           sigma=2, gamma=2, c=25000)

# maybe place bounds on some parameters
params['center'].min = 2
params['center'].max = 12
params['amplitude'].max = 0. 

# do the fit, print out report with results 
result = model.fit(y, params, x=x)
print(result.fit_report())

# plot data, best fit, fit interpolated to `xx`
plt.plot(x, y, 'b+:', label='data')
plt.plot(x, result.best_fit, 'ko', label='fitted points')
plt.plot(xx, result.eval(x=xx), 'r-', label='interpolated fit')
plt.legend()
plt.show()

是的，您只需将VoigtModel()替换为GaussianModel()或LorentzianModel()，然后重新进行拟合并比较拟合统计数据，以确定哪个模型更好。

对于Voigt型号，打印的报告将为

[[Model]]
    (Model(voigt) + Model(constant))
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 41
    # data points      = 13
    # variables        = 4
    chi-square         = 17548672.8
    reduced chi-square = 1949852.54
    Akaike info crit   = 191.502014
    Bayesian info crit = 193.761811
[[Variables]]
    amplitude: -173004.338 +/- 30031.4068 (17.36%) (init = -100000)
    center:     8.06574198 +/- 0.16209266 (2.01%) (init = 8)
    sigma:      1.96247322 +/- 0.23522096 (11.99%) (init = 2)
    c:          23800.6655 +/- 1474.58991 (6.20%) (init = 25000)
    gamma:      1.96247322 +/- 0.23522096 (11.99%) == 'sigma'
    fwhm:       7.06743644 +/- 0.51511574 (7.29%) == '1.0692*gamma+sqrt(0.8664*gamma**2+5.545083*sigma**2)'
    height:    -18399.0337 +/- 2273.61672 (12.36%) == '(amplitude/(max(2.220446049250313e-16, sigma*sqrt(2*pi))))*wofz((1j*gamma)/(max(2.220446049250313e-16, sigma*sqrt(2)))).real'
[[Correlations]] (unreported correlations are < 0.100)
    C(amplitude, c)     = -0.957
    C(amplitude, sigma) = -0.916
    C(sigma, c)         =  0.831
    C(center, c)        = -0.151

请注意，默认情况下，gamma被约束为与sigma相同的值。可以取消该约束，并使gamma与params['gamma'].set(expr=None, vary=True, min=1.e-9)独立变化。我认为您在此数据集中可能没有足够的数据点来稳健而独立地确定gamma。

该拟合的曲线图如下所示：