Numerical integration with scipy.integrate: the area under a measured peak

The problem: how much light is in one spectral line?

A spectrometer records an emission line as 401 readings of intensity, in counts, from 500 to 540 nm, one every 0.1 nm. The line sits on a sloping background, and the number you want is the area of the line alone, in counts·nm. That area is proportional to the number of emitters: to a chemist it is a concentration, to a physicist an oscillator strength. The same steps give the area of any peak on a background, whether it is a peak in a chromatogram, an X-ray diffraction peak, or a resonance in a vibration spectrum. Integrate the whole record and you get 2855 counts·nm, of which more than half is background.

scipy.integrate covers two situations. When you have samples, trapezoid and simpson integrate them as they are. When you have a formula, quad chooses its own points, returns an error estimate with the value, and copes with an infinite limit or a point where the integrand blows up. The formulas here come from the same peak and from the pendulum of solve_ivp from the ground up: the pendulum beyond small angles, whose period is the integral of a function that is infinite at one end.

Simulated emission spectrum, counts against wavelength from 500 to 540 nm. Dots: measured points. Thin line: straight baseline fitted to the flanks. Shaded: the line area above the baseline between the dashed window edges, 1253 ± 4 counts·nm.

This is where we end up: a line area of 1253 ± 4 counts·nm, against a true 1250 that the simulation knows and a real spectrometer would not. The uncertainty comes from the noise and from the baseline, and Step 5 builds the figure.

Setup

The spectrum is simulated, so the truth is known. Its constants sit at the top of the cell; after this cell nothing uses them except the reruns in Step 5, which only a simulation can do.

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import trapezoid, simpson, quad
from scipy.special import ellipk

# the truth, which a real spectrometer never shows you; only the Step 5 reruns use it
AREA = 1250.0            # line area, counts·nm
CENTER, WIDTH = 520.0, 0.8   # nm, center and standard deviation of the line
B0, B1 = 40.0, 1.5       # baseline: counts at the center, counts/nm
NOISE = 4.0              # standard deviation of the noise, counts

rng = np.random.default_rng(4)
lam = np.linspace(500, 540, 401)    # wavelength, nm
counts = (AREA / (WIDTH * np.sqrt(2 * np.pi)) * np.exp(-0.5 * ((lam - CENTER) / WIDTH) ** 2)
          + B0 + B1 * (lam - CENTER) + rng.normal(0, NOISE, lam.size))

plt.rcParams.update({
    "figure.figsize": (7, 3.6), "figure.dpi": 110,
    "axes.spines.top": False, "axes.spines.right": False,
    "axes.grid": True, "grid.alpha": 0.25,
    "font.size": 11, "lines.linewidth": 1.8,
})
INK, ACCENT, SECOND, MUTED = "#1f2a44", "#c8553d", "#2a7f9e", "#8a8f98"

print(f"{lam.size} points, spacing {lam[1] - lam[0]:.1f} nm, maximum {counts.max():.0f} counts")

401 points, spacing 0.1 nm, maximum 662 counts

Step 1: Integrate sampled data with trapezoid and simpson

Look at the record before you integrate it:

fig, ax = plt.subplots(figsize=(7, 3.2))
ax.plot(lam, counts, "o", color=INK, ms=3)
ax.text(528, 140, "background", color=MUTED)
ax.set(xlabel="λ / nm", ylabel="intensity / counts", xlim=(499.5, 540.5))
plt.show()

The line is a spike near 520 nm, a couple of nanometers wide, that reaches 662 counts with the background included. Under it the background rises from left to right, and the outer 5 nm at each end hold nothing but background.

The trapezoid rule joins neighboring points with straight lines and adds up the areas under them. Simpson's rule fits a parabola through each three neighboring points instead. Both take the samples and their positions:

print(f"trapezoid {trapezoid(counts, x=lam):8.2f} counts·nm")
print(f"simpson   {simpson(counts, x=lam):8.2f} counts·nm")

trapezoid  2855.09 counts·nm
simpson    2853.54 counts·nm

The two rules agree to 0.05 %, so the rule is not the problem. What it integrates is: the background is in there too, roughly 40 counts across 40 nm, or 1600 of the 2855 counts·nm.

Step 2: Find the peak, choose a window, subtract the baseline

The background has to come off, and the window has to come from the data, because on a real measurement you know neither where the line is nor how wide. Three moves, a line or two of code each. First, a rough baseline: a straight line through the outer 5 nm at each end. Second, the peak position, where the counts minus that rough baseline are largest, and the full width at half maximum (FWHM), the width of the peak at half its height. np.ptp is maximum minus minimum, so here it is the distance between the outermost points above half height. Third, a window of two and a half widths on each side of the peak:

edge = (lam < 505) | (lam > 535)
net = counts - np.polyval(np.polyfit(lam[edge], counts[edge], 1), lam)
peak, fwhm = lam[np.argmax(net)], np.ptp(lam[net >= net.max() / 2])
half = 2.5 * fwhm
print(f"peak at {peak:.1f} nm, FWHM {fwhm:.1f} nm, window ±{half:.1f} nm")

peak at 520.1 nm, FWHM 1.8 nm, window ±4.5 nm

The FWHM is good to one sample spacing, 0.1 nm, which is plenty for choosing a window. Whether ±2.5 FWHM keeps the whole line, Step 3 checks.

The flanks, everything outside the window, get the real baseline from np.polyfit. Measure the wavelength from the peak, so that the intercept is the background under the peak. cov=True also returns the covariance matrix C of the two coefficients, their variances and how they move together; Step 5 needs it. Subtract the baseline and integrate only the window:

win = np.abs(lam - peak) <= half + 1e-9     # the slack keeps the edge points despite rounding
flank = ~win
p, C = np.polyfit(lam[flank] - peak, counts[flank], 1, cov=True)
slope, intercept = p
d_slope, d_intercept = np.sqrt(np.diag(C))
line = counts[win] - np.polyval(p, lam[win] - peak)
area = trapezoid(line, x=lam[win])

print(f"{win.sum()} points in the window")
print(f"baseline under the peak {intercept:.2f} ± {d_intercept:.2f} counts, slope {slope:.3f} ± {d_slope:.3f} counts/nm")
print(f"line area: trapezoid {area:.2f}, simpson {simpson(line, x=lam[win]):.2f} counts·nm")

91 points in the window
baseline under the peak 40.18 ± 0.23 counts, slope 1.483 ± 0.018 counts/nm
line area: trapezoid 1252.92, simpson 1252.44 counts·nm

The background is gone, and the area is within 3 counts·nm of the true 1250. The two rules now differ by 0.48 counts·nm. Whether that is large depends on the uncertainty of the area, which Step 5 computes and Pitfall 2 puts to use.

Step 3: Integrate a formula with quad and read its error estimate

When the integrand is a formula, quad picks its own points. It cuts the interval into pieces, more of them where the function is hard, and on each piece evaluates the function at the 21 points of a Gauss-Kronrod rule, all inside the piece and none at its ends.

Step 2 left open whether a window of ±2.5 FWHM keeps the whole line. The answer is a fraction, the area inside the window over the total, so the height of the formula does not matter. Take a Gaussian of unit height with the FWHM from Step 2, which for a Gaussian is 2.355 standard deviations, and center it at zero. The total needs infinite limits:

sigma = fwhm / (2 * np.sqrt(2 * np.log(2)))

def shape(u):
    return np.exp(-0.5 * (u / sigma) ** 2)

total, err = quad(shape, -np.inf, np.inf)
print(f"total {total:.12f} nm, error estimate {err:.1e}")

total 1.916040634976 nm, error estimate 1.1e-08

On an infinite range quad maps the variable onto a finite interval and uses a 15-point rule there instead of the 21. The mapped samples crowd around zero, which is why the shape sits there; Pitfall 1 shows what happens when it does not. Now the three windows:

for k in [1, 1.5, 2.5]:
    inside, err = quad(shape, -k * fwhm, k * fwhm)
    print(f"±{k:3.1f} FWHM: lost fraction {1 - inside / total:.1e}, error estimate {err:.0e}")

±1.0 FWHM: lost fraction 1.9e-02, error estimate 3e-09
±1.5 FWHM: lost fraction 4.1e-04, error estimate 5e-14
±2.5 FWHM: lost fraction 3.9e-09, error estimate 2e-14

A window of ±1 FWHM loses 1.9 % of the area, ±1.5 FWHM four parts in ten thousand, ±2.5 FWHM four parts in a billion: the window from Step 2 keeps the whole line for any Gaussian. A line of another shape goes into shape and gets the same question. The error estimates, 3e-09 and below, measure how well quad integrated the formula, not the noise in your data.

The tolerances epsabs and epsrel both default to 1.49e-8; quad stops when its error estimate falls below the larger of epsabs and epsrel times the value. full_output=1 adds a dictionary whose neval is the cost:

for tol in [1.49e-8, 1e-3]:
    value, err, info = quad(shape, -half, half, epsabs=tol, epsrel=tol, full_output=1)
    print(f"tolerance {tol:.2g}: {info['neval']:3d} evaluations, value {value:.12f} nm")

tolerance 1.5e-08: 147 evaluations, value 1.916040627443 nm
tolerance 0.001:  63 evaluations, value 1.916040627443 nm

Loosening to 1e-3 cuts the work from 147 to 63 evaluations and leaves all twelve printed digits alone. Keep the defaults unless one evaluation of your function is expensive.

Step 4: Handle singular points: the pendulum period

The period of a pendulum released from rest at \(\theta_0\) follows from energy conservation per unit mass, the mass having canceled:

\[\tfrac12 L^2\dot\theta^2 = gL(\cos\theta - \cos\theta_0) .\]

A quarter period is the time from \(\theta = 0\) to \(\theta_0\), the integral of \(d\theta/\dot\theta\), so

\[T = 4\sqrt{\frac{L}{2g}}\int_0^{\theta_0}\frac{d\theta}{\sqrt{\cos\theta - \cos\theta_0}} .\]

The integrand is infinite at \(\theta_0\), where the bob stops. The integral is still finite: near \(\theta_0\) the integrand behaves like \(1/\sqrt{\theta_0 - \theta}\), whose area is finite. quad takes it because it never evaluates the endpoints, and its subdivision crowds pieces toward the trouble. The check is the exact period, \(4\sqrt{L/g}\,K(m)\) with \(m = \sin^2(\theta_0/2)\), where \(K\) is the complete elliptic integral of the first kind and scipy.special.ellipk takes \(m\) as its argument:

g, L = 9.81, 1.0                  # m/s², m
theta0 = np.radians(90)
I, err, info = quad(lambda theta: 1 / np.sqrt(np.cos(theta) - np.cos(theta0)), 0, theta0, full_output=1)
T = 4 * np.sqrt(L / (2 * g)) * I
T0 = 2 * np.pi * np.sqrt(L / g)
T_exact = 4 * np.sqrt(L / g) * ellipk(np.sin(theta0 / 2) ** 2)

print(f"T = {T:.6f} s, small-angle T0 = {T0:.3f} s, T/T0 = {T / T0:.4f}, {info['neval']} evaluations")
print(f"error estimate {4 * np.sqrt(L / (2 * g)) * err:.1e} s, actual error {T - T_exact:.1e} s")

T = 2.367842 s, small-angle T0 = 2.006 s, T/T0 = 1.1803, 231 evaluations
error estimate 2.5e-09 s, actual error -8.6e-14 s

That is the 2.367842 s that solve_ivp found with an event in solve_ivp from the ground up, 18 % above the small-angle formula, from one call. The error estimate is an upper bound, and here it overstates the actual error thirty thousand times.

A singular point inside the interval is a different matter. The function \(1/\sqrt{|x|}\) on \((-1, 1)\) has an area of exactly 4, and the middle of the interval is one of the 21 points of the rule:

f = lambda x: 1 / np.sqrt(np.abs(x))
for pts in [None, [0]]:
    with np.errstate(divide="ignore"):      # 1/0 becomes inf without a NumPy warning
        value, err = quad(f, -1, 1, points=pts)
    print(f"points={pts}: value {value}, error estimate {err:.0e}")

points=None: value inf, error estimate inf
points=[0]: value 3.9999999999999813, error estimate 6e-14

Without points, the rule evaluates the function at zero, NumPy turns 1/0 into infinity, and quad hands back infinity as both value and error estimate. With points=[0] quad splits the interval at zero, the singularity sits at an endpoint again, and the answer is 4 to 14 digits. Two calls on \((-1, 0)\) and \((0, 1)\) do the same. quad does not always land on an interior singularity, but you cannot rely on it missing: when you know where the integrand blows up, say so.

Step 5: Put an uncertainty on the peak area

The trapezoid area is a weighted sum of the counts, \(A = \sum_i w_i y_i\), with weights fixed by the wavelengths. Row \(i\) of an identity matrix is a spectrum that is 1 at point \(i\) and 0 elsewhere, so its area is \(w_i\):

w = trapezoid(np.eye(win.sum()), x=lam[win])
print("weights / nm:", w[:3], "...", w[-2:])

weights / nm: [0.05 0.1  0.1 ] ... [0.1  0.05]

The weights are 0.1 nm inside the window and 0.05 nm at its ends. Two things make \(A\) uncertain. First, the noise. Its standard deviation \(\sigma_y\) comes from the scatter of the flank points around the baseline, with ddof=2 because two parameters were fitted. The variance of a weighted sum of independent values is \(\sigma_y^2\sum_i w_i^2\), so the noise contributes \(\sigma_y\sqrt{\sum_i w_i^2}\).

Second, the baseline. Its area under the window is slope × \(\sum_i w_i(\lambda_i - \lambda_\text{peak})\) + intercept × \(\sum_i w_i\), a weighted sum of the two coefficients with weights \(\mathbf g = \big(\sum_i w_i(\lambda_i - \lambda_\text{peak}),\ \sum_i w_i\big)\). Its variance is \(\mathbf g^\mathsf{T} C\,\mathbf g\), where the off-diagonal entries of \(C\) carry the correlation of slope and intercept:

resid = counts[flank] - np.polyval(p, lam[flank] - peak)
sigma_y = np.std(resid, ddof=2)
dA_noise = sigma_y * np.sqrt(np.sum(w ** 2))

gvec = np.array([np.sum(w * (lam[win] - peak)), np.sum(w)])
dA_base = np.sqrt(gvec @ C @ gvec)
dA = np.hypot(dA_noise, dA_base)

print(f"noise sigma_y {sigma_y:.2f} counts -> noise term {dA_noise:.2f} counts·nm")
print(f"g = ({gvec[0]:.1f}, {gvec[1]:.1f}) nm -> baseline term {dA_base:.2f} counts·nm")
print(f"A = {area:.1f} ± {dA:.1f} counts·nm")
w_simpson = simpson(np.eye(win.sum()), x=lam[win])
print(f"noise term with Simpson weights {sigma_y * np.sqrt(np.sum(w_simpson ** 2)):.2f} counts·nm")

noise sigma_y 4.07 counts -> noise term 3.85 counts·nm
g = (0.0, 9.0) nm -> baseline term 2.08 counts·nm
A = 1252.9 ± 4.4 counts·nm
noise term with Simpson weights 4.07 counts·nm

The vector \(\mathbf g\) is (0, 9.0 nm) because the window is symmetric about the point the baseline was measured from, so the baseline term is the window width times the intercept's uncertainty, 9.0 nm × 0.23 counts. Flank and window share no points, so the two terms are independent and add in squares: 1252.9 ± 4.4 counts·nm, 0.7 standard deviations from the true 1250.

A simulation can check the 4.4 by rerunning the measurement 2000 times, each with its own noise, baseline fit, and area over the same window. A real measurement cannot:

R = 2000
reruns = (AREA / (WIDTH * np.sqrt(2 * np.pi)) * np.exp(-0.5 * ((lam - CENTER) / WIDTH) ** 2)
          + B0 + B1 * (lam - CENTER) + rng.normal(0, NOISE, (R, lam.size)))
P = np.polyfit(lam[flank] - peak, reruns[:, flank].T, 1)         # one fit per rerun
baselines = P[0][:, None] * (lam[win] - peak) + P[1][:, None]
areas = trapezoid(reruns[:, win] - baselines, x=lam[win], axis=1)
print(f"{R} reruns: mean {areas.mean():.2f}, standard deviation {areas.std():.2f} counts·nm")

2000 reruns: mean 1249.97, standard deviation 4.25 counts·nm

The reruns scatter by 4.25 counts·nm around 1249.97, so the propagated 4.4 is right to 3 % and the method has no bias. The figure from the top puts the result on the spectrum:

baseline = np.polyval(p, lam - peak)
fig, ax = plt.subplots()
ax.plot(lam, counts, "o", color=INK, ms=3, mew=0, alpha=0.5)
ax.plot(lam, baseline, color=SECOND, lw=1.2, zorder=3)    # drawn over the dots
ax.fill_between(lam[win], baseline[win], counts[win], color=ACCENT, alpha=0.35, lw=0)
for edge_lam in [peak - half, peak + half]:
    ax.axvline(edge_lam, color=MUTED, ls="--", lw=1)
ax.text(502, 80, "baseline", color=SECOND)
ax.text(peak + half + 0.8, 400, f"area {area:.0f} ± {dA:.0f} counts·nm", color=ACCENT)
ax.set(xlabel="λ / nm", ylabel="intensity / counts", xlim=(499.5, 540.5))
plt.show()

Pitfalls

quad returns zero for a peak it never saw. Move the Step 3 shape to the peak position and integrate it over infinite limits, or over 0 to 2000 nm:

for label, a, b in [("-inf, inf", -np.inf, np.inf), ("0, 2000", 0, 2000), ("peak ± 8", peak - 8, peak + 8)]:
    value, err = quad(lambda x: shape(x - peak), a, b)
    print(f"{label:10s} {value:.5f} nm, error estimate {err:.0e}")

-inf, inf  0.00000 nm, error estimate 0e+00
0, 2000    0.00000 nm, error estimate 0e+00
peak ± 8   1.91604 nm, error estimate 7e-10

Zero, with an error estimate of zero. quad knows the function only through its samples, and on a range a thousand times wider than the line none of them lands in the 2 nm where the line lives. The error estimate comes from the same samples, so it agrees. Give limits that hug the peak, as in the third line, or shift the variable so that the peak sits at zero before you use infinite limits. An estimate of zero for an integral you expect to be nonzero means quad never saw your function.

Reporting an integration error as the uncertainty. The area gets written as 1252.9 ± 0.5 counts·nm, with 0.5 the difference between trapezoid and Simpson, or with the even smaller error estimate of quad run on a spline drawn through the points (CubicSpline from scipy.interpolate). Either understates the uncertainty by a factor of nine or more. Both numbers measure how well a rule integrates the points you have, and neither knows that the points themselves would come out differently on a rerun, which is what Step 5's 4.4 counts·nm measures. A spline passes through every noisy point, so integrating it precisely integrates the noise precisely. Report the propagated uncertainty and compare the rule difference with it; the choice of rule matters only when the two are comparable. A higher-order rule does not shrink the uncertainty either: Simpson's alternating weights raise the noise term from 3.85 to 4.07 counts·nm.

Leaving out x. Integrate the counts on their own and the area is ten times too large:

print(f"{trapezoid(counts):.1f}")

28550.9

Without x or dx the spacing is taken as 1, so the result is in counts·sample. Always pass x. The unit of an integral is the unit of \(y\) times the unit of \(x\), counts·nm here, and the same line on a wavenumber or energy axis has a different area.

Variations

A running integral. cumulative_trapezoid(y, x=x, initial=0) gives the area up to every point: the energy delivered, from a power log, or the wavelength where half of a line's area is reached.
A Lorentzian or Voigt line. A Voigt line is a Gaussian convolved with a Lorentzian. The wings reach far beyond any window. Run the Step 3 loop on a Lorentzian shape and ±2.5 FWHM keeps only \((2/\pi)\arctan 5\), 87 %, of the area. Fit the line shape instead, as in Fit a curve to data with error bars and draw a confidence band, and integrate the fitted model with quad on infinite limits, centered as in Step 3.
Two dimensions. dblquad(f, a, b, c, d) for the power in a beam spot whose intensity is a formula, with f(y, x) taking the inner variable first. For an image, apply trapezoid twice, once along each axis.
Oscillating integrands. quad(f, 0, np.inf, weight="cos", wvar=omega) for Fourier-type integrals, where a plain quad struggles with endless sign changes.

Cheat sheet

trapezoid(y, x=x)                       # always pass x; unit is unit(y) * unit(x)
simpson(y, x=x)                         # smooth data, few points
fwhm = np.ptp(x[net >= net.max() / 2])  # width from the data; window ±2.5 FWHM for a Gaussian
p, C = np.polyfit(x[flank] - peak, y[flank], 1, cov=True)   # baseline; sigma_y from its residuals, ddof=2
value, err = quad(f, a, b)              # err is numerical, not your uncertainty
quad(f, -np.inf, np.inf)                # center the peak at zero first
quad(f, a, b, points=[c])               # singular point c inside (a, b)
w = trapezoid(np.eye(n), x=x[win])      # n = win.sum(); the area is sum(w * y[win])
g = np.array([np.sum(w * (x[win] - peak)), np.sum(w)])   # baseline area is g @ p
dA = np.hypot(sigma_y * np.sqrt(np.sum(w**2)), np.sqrt(g @ C @ g))   # noise, baseline