I am writing a script for fitting peak shapes to spectroscopic data in Python, using Scipy, Numpy and Matplotlib. It can fit multiple peaks at once. The peak profile (for now) is Pseudo-Voigt, which is a linear combination of a Gaussian (aka Normal) and Lorentzian (aka Cauchy) distribution.
I have an option switch with which I can either let the software optimize the contribution of Gaussian and Lorentzian or set it to a fixed value (where 0 = pure Gaussian and 1 = pure Lorentzian). Works as it should, plotting the fitted peaks looks as expected. The problem starts, when I try to calculate the integrals of the peaks using scipy.integrate.
So far I tried scipy.integrate.quad, scipy.integrate.quadrature, scipy.integrate.fixed_quad and scipy.integrate.romberg. The integral becomes something like 1.73476E-34 (not always the same number) when the peak is pure Gaussian, even for peaks that clearly have larger areas than neighbouring peaks that are not pure Gaussian but return a finite integral of something on the order of 10 to 1000. Here's what the relevant parts look like:
# 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
This is the called by the plot function via:
Fit = PseudoVoigtFunction(WavNr, Pos, Amp, GammaL, FracL)
which works fine. It is also called by the integrator via:
PeakArea, PeakAreaError = integrate.quad(PseudoVoigtFunction, -np.inf, np.inf, args=(Pos, Amp, GammaL, FracL))
or any of the other variants supplied by scipy.integrate, all with the same result that if FracL = 0, then PeakArea = (almost) 0.
I'm sure the problem is just me being too stupid to figure out how scipy.integrate works with slightly more complicated functions than I can find examples for. Hoping someone sees the obvious error that I don't. Two days of searching stackoverflow and Scipy Docs and rearranging and completely rewriting my code have taken me nowhere. I'm suspecting that the args in the scipy.integrate are somehow related to the problem but for all I could find out, they seem to be arranged correctly.
Thanks in advance, Os
As I'm sure you are aware, the interval (-inf, inf) is pretty big. :) A Gaussian decays very fast, so except for an interval near the peak, a Gaussian is numerically indistinguishable from 0. I suspect quad is simply not seeing your peak.
A simple work-around is to split the integral into two intervals, (-inf, pos) and (pos, inf). (Your function is symmetric about Pos, so you really only need twice the integral over (-inf, pos).)
Here's an example. I don't know if these parameter values are anywhere near the typical values you use, but they illustrate the point.
In [259]: pos = 1500.0
In [260]: amp = 4.0
In [261]: gammal = 0.5
In [262]: fracl = 0 # Pure Gaussian
quad thinks the integral is 0:
In [263]: quad(PseudoVoigtFunction, -np.inf, np.inf, args=(pos, amp, gammal, fracl))
Out[263]: (0.0, 0.0)
Instead, integrate over (-inf, pos) and (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)
So the integral over (-inf, inf) is approximately 3.010767.