将 Voigt 函数拟合到 Python 中的数据答案

【问题标题】：Fitting Voigt function to data in Python将 Voigt 函数拟合到 Python 中的数据
【发布时间】：2019-11-27 06:22:40
【问题描述】：

我最近运行了一个脚本，用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()

【问题讨论】：

请您添加 x 和 y 的值，好吗？还有，xx。
对不起！现在完成
实际上，线形应该是洛伦兹的，它们会被扩展为高斯。所以，现在你正在做的不是扩大你的高斯，而是让它更紧。
代码中缺少导入，请提供完整的工作示例！
添加了有关问题和导入的更多信息。

标签： python curve-fitting

【解决方案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。

适合的情节如下所示：

【讨论】：

这看起来不错！最后3-4点为背景，其余为线。（它非常次优，但不幸的是，它是我们唯一拥有的数据集。是否可以告诉它使用前 3-4 个点作为连续统一体并告诉它不要超过它？我也玩过权重，因为中心点在光学上很厚，不能用 voigt 建模。我个人最感兴趣的是轮廓的偏移和宽度，所以我认为在这种情况下是伽玛。
我上面说了，我可能听不懂你。最后 3-4 个点在后台？您的意思是峰值以 1 或 2 为中心并且具有正幅度？或者 first 3-4 个峰值是背景，峰值以 8 为中心并且具有负幅值？当然，这表明您在理解您正在建模的信号方面还有一些工作要做。与理解模型相比，拟合的机制相对简单。另外：对于 Voigt 轮廓，gamma 给出了类似洛伦兹尾巴的宽度。 “宽度的偏移”不清楚。

【解决方案2】：

我设法得到了一些东西，但不是很满意。如果您将偏移量作为参数删除并直接在Voigt 函数中添加20000，起始值[8, 126000, 0.71, 2]（特定值无关紧要）你会得到类似

现在，拟合产生了 gamma 的值，这是我无法证明的负数。我希望gamma 总是积极的，但也许我错了，这完全没问题。

您可以尝试的一件事是镜像您的数据，使其成为“正”峰值（同时消除背景）和/或标准化值。这可能会帮助您实现融合。

我不知道为什么当使用偏移量作为参数时，求解器在寻找最优值时会遇到问题。也许您需要不同的优化器例程。

也许使用 lmfit 包会是一个更好的选择，它是 scipy 的一个包装器，可以用许多预建的线形拟合非线性函数。甚至还有一个example 适合Voigt 配置文件。

【讨论】：