Numerical Integration

Why Does This Exist?

Most integrals in ML have no closed form: Bayesian posteriors, expected losses over continuous noise, partition functions. Analytic integration, taught in integration and area, handles the textbook cases; everything else needs numerical integration (quadrature): approximate $\int_a^b f$ by a weighted sum $\sum_i w_i f(x_i)$ of function values at chosen points.

The design space is one trade: where to put the points and how to weight them. Trapezoids join values with straight lines (cheap, robust). Simpson fits parabolas (exact for cubics and below). Gaussian quadrature optimises both points and weights (maximum exactness per evaluation). Monte Carlo scatters points randomly (dimension-proof, slow $1/\sqrt{N}$ convergence). Expectations over distributions connect in Monte Carlo methods.

The one idea you need first

A Riemann sum slices $[a, b]$ into panels and totals rectangle areas: the integral is its limit as panels thin. Quadrature rules are smarter panels: trapezoid averages the two endpoint heights per panel, Simpson weights endpoint-midpoint-endpoint as $1$-$4$-$1$ per double panel. Same limit, far fewer panels.

Think of It Like This

How It Actually Works

On $[0, 1]$ with $f(x) = x^2$ and two panels ($h = 0.5$, points $0, 0.5, 1$):

• Trapezoid: $h \cdot (f_0/2 + f_1 + f_2/2) = 0.5 \cdot (0 + 0.25 + 0.5) = 0.375$. Error $O(h^2)$: overshoots because straight tops sit above the convex curve.
• Simpson: $(f_0 + 4f_1 + f_2)/6 = (0 + 1 + 1)/6 = 1/3 \approx 0.3333$. Exact: parabolas reproduce quadratics, so the only error is rounding.

Same three function values, 0.375 versus 0.3333 against truth $1/3$. The weights did all the work. Gaussian quadrature with $n$ points integrates degree-$(2n-1)$ polynomials exactly; Monte Carlo's $1/\sqrt{N}$ rate ignores dimension, which is why hundred-dimensional expectations sample rather than grid.

Worked example

Integrate $f(x) = x^2$ on $[0, 1]$ (truth $\int_0^1 x^2 dx = 1/3$). The trapezoid computation above gives $0.375$, off by $0.0417$. Halve $h$ to $0.25$ (points $0, 0.25, 0.5, 0.75, 1$): $0.25 \cdot (0/2 + 0.0625 + 0.25 + 0.5625 + 1/2) = 0.25 \cdot 1.375 = 0.34375$. Error falls $0.0417 \to 0.0104$, a factor of 4 for halved $h$: the $O(h^2)$ rate made visible. Simpson needed no refinement at all.

Code

The point is the check: implement both rules and confirm the error ratio on refinement.

1def trapezoid(f: list[float], h: float) -> float:2    return h * (f[0] / 2 + sum(f[1:-1]) + f[-1] / 2)3
4
5def simpson(f0: float, f1: float, f2: float, h: float) -> float:6    return h / 3 * (f0 + 4 * f1 + f2)  # one double panel7
8
9vals = [0.0, 0.25, 1.0]  # x squared at 0, 0.5, 110print(trapezoid(vals, 0.5))  # -> 0.37511print(round(simpson(0.0, 0.25, 1.0, 0.5), 4))  # -> 0.3333

Watch Out For

The Quick Version

• Quadrature replaces integrals with weighted sums of function values at chosen points.
• Trapezoid on $x^2$ gives 0.375; Simpson gives exactly $1/3$ from the same three values.
• Halving trapezoid panels quarters the error: the $O(h^2)$ rate in action.
• Gaussian quadrature maximises exactness per evaluation; Monte Carlo ignores dimension at $1/\sqrt{N}$.
• Composite low-order rules beat single high-order ones; grids die in high dimensions.

What to Read Next

• Integration and Area Under a Curve
• Monte Carlo Methods
• Expectation, Variance and Covariance