python - 将 Voigt 函数拟合到 Python 中的数据

Question

我最近运行了一个脚本，用 help of SO 将高斯拟合到我的吸收曲线中。 .我希望如果我简单地用 Voigt 函数替换高斯函数，一切都会很好，但事实似乎并非如此。我认为主要是因为它是一个转移的 voigt。

编辑:分布图是光学厚度不同的吸收线。实际上，它们将是光学厚特征和薄特征的混合体。就像这张图中的底部一样。当前数据将更像顶部图像，但底部图像可能已经变平了一点。 (而且我们只看到侧面的左侧，有点超出中心)

对于 Gauss，它看起来像这样，正如预测的那样，底部似乎没有合身所希望的那么深，但仍然非常接近。配置文件本身应该仍然是一个 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()

Answer 1

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

要使用 lmfit(它内置了 Voigt、Gaussian 和许多其他模型，并且非常努力地使它们可以互换)，我建议从这样的事情开始:

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。

该拟合的图如下所示: