xs <- seq(0, 1, length.out = 200)
draw_line(
x = xs,
y = xs^2,
area = TRUE,
points = FALSE,
xlab = "x",
ylab = "f(x)"
)14 Integration
Differentiation takes a function apart to find its rate of change. Integration puts it back together: given a rate, recover the total.
For data science the payoff is probability. A continuous random variable has a density rather than a list of probabilities, and every question you can ask about it — what is the chance of landing in this range, what is the mean, what is the variance — is an integral (Section 14.11).
14.1 Area under a curve
The definite integral of \(f\) from \(a\) to \(b\),
\[ \int_a^b f(x)\,dx \]
is the signed area between the curve and the \(x\)-axis. Signed: area below the axis counts as negative.
The pieces of the notation each mean something. \(\int\) is an elongated S for “sum”; \(f(x)\,dx\) is the area of a sliver of width \(dx\) and height \(f(x)\); \(a\) and \(b\) are the limits of integration. Read the whole thing as sum up \(f(x)\,dx\) from \(a\) to \(b\).
14.2 Riemann sums
How do you compute an area under a curve? Approximate it with rectangles, then let them get thin.
Divide \([a,b]\) into \(n\) pieces of width \(\Delta x = (b-a)/n\), pick a point \(x_i^*\) in each, and add up the rectangles:
\[ \sum_{i=1}^{n} f(x_i^*)\,\Delta x \]
That is a Riemann sum. The integral is its limit:
\[ \int_a^b f(x)\,dx = \lim_{n\to\infty}\sum_{i=1}^{n} f(x_i^*)\,\Delta x \tag{14.1}\]
Compare Equation 14.1 with Equation 11.1. Both define an operation as a limit of something finite and computable, and in both the limit is what makes it exact.
n <- 8
edges <- seq(0, 1, length.out = n + 1)
mids <- (edges[-1] + edges[-(n + 1)]) / 2
hts <- mids^2
sx <- c(0, as.vector(rbind(edges[-(n + 1)], edges[-1])), 1)
sy <- c(0, rep(hts, each = 2), 0)
draw_plane(
curves = list(
`x^2` = cbind(xs, xs^2),
rectangles = cbind(sx, sy)
),
curve_color = c(amds_colors[1], amds_gray),
xlab = "x", ylab = "f(x)"
)riemann <- function(f, a, b, n, rule = "mid") {
edges <- seq(a, b, length.out = n + 1)
pts <- switch(rule,
left = edges[-(n + 1)],
right = edges[-1],
mid = (edges[-1] + edges[-(n + 1)]) / 2
)
sum(f(pts)) * (b - a) / n
}
sq <- function(x) x^2
sapply(c(4, 8, 50, 1000), function(n) riemann(sq, 0, 1, n))[1] 0.3281250 0.3320312 0.3333000 0.3333333Converging on \(1/3\). The left and right rules bracket the answer, since \(x^2\) is increasing:
14.3 The antiderivative
An antiderivative of \(f\) is any function \(F\) with \(F' = f\). Differentiation run backwards.
\[ \int f(x)\,dx = F(x) + C \]
This is the indefinite integral — no limits, and the answer is a family of functions. The constant \(C\) is there because differentiating kills constants, so \(x^3/3\) and \(x^3/3 + 7\) are both antiderivatives of \(x^2\).
Reversing the table from Section 11.4:
| \(f(x)\) | \(\int f(x)\,dx\) |
|---|---|
| \(x^n\), \(n \neq -1\) | \(\dfrac{x^{n+1}}{n+1} + C\) |
| \(1/x\) | \(\log\lvert x\rvert + C\) |
| \(e^x\) | \(e^x + C\) |
| \(e^{ax}\) | \(\dfrac{1}{a}e^{ax} + C\) |
| \(\cos x\) | \(\sin x + C\) |
| \(\sin x\) | \(-\cos x + C\) |
The \(n \neq -1\) exclusion is not a technicality — the power rule would divide by zero, and the missing case is exactly the one that produces the logarithm.
WarningWatch out
Not every function has an antiderivative you can write down. \(e^{-x^2}\) is the famous case: it is perfectly integrable, but no combination of elementary functions differentiates to it.
This is why the normal distribution has no closed-form CDF and why pnorm() is a numerical routine rather than a formula. Being unable to find \(F\) is normal, and it is what Section 14.10 exists for.
14.4 The fundamental theorem of calculus
Two operations defined completely differently — one a limit of difference quotients, the other a limit of sums of rectangles — turn out to be inverses.
\[ \int_a^b f(x)\,dx = F(b) - F(a) \qquad\text{where } F' = f \tag{14.2}\]
This is extraordinary and easy to take for granted. It says you can compute an area without ever summing a rectangle: find any antiderivative, evaluate at both ends, subtract. Note the \(C\) cancels, which is why the indefinite integral’s ambiguity does no harm here.
Worked example. \(\int_0^1 x^2\,dx\). An antiderivative is \(F(x) = x^3/3\), so
\[ \int_0^1 x^2\,dx = \frac{1^3}{3} - \frac{0^3}{3} = \frac{1}{3} \]
matching the Riemann sums above. Another: \(\int_0^\pi \sin x\,dx = [-\cos x]_0^\pi = 1 + 1 = 2\).
exact riemann sine
0.3333333 0.3333333 2.0000000 The bracket notation \([F(x)]_a^b\) is standard shorthand for \(F(b) - F(a)\).
14.5 Techniques of integration
Differentiation is mechanical: apply the rules and you are done. Integration is not — it is pattern recognition, and two patterns cover most cases.
Substitution reverses the chain rule. If the integrand contains a function and its own derivative, substitute \(u = g(x)\), \(du = g'(x)\,dx\):
\[ \int f(g(x))\,g'(x)\,dx = \int f(u)\,du \]
Worked example. \(\int_0^1 2x\,e^{x^2}dx\). Take \(u = x^2\), so \(du = 2x\,dx\):
\[ \int_{u=0}^{u=1} e^u\,du = e^1 - e^0 = e - 1 \approx 1.718 \]
Integration by parts reverses the product rule:
\[ \int u\,dv = uv - \int v\,du \]
Use it when the integrand is a product where one factor gets simpler on differentiation.
Worked example. \(\int_0^1 x e^x dx\). Take \(u = x\) (which differentiates to 1) and \(dv = e^x dx\):
\[ \int_0^1 xe^x dx = [xe^x]_0^1 - \int_0^1 e^x dx = e - (e - 1) = 1 \]
14.6 Definite integrals
Properties worth knowing, all following from Equation 14.2:
\[ \begin{aligned} \int_a^b (f + g) &= \int_a^b f + \int_a^b g &&\text{linearity} \\ \int_a^b cf &= c\int_a^b f &&\text{constants factor out} \\ \int_a^b f &= -\int_b^a f &&\text{swapping limits flips the sign} \\ \int_a^b f &= \int_a^c f + \int_c^b f &&\text{splitting the interval} \end{aligned} \]
Linearity is the one that matters most: it is why the expectation of a sum is the sum of expectations, no independence required.
14.7 Improper integrals
An integral is improper when a limit is infinite, or the integrand blows up inside the range. Define it as a limit:
\[ \int_1^\infty f(x)\,dx = \lim_{b\to\infty}\int_1^b f(x)\,dx \]
If the limit is finite the integral converges; otherwise it diverges.
Two functions that look similar and behave completely differently:
\[ \int_1^\infty \frac{1}{x^2}dx = \left[-\frac{1}{x}\right]_1^\infty = 1 \qquad \int_1^\infty \frac{1}{x}dx = [\log x]_1^\infty = \infty \]
[,1] [,2] [,3]
one_over_x2 0.900000 0.99000 0.99990
one_over_x 2.302585 4.60517 9.21034The first settles down to 1; the second grows without bound — slowly, but forever. Both curves go to zero, so going to zero is not enough: what matters is how fast.
draw_fun(
list(
`1/x` = function(x) 1 / x,
`1/x^2` = function(x) 1 / x^2
),
from = 1, to = 8,
ylab = "f(x)"
)
NoteIn machine learning
Convergence of an improper integral is exactly the condition for a probability distribution to exist. A density must satisfy \(\int_{-\infty}^{\infty} p(x)\,dx = 1\), so it has to decay fast enough. Heavy-tailed distributions decay slowly, which is why some of them have no finite mean or variance — the defining integral diverges (Section 19.5).
14.8 Multiple integrals
For \(f(x,y)\), integrating over a region gives volume under the surface. A double integral over a rectangle evaluates as two single integrals, innermost first:
\[ \int_c^d\!\!\int_a^b f(x,y)\,dx\,dy \]
Worked example. \(\int_0^1\!\!\int_0^1 xy\,dx\,dy\). Inner integral, treating \(y\) as constant:
\[ \int_0^1 xy\,dx = y\left[\frac{x^2}{2}\right]_0^1 = \frac{y}{2} \]
Then the outer:
\[ \int_0^1 \frac{y}{2}dy = \left[\frac{y^2}{4}\right]_0^1 = \frac{1}{4} \]
[1] 0.25Over a rectangle with a separable integrand the order does not matter. Over an awkward region it very much does, and choosing well is most of the work.
14.9 Change of variables
Substitution in several dimensions. Changing coordinates stretches area, and the stretch factor is the absolute determinant of the Jacobian (Section 12.7):
\[ \int_R f(\mathbf{x})\,d\mathbf{x} = \int_{R'} f(\mathbf{g}(\mathbf{u}))\,\lvert\det\mathbf{J}_{\mathbf{g}}\rvert\,d\mathbf{u} \tag{14.3}\]
This is Section 5.13 cashing in: the determinant was defined as an area scaling factor, and here it is, scaling areas.
The standard case is polar coordinates, \(x = r\cos\theta\), \(y = r\sin\theta\), where \(\det\mathbf{J} = r\) — so \(dx\,dy\) becomes \(r\,dr\,d\theta\). That extra \(r\) is what makes the Gaussian integral tractable:
\[ \int_{-\infty}^{\infty} e^{-x^2/2}dx = \sqrt{2\pi} \tag{14.4}\]
There is no elementary antiderivative for \(e^{-x^2/2}\), but squaring the integral turns it into a double integral over the plane, and in polar coordinates the stray \(r\) from the Jacobian is exactly what a substitution needs.
numeric sqrt_2pi
2.506628 2.506628 That \(\sqrt{2\pi}\) is the normalizing constant in the normal density, and it comes from here.
14.10 Numerical integration
Since most integrands have no elementary antiderivative, numerical quadrature is the normal case, not the fallback.
Riemann sums work but converge slowly. Better rules approximate the curve rather than stepping it:
| Rule | Fits each panel with | Error |
|---|---|---|
| Midpoint | a rectangle | \(O(h^2)\) |
| Trapezoid | a straight line | \(O(h^2)\) |
| Simpson | a parabola | \(O(h^4)\) |
trapezoid <- function(f, a, b, n) {
x <- seq(a, b, length.out = n + 1)
y <- f(x)
(b - a) / n * (sum(y) - (y[1] + y[n + 1]) / 2)
}
simpson <- function(f, a, b, n) {
# Panels are paired into parabolas, so n must be even.
# Without this the weights come out wrong and the answer
# is plausible but incorrect — the worst kind of failure.
stopifnot(n %% 2 == 0)
x <- seq(a, b, length.out = n + 1)
y <- f(x)
# Weights run 1, 4, 2, 4, ..., 4, 1: the interior nodes
# alternate, and n being even is what ends them on a 4.
w <- rep(c(2, 4), length.out = n + 1)
w[c(1, n + 1)] <- 1
(b - a) / (3 * n) * sum(w * y)
}
nq <- c(4, 8, 16)
rbind(
trapezoid = sapply(nq, trapezoid, f = sq, a = 0, b = 1),
simpson = sapply(nq, simpson, f = sq, a = 0, b = 1)
) [,1] [,2] [,3]
trapezoid 0.3437500 0.3359375 0.3339844
simpson 0.3333333 0.3333333 0.3333333Simpson’s rule is exact here at every \(n\), because it fits each panel with a parabola and \(x^2\) is a parabola. Nothing is being approximated.
Notice that \(n\) must be even. Simpson’s rule works by taking panels in pairs and passing one parabola through each trio of points, so an odd count leaves a panel with no partner. That pairing is also where the accuracy comes from: fitting a parabola rather than a line is what buys \(O(h^4)\) instead of \(O(h^2)\), and it is why Simpson integrates cubics exactly as well (exercise 4).
The stopifnot() in the code above enforces it. That is not defensive clutter — with an odd \(n\) the weight pattern silently comes out wrong and the function returns a plausible number rather than an error, which is precisely the failure mode this book keeps warning about.
In practice use R’s adaptive routine, which subdivides where the integrand is difficult and returns an error estimate:
value abs_error
3.333333e-01 3.700743e-15
WarningWatch out
Quadrature rules like these are excellent in one or two dimensions and hopeless beyond about five: a grid with \(k\) points per axis needs \(k^d\) evaluations, so the cost is exponential in dimension.
This is why high-dimensional integrals — which is to say, most Bayesian computation — use Monte Carlo instead. Its error falls as \(1/\sqrt{N}\) regardless of dimension. That rate is poor in one dimension and unbeatable in fifty (Section 22.4).
14.11 Integrals in probability
For a continuous random variable with density \(p(x)\), every quantity of interest is an integral.
\[ \begin{aligned} \int_{-\infty}^{\infty} p(x)\,dx &= 1 &&\text{total probability} \\ P(a \leq X \leq b) &= \int_a^b p(x)\,dx &&\text{probability of a range} \\ \mathbb{E}[X] &= \int_{-\infty}^{\infty} x\,p(x)\,dx &&\text{mean} \\ \mathbb{E}[g(X)] &= \int_{-\infty}^{\infty} g(x)\,p(x)\,dx &&\text{any function} \end{aligned} \]
Each mirrors a sum from the discrete case, with \(\int\) in place of \(\sum\) and a density in place of a probability.
The standard normal density is
\[ p(x) = \frac{1}{\sqrt{2\pi}}e^{-x^2/2} \]
whose \(\sqrt{2\pi}\) is precisely Equation 14.4, dividing through so the total comes to 1.
total within_1sd within_2sd mean
1.0000000 0.6826895 0.9544997 0.0000000 The familiar 68% and 95% are these integrals, and the mean is zero by symmetry — the positive and negative halves cancel exactly.
14.12 Summary
| Idea | Statement |
|---|---|
| Definite integral | signed area under the curve |
| Riemann sum | \(\sum f(x_i^*)\Delta x\), exact in the limit |
| Antiderivative | any \(F\) with \(F' = f\) |
| Fundamental theorem | \(\int_a^b f = F(b) - F(a)\) |
| Substitution | reverses the chain rule |
| By parts | \(\int u\,dv = uv - \int v\,du\) |
| Change of variables | multiply by \(\lvert\det\mathbf{J}\rvert\) |
| In R | integrate(f, a, b) |
14.13 Exercises
1. Compute \(\int_0^2 (3x^2 + 1)\,dx\) by hand and check with integrate().
TipSolution
2. Use substitution for \(\int_0^1 x\,e^{x^2}dx\). How does it relate to the worked example?
TipSolution
With \(u = x^2\), \(du = 2x\,dx\), so \(x\,dx = du/2\):
\[ \int_0^1 xe^{x^2}dx = \frac{1}{2}\int_0^1 e^u du = \frac{e-1}{2} \approx 0.859 \]
Exactly half the worked example, because the integrand is half of it. Linearity, not a coincidence.
3. Does \(\int_1^\infty x^{-1.5}dx\) converge? What about \(\int_1^\infty x^{-0.5}dx\)?
TipSolution
For \(\int_1^\infty x^{-p}dx\) the antiderivative is \(x^{1-p}/(1-p)\), which stays finite as \(x \to \infty\) exactly when \(p > 1\).
So \(p = 1.5\) converges, to \(1/0.5 = 2\). And \(p = 0.5\) diverges.
The threshold is \(p = 1\), which is the \(1/x\) case — the boundary between the two, and divergent itself.
4. Verify that Simpson’s rule is exact for any cubic, using \(\int_0^2 x^3 dx = 4\).
TipSolution
cube <- function(x) x^3
c(
exact = 4,
simpson_n2 = simpson(cube, 0, 2, 2),
simpson_n4 = simpson(cube, 0, 2, 4),
trapezoid_n4 = trapezoid(cube, 0, 2, 4)
) exact simpson_n2 simpson_n4 trapezoid_n4
4.00 4.00 4.00 4.25 Simpson is exact with just two panels, while the trapezoid rule is not. Fitting each panel with a parabola gets cubics right too — the error terms cancel by symmetry, which is why Simpson is \(O(h^4)\) rather than the \(O(h^3)\) you might expect.
5. For the standard normal, compute \(P(-1.96 \leq X \leq 1.96)\) by integration and compare with pnorm().
TipSolution
integral pnorm
0.9500042 0.9500042 The familiar 95%. pnorm() is this integral, computed by a specialized routine — recall from Section 14.3 that no elementary antiderivative exists, so there is no formula to evaluate instead.
6. Compute \(\mathbb{E}[X^2]\) for the standard normal by integration. What does it tell you about the variance?
TipSolution
E_X2 E_X variance
1 0 1 \(\mathbb{E}[X^2] = 1\) and \(\mathbb{E}[X] = 0\), so the variance \(\mathbb{E}[X^2] - (\mathbb{E}[X])^2 = 1\) — as the name standard normal promises.
Note that \(\mathbb{E}[X^2] \neq (\mathbb{E}[X])^2\) here, and the gap between them is the variance. That identity gets used constantly in Section 19.6.