scipy.integrate伪Voigt函数，积分变为0

Question

我正在使用 Scipy、Numpy 和 Matplotlib 编写一个脚本，用于将峰形拟合到 Python 中的光谱数据。它可以同时拟合多个峰。峰值轮廓（目前）是 Pseudo-Voigt，它是高斯（又名正态）和洛伦兹（又名柯西）分布的线性组合。

我有一个选项开关，我可以通过它让软件优化高斯和洛伦兹的贡献或将其设置为固定值（其中 0 = 纯高斯和 1 = 纯洛伦兹）。正常工作，绘制拟合峰看起来符合预期。当我尝试使用 scipy.integrate 计算峰的积分时，问题就开始了。

到目前为止，我尝试了 scipy.integrate.quad、scipy.integrate.quadrature、scipy.integrate.fixed_quad 和 scipy.integrate.romberg。当峰是纯高斯时，积分变成类似于 1.73476E-34（不总是相同的数字），即使峰的面积明显大于非纯高斯但 return 的相邻峰的有限积分大约在 10 到 1000 之间。相关部分如下所示：

# Function defining the peak functions for plotting and integration
# WavNr: Wave number, the x-axis over which shall be integrated
# Pos: Peak center position
# Amp: Amplitude of the peak
# GammaL: Gamma parameter of the Lorentzian distribution
# FracL: Fraction of Lorentzian distribution
def PseudoVoigtFunction(WavNr, Pos, Amp, GammaL, FracL):
    SigmaG = GammaL / np.sqrt(2*np.log(2)) # Calculate the sigma parameter  for the Gaussian distribution from GammaL (coupled in Pseudo-Voigt)
    LorentzPart = Amp * (GammaL**2 / ((WavNr - Pos)**2 + GammaL**2)) # Lorentzian distribution
    GaussPart = Amp * np.exp( -((WavNr - Pos)/SigmaG)**2) # Gaussian distribution
    Fit = FracL * LorentzPart + (1 - FracL) * GaussPart # Linear combination of the two parts (or distributions)
    return Fit

这是绘图函数通过以下方式调用的：

Fit = PseudoVoigtFunction(WavNr, Pos, Amp, GammaL, FracL)

效果很好。集成商也通过以下方式调用它：

PeakArea, PeakAreaError = integrate.quad(PseudoVoigtFunction, -np.inf, np.inf, args=(Pos, Amp, GammaL, FracL))

或 scipy.integrate 提供的任何其他变体，都具有相同的结果，如果 FracL = 0，则 PeakArea =（几乎）0。

我敢肯定问题出在我太笨了，无法弄清楚 scipy.integrate 如何处理比我能找到的示例稍微复杂的函数。希望有人看到我没有看到的明显错误。搜索 Whosebug 和 Scipy 文档以及重新排列和完全重写我的代码的两天时间让我无所适从。我怀疑 scipy.integrate 中的参数与问题有某种关系，但据我所知，它们似乎排列正确。

提前致谢， Os

Answer 1

我相信您已经知道，间隔 (-inf, inf) 相当大。 :) 高斯衰减得非常快，所以除了峰值附近的间隔外，高斯在数值上与 0 没有区别。我怀疑 quad 根本没有看到你的峰值。

一个简单的解决方法是将积分分成两个区间，(-inf, pos) 和 (pos, inf)。（你的函数关于 Pos 是对称的，所以你实际上只需要两倍于 (-inf, pos) 的积分。）

举个例子。我不知道这些参数值是否接近您使用的典型值，但它们说明了这一点。

In [259]: pos = 1500.0

In [260]: amp = 4.0

In [261]: gammal = 0.5

In [262]: fracl = 0  # Pure Gaussian

quad认为积分为0:

In [263]: quad(PseudoVoigtFunction, -np.inf, np.inf, args=(pos, amp, gammal, fracl))
Out[263]: (0.0, 0.0)

相反，整合 (-inf, pos) 和 (pos, inf)：

In [264]: quad(PseudoVoigtFunction, -np.inf, pos, args=(pos, amp, gammal, fracl))
Out[264]: (1.5053836955785238, 3.616268258191726e-11)

In [265]: quad(PseudoVoigtFunction, pos, np.inf, args=(pos, amp, gammal, fracl))
Out[265]: (1.5053836955785238, 3.616268258191726e-11)

因此 (-inf, inf) 的积分约为 3.010767。

scipy.integrate伪Voigt函数，积分变为0

scipy.integrate Pseudo-Voigt function, integral becomes 0

python

numpy

integrate

scipy

data-fitting