我正在使用 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
(并不总是相同的数字),即使峰值的面积明显大于非纯高斯分布的相邻峰值,但会返回 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 如何与比我能找到的示例稍微复杂的函数一起工作。希望有人看到我没有看到的明显错误。两天的搜索 stackoverflow 和 Scipy Docs 以及重新排列和完全重写我的代码却毫无进展。我怀疑 scipy.integrate 中的参数在某种程度上与问题有关,但就我所知,它们似乎排列正确。
提前致谢, 操作系统
我相信您知道,间隔(-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。