Integrals, expectations and the change of variables

Before this lesson: lesson 01 (derivatives in one variable), lesson 03 (the Jacobian $Df$ and the linear approximation), linear algebra lesson 03 (the determinant as a volume factor). It helps to know random variables and densities from the probability and statistics track; we repeat what we need. Time: about 85 minutes, plus the notebook.

---

Why this matters #

Expectations are everywhere in ML. The training loss is an average that estimates an expected loss. The evidence lower bound (ELBO) of a variational autoencoder, a policy's expected return in reinforcement learning and the marginal likelihood of a Bayesian model are all integrals. Most of them have no closed form, so we estimate them with Monte Carlo averages, and we need to know how fast these averages converge.

The change-of-variables rule is the other main tool. It gives the density of a transformed random variable. This is the basis of normalizing flows (exact likelihoods for generative models), of the likelihood of continuous flows and of the probability-flow ODE of diffusion models, and of sampling by the inverse CDF, by Box–Muller and by reparameterization. In science, many models predict a log-quantity (log-precipitation, log-permeability, the log of a ground-motion amplitude). To report a likelihood for the quantity itself, you need the Jacobian factor of the logarithm.

1. Predict first #

2. The picture #

An integral is a limit of sums #

Cut the interval $[a, b]$ into $n$ strips of width $h = (b - a)/n$. On each strip, put a rectangle whose height is the value of $f$ at one point of the strip. The sum of the rectangle areas is a Riemann sum. When $n$ grows, the sum approaches the area under the curve: the integral $\int_a^b f(x)\,dx$. Area below the $x$-axis counts as negative.

The choice of the point changes how fast the sum converges. The left sum uses the left end of each strip. The midpoint rule uses the middle. The trapezoid rule joins the two end values by a straight line. A completely different idea is Monte Carlo: pick $n$ random points $u_1, \dots, u_n$, uniform on $[a, b]$, and use $(b - a)$ times the average of $f(u_i)$.

Interactive figure: loading…

With $n = 4$, the midpoint rule for $\int_0^1 x^2\,dx = 1/3$ uses the midpoints $0.125, 0.375, 0.625, 0.875$:

so the error is $1/3 - 0.328125 = 0.0052$.

A smooth map scales small areas by $|\det DT|$ #

Now go to two dimensions. A map $T$ sends a point $y$ to a point $x = T(y)$. On a small square around $y_0$, a smooth map is almost linear: $T(y_0 + h) \approx T(y_0) + DT(y_0)\,h$, where $DT(y_0)$ is the Jacobian matrix (lesson 03). A linear map multiplies every area by $|\det|$ of its matrix (linear algebra lesson 03). So the small square goes to a small curved region with area about $|\det DT(y_0)|$ times the area of the square.

Interactive figure: loading…

This is why an integral in new coordinates needs the factor $|\det DT|$. Each small box of volume $dy$ in the new coordinates covers a region of volume $|\det DT(y)|\,dy$ in the old ones.

Probability moves, so the density must change #

A random variable $X$ with a density $p_X$ puts probability $p_X(x)\,dx$ on a small interval of width $dx$ around $x$. Take the simplest case: $X$ uniform on $(0, 1)$ and $Y = 2X$. Then $Y$ is uniform on $(0, 2)$. The interval got twice as long, and the total probability is still 1, so the density must be half as high: $p_Y = 1/2$ on $(0, 2)$.

For a general smooth increasing map $y = g(x)$, a small interval $[x, x + dx]$ goes to $[y, y + dy]$ with $dy = g'(x)\,dx$. Both intervals hold the same probability:

where $x = g^{-1}(y)$ is the point that $g$ sends to $y$.

Where $g$ stretches space ($|g'| > 1$), the density goes down. Where $g$ compresses space ($|g'| < 1$), the density goes up.

Interactive figure: loading…

3. Definitions and main results #

3.1 The definite integral and the fundamental theorem #

Let $f$ be a function on $[a, b]$. Cut $[a, b]$ into $n$ pieces $[x_{i-1}, x_i]$ of width $h = (b - a)/n$, with $x_i = a + ih$, and choose a point $\xi_i$ in each piece. The Riemann sum is

$f$ is (Riemann) integrable if $S_n$ tends to the same limit for every choice of the points $\xi_i$ as $n \to \infty$. This limit is the definite integral $\int_a^b f(x)\,dx$. Every continuous function on $[a, b]$ is integrable.

How fast the rules converge. Let $f$ have a continuous second derivative. Then, for some points $\xi \in [a, b]$:

Rule	Height of strip $i$	Error $I - S_n$	Order
left sum	$f(x_{i-1})$	$\frac h2\big(f(b) - f(a)\big) + O(h^2)$	$O(h)$
midpoint	$f\big(\frac{x_{i-1} + x_i}{2}\big)$	$\frac{(b - a)h^2}{24} f''(\xi)$	$O(h^2)$
trapezoid	$\frac12\big(f(x_{i-1}) + f(x_i)\big)$	$-\frac{(b - a)h^2}{12} f''(\xi)$	$O(h^2)$

For $f(x) = x^2$ on $[0, 1]$, $f'' = 2$, so the midpoint error is exactly $h^2/12$; with $h = 1/4$ this is $0.0052$, as in section 2. The midpoint and trapezoid errors have opposite signs, and the midpoint error is about half as large. (If $f(a) = f(b)$, the left sum equals the trapezoid sum, so it is also $O(h^2)$.)

Theorem (fundamental theorem of calculus, FTC). Let $f$ be continuous on $[a, b]$.

1. The function $F(x) = \int_a^x f(t)\,dt$ is differentiable, and $F'(x) = f(x)$.
2. If $G$ is any antiderivative of $f$ (that is, $G' = f$ on $[a, b]$), then $\int_a^b f(x)\,dx = G(b) - G(a)$.

Proof.

1. For $h > 0$, $F(x + h) - F(x) = \int_x^{x + h} f(t)\,dt$. On the small interval $[x, x + h]$, $f$ lies between its minimum $m_h$ and its maximum $M_h$, so $m_h\,h \le F(x + h) - F(x) \le M_h\,h$.
2. Divide by $h$. Because $f$ is continuous, $m_h$ and $M_h$ both tend to $f(x)$ as $h \to 0$. So $(F(x + h) - F(x))/h \to f(x)$. The case $h < 0$ is the same.
3. For part 2: $(F - G)' = f - f = 0$, so $F - G$ is a constant $C$ (mean value theorem, lesson 01).
4. Then $G(b) - G(a) = F(b) - F(a) = \int_a^b f - 0$. $\blacksquare$

Substitution (the 1-D change of variables). Let $g$ have a continuous derivative on $[c, d]$, and let $f$ be continuous. Then

Proof. Let $G$ be an antiderivative of $f$. By the chain rule (lesson 04), $\frac{d}{dy}G(g(y)) = f(g(y))\,g'(y)$. Apply the FTC to both sides: each side equals $G(g(d)) - G(g(c))$. $\blacksquare$

If $g$ is one-to-one, write the image interval as $I = g([c, d])$. When $g$ is decreasing, $g(d) < g(c)$: the limits come out swapped, and swapping them back changes the sign. Both cases give one formula with an absolute value:

This is the form that generalizes to $d$ dimensions. Example: with $x = g(y) = y^2$, $\int_0^1 2y\,e^{y^2}\,dy = \int_0^1 e^x\,dx = e - 1 = 1.7183$.

Integration by parts. From the product rule $(uv)' = u'v + uv'$ and the FTC:

Lessons 12 and 13 use it all the time. One example that appears in ML is Stein's lemma. Let $X \sim N(0, 1)$ with density $\varphi(x) = e^{-x^2/2}/\sqrt{2\pi}$, so $\varphi'(x) = -x\,\varphi(x)$. Take $u = f$ and $v = -\varphi$. If $f\varphi \to 0$ at $\pm\infty$, the boundary term vanishes, and

For $f = \sin$: $\mathbb E[X\sin X] = \mathbb E[\cos X] = e^{-1/2} = 0.6065$.

3.2 Expectations are integrals #

A random variable $X$ in $\R^d$ has a density $p_X$ if $p_X \ge 0$, $\int p_X(x)\,dx = 1$, and $P(X \in B) = \int_B p_X(x)\,dx$ for every (measurable) set $B$. A density is not a probability: it is probability per unit volume. So it can be larger than 1. The uniform density on $(0, 1/2)$ is 2.

The expectation of $f(X)$ is the integral

defined when $\int |f(x)|\,p_X(x)\,dx < \infty$. You do not need the density of $f(X)$ for this; the density of $X$ is enough. A 1-D integral is also an expectation: $\int_a^b f(x)\,dx = (b - a)\,\mathbb E[f(U)]$ with $U$ uniform on $[a, b]$. This is how the Monte Carlo mode of the figure works.

Monte Carlo. Draw $x_1, \dots, x_N$ independently from $p_X$, and take the average

1. Unbiased. $\mathbb E[\hat\mu_N] = \frac1N \sum_i \mathbb E[f(x_i)] = \mu$, where $\mu = \mathbb E[f(X)]$.
2. Variance. The terms are independent, so their variances add: $\operatorname{Var}(\hat\mu_N) = \frac{1}{N^2}\cdot N\sigma^2 = \sigma^2/N$, where $\sigma^2 = \operatorname{Var} f(X)$.
3. Standard error. The typical error is $\sigma/\sqrt N$. By the central limit theorem, $\hat\mu_N$ is within $\pm 2\sigma/\sqrt N$ of $\mu$ in about 95% of runs. In practice you estimate $\sigma$ by the sample standard deviation.

The rate $1/\sqrt N$ does not depend on the dimension $d$. A grid rule with $n$ points in $d$ dimensions has only $n^{1/d}$ points per axis, so the midpoint error is $O(n^{-2/d})$. For $d > 4$, Monte Carlo wins. This is why high-dimensional expectations in ML are estimated by sampling. Variance reduction and importance sampling are in the probability and statistics track.

Example. $X \sim N(0, 1)$ and $f(x) = x^2$: $\mu = 1$ and $\sigma^2 = \mathbb E[X^4] - 1 = 3 - 1 = 2$. With $N = 10\,000$ samples, the standard error is $\sqrt{2/10\,000} = 0.0141$.

3.3 Multiple integrals and Fubini's theorem #

A double integral $\int_R f(x, y)\,d(x, y)$ over a rectangle $R$ is the limit of sums over small rectangles of area $\Delta x\,\Delta y$, in the same way as in 1-D. To compute it, we integrate one variable at a time: an iterated integral. The theorems below say when the order does not matter. Write a point of $\R^{m+k}$ as $(x, y)$ with $x \in \R^m$, $y \in \R^k$.

• Tonelli's theorem. If $f \ge 0$, then $\int f(x, y)\,d(x, y) = \int \big(\int f(x, y)\,dy\big)\,dx = \int \big(\int f(x, y)\,dx\big)\,dy$. All three may be $+\infty$, but they are always equal.
• Fubini's theorem. The same equalities hold for every $f$ with $\int |f(x, y)|\,d(x, y) < \infty$.
• Product case. $\int f(x)\,g(y)\,d(x, y) = \big(\int f(x)\,dx\big)\big(\int g(y)\,dy\big)$.

Example: $\int_0^1\!\int_0^2 x y^2\,dy\,dx = \int_0^1 x \cdot \tfrac83\,dx = \tfrac43$, and in the other order $\int_0^2 y^2 \cdot \tfrac12\,dy = \tfrac43$. In probability, Fubini's theorem lets you compute an expectation over $(X, Y)$ one variable at a time, and swap a sum and an expectation.

When the order matters: a counterexample

Take $f(x, y) = \dfrac{x^2 - y^2}{(x^2 + y^2)^2}$ on $(0, 1)^2$.

1. Since $\frac{\partial}{\partial y}\frac{y}{x^2 + y^2} = \frac{x^2 - y^2}{(x^2 + y^2)^2}$, the inner integral is $\int_0^1 f(x, y)\,dy = \frac{1}{1 + x^2}$.
2. So $\int_0^1\!\int_0^1 f\,dy\,dx = \int_0^1 \frac{dx}{1 + x^2} = \arctan 1 = \frac\pi4$.
3. But $f(y, x) = -f(x, y)$, so the other order gives $-\frac\pi4$.

The two orders disagree. Fubini does not apply here, because $\int |f| = \infty$. In polar coordinates $f = \cos(2\theta)/r^2$, so near the corner $(0, 0)$ the integral of $|f|$ contains $\int_0 \frac{1}{r^2}\,r\,dr = \int_0 \frac{dr}{r} = \infty$.

3.4 Polar coordinates #

Polar coordinates describe a point of the plane by its distance $r$ from the origin and its angle $\theta$: $T(r, \theta) = (r\cos\theta,\ r\sin\theta)$. The Jacobian is

so $\det DT = r\cos^2\theta + r\sin^2\theta = r$.

So the area element is $dx\,dy = r\,dr\,d\theta$. The picture: a small polar cell has a radial side $\Delta r$ and an arc side $r\,\Delta\theta$ (an arc of angle $\Delta\theta$ on a circle of radius $r$). It is almost a rectangle, so its area is about $r\,\Delta r\,\Delta\theta$. Cells far from the origin are bigger.

Interactive figure: loading…

For an integrable $f$ on $\R^2$:

(The origin and the ray $\theta = 0$ are missing from the image of $(0, \infty) \times (0, 2\pi)$, but they have area 0, so they do not change the integral.) Worked example (a) uses this formula, and lesson 11 uses it to compute $\int e^{-x^2}dx = \sqrt\pi$.

3.5 Change of variables in $d$ dimensions #

A diffeomorphism $T: V \to U$ between open sets of $\R^d$ is a one-to-one map onto $U$ such that $T$ and $T^{-1}$ both have continuous derivatives. For a one-to-one map with a continuous derivative, this holds exactly when $\det DT(y) \ne 0$ everywhere (inverse function theorem). Then $D(T^{-1})(x) = [DT(y)]^{-1}$ at $x = T(y)$.

Theorem (change of variables). Let $T: V \to U$ be a diffeomorphism and $f$ integrable on $U$. Then

In words: substitute $x = T(y)$, and replace $dx$ by $|\det DT(y)|\,dy$. In 1-D this is the substitution rule with $|g'|$. For polar coordinates it gives $r\,dr\,d\theta$.

Why it is true (proof sketch)

1. Linear maps. A linear map $A$ sends a box $Q$ to a parallelepiped of volume $|\det A|\,\mathrm{vol}(Q)$ (linear algebra lesson 03).
2. Small cubes. Let $Q$ be a cube of side $\varepsilon$ centred at $y$. By the linear approximation (lesson 03), $T(y + h) = T(y) + DT(y)\,h + o(\|h\|)$. So $T(Q)$ differs from the parallelepiped $T(y) + DT(y)(Q - y)$ only by $o(\varepsilon)$ near its boundary, and $\mathrm{vol}(T(Q)) = |\det DT(y)|\,\varepsilon^d + o(\varepsilon^d)$.
3. Sum over cubes. Cut $V$ into small cubes $Q_i$ with centres $y_i$. Their images $T(Q_i)$ cut $U$ into small pieces, so $\int_U f\,dx = \sum_i \int_{T(Q_i)} f\,dx \approx \sum_i f(T(y_i))\,|\det DT(y_i)|\,\varepsilon^d$.
4. The right side is a Riemann sum of $f(T(y))\,|\det DT(y)|$ over $V$. As $\varepsilon \to 0$ it tends to $\int_V f(T(y))\,|\det DT(y)|\,dy$.

A full proof must show that the errors in steps 2–3 are small uniformly over all cubes. See Spivak, Calculus on Manifolds, theorem 3-13, or Munkres, Analysis on Manifolds, §17–19.

3.6 The density of a transformed random variable #

Theorem (change of variables for densities). Let $X$ have the density $p_X$ on an open set $U \subseteq \R^d$, let $g: U \to V$ be a diffeomorphism, and let $Y = g(X)$. Then for $y \in V$, with $x = g^{-1}(y)$:

and in log form

Read it as conservation of probability: "density times volume" is the same before and after the map, $p_Y(y)\,|dy| = p_X(x)\,|dx|$.

Proof.

1. 1-D, $g$ increasing. The cumulative distribution function (CDF) of $Y$ is $F_Y(y) = P(g(X) \le y) = P(X \le g^{-1}(y)) = F_X(g^{-1}(y))$. Differentiate with the chain rule: $p_Y(y) = p_X(g^{-1}(y))\,(g^{-1})'(y)$, and $(g^{-1})' > 0$.
2. 1-D, $g$ decreasing. Now $g(X) \le y$ exactly when $X \ge g^{-1}(y)$, so $F_Y(y) = 1 - F_X(g^{-1}(y))$ and $p_Y(y) = -p_X(g^{-1}(y))\,(g^{-1})'(y)$, with $(g^{-1})' < 0$. Both cases give $p_Y(y) = p_X(g^{-1}(y))\,|(g^{-1})'(y)|$.
3. $d$ dimensions. A CDF is not convenient in $\R^d$, so we use sets. For every set $B \subseteq V$: $P(Y \in B) = P(X \in g^{-1}(B)) = \int_{g^{-1}(B)} p_X(x)\,dx$.
4. Change variables with $T = g^{-1}$, which maps $B$ onto $g^{-1}(B)$ (section 3.5): $P(Y \in B) = \int_B p_X(g^{-1}(y))\,|\det Dg^{-1}(y)|\,dy$.
5. This holds for every $B$. A density is fixed by its integrals over all sets (up to a set of volume zero), so the integrand is $p_Y$.
6. By the inverse function theorem, $Dg^{-1}(y) = [Dg(x)]^{-1}$, and $\det(M^{-1}) = 1/\det M$. Take logs. $\blacksquare$

Which way does the Jacobian go? The two maps go in opposite directions. In the integral formula of section 3.5, $T$ maps the new variable $y$ to the old variable $x$. Here $g$ maps $X$ to $Y$, and the proof uses the integral formula with $T = g^{-1}$. A quick check: for $Y = 2X$ the density must go down, so you divide by $|g'| = 2$.

Maps that are not one-to-one. Suppose $U$ splits into pieces $U_1, \dots, U_K$ (up to a set of volume zero), and $g$ is a diffeomorphism from each $U_k$ onto its image. Each piece is a branch. Then $P(Y \in B) = \sum_k P(X \in U_k,\ g(X) \in B)$. Apply the theorem on each piece:

Example: $Y = X^2$ with $X \sim N(0, 1)$. The pieces are $x < 0$ and $x > 0$. The preimages of $y > 0$ are $\pm\sqrt y$, and $|g'(x)| = 2|x| = 2\sqrt y$. So

the $\chi^2_1$ density. With one branch only, you get half of it, and it integrates to $1/2$.

3.7 Compositions, triangular Jacobians and coupling layers #

Let $g = g_K \circ \dots \circ g_1$, with $z_0 = x$, $z_k = g_k(z_{k-1})$ and $y = z_K$. By the chain rule, $Dg(x) = Dg_K(z_{K-1}) \cdots Dg_1(z_0)$, and the determinant of a product is the product of the determinants. So the log-determinants add:

• Density evaluation (training by maximum likelihood): start from $y$, invert the layers in reverse order ($z_{K-1} = g_K^{-1}(y)$, and so on), and add up the log-determinants on the way.
• Sampling: draw $z_0 \sim p_X$ and apply $g_1, \dots, g_K$.
• You never form the Jacobian of the whole composition. Each layer needs a cheap inverse (for density evaluation) and a cheap $\log|\det|$.

Why "cheap" matters. A dense $d \times d$ determinant costs $O(d^3)$ operations (an LU factorization). For a $64 \times 64$ RGB image, $d = 12\,288$, for every layer and every sample. This is far too slow inside a training loop.

Triangular Jacobians. The determinant of a triangular matrix is the product of its diagonal entries. So if $Dg$ is triangular, $\log|\det Dg| = \sum_i \log|\partial g_i/\partial x_i|$: only $O(d)$ work.

Affine coupling layer (RealNVP). Split $x = (x_A, x_B)$, where $x_A$ holds the first $k$ coordinates. Let $s, t: \R^k \to \R^{d-k}$ be any functions (in practice, neural networks), and let $\odot$ be the elementwise product:

1. Jacobian. $y_A$ does not depend on $x_B$, so the upper-right block is $0$. Among the coordinates of $x_B$, each $y_{B,i}$ depends only on $x_{B,i}$, so the lower-right block is diagonal:

2. Log-determinant. The matrix is block lower-triangular, so $\log|\det Dg| = \sum_i s_i(x_A)$. The block $\partial y_B/\partial x_A$ holds the derivatives of the networks $s$ and $t$, but it never enters the determinant.
3. Inverse. $x_A = y_A$ and $x_B = (y_B - t(y_A)) \odot \exp(-s(y_A))$. This needs only forward passes of $s$ and $t$, never their inverses. So $s$ and $t$ can be any networks.
4. Price. One layer does not change $x_A$, and it acts on $x_B$ as an affine map. So we stack layers and swap (or permute) the roles of $A$ and $B$ between layers.

Interactive figure: loading…

Autoregressive flows (MAF, IAF: $y_i$ depends on $x_1, \dots, x_i$) also have triangular Jacobians. Continuous flows $dz/dt = f(z, t)$ use the continuous version $\frac{d}{dt}\log p(z(t)) = -\operatorname{tr}\big(\partial f/\partial z\big)$; the trace is estimated without forming the Jacobian (lesson 08). Both are treated in the generative models track (normalizing flows, flow matching).

Answers to section 1. (1) $p_Y(1) = 1/2$: $g^{-1}(y) = y/2$, so $p_Y(y) = 1 \cdot \frac12$ on $(0, 2)$. (2) No. The mode of $Y = e^X$ is at $1/e \approx 0.368$ (computed in section 4). A nonlinear map does not carry the mode along, because the Jacobian factor $1/y$ is not constant. The mean does not move that way either: $\mathbb E[e^X] = e^{1/2} = 1.6487$, not $e^0$. (3) No. Its area is about $r\,\Delta r\,\Delta\theta$, because the arc side has length $r\,\Delta\theta$. Only at $r = 1$ is $\Delta r\,\Delta\theta$ right.

4. Worked example #

(a) A Gaussian bump over a disk, in polar coordinates. Compute $I(R) = \iint_{x^2 + y^2 \le R^2} e^{-(x^2 + y^2)}\,dx\,dy$.

1. In polar coordinates the disk is $0 \le r \le R$, $0 \le \theta < 2\pi$, the integrand is $e^{-r^2}$, and $dx\,dy = r\,dr\,d\theta$.
2. The integrand is a product of a function of $\theta$ (the constant 1) and a function of $r$, so by Fubini: $I(R) = \int_0^{2\pi} d\theta \cdot \int_0^R e^{-r^2}\,r\,dr = 2\pi \int_0^R e^{-r^2}\,r\,dr$.
3. Substitute $u = r^2$, $du = 2r\,dr$: $\int_0^R e^{-r^2}\,r\,dr = \frac12\int_0^{R^2} e^{-u}\,du = \frac12\big(1 - e^{-R^2}\big)$.
4. So $I(R) = \pi\big(1 - e^{-R^2}\big)$. For $R = 1$: $\pi(1 - e^{-1}) = 1.9859$. As $R \to \infty$, $I \to \pi$.

In Cartesian coordinates, $e^{-x^2}$ has no elementary antiderivative. The extra factor $r$ from polar coordinates is exactly what makes the substitution possible. Lesson 11 uses $I(\infty) = \pi$ to get $\int e^{-x^2}dx = \sqrt\pi$.

(b) The log-normal density. $X \sim N(0, 1)$ and $Y = e^X$. Then $g^{-1}(y) = \log y$ and $(g^{-1})'(y) = 1/y > 0$, so for $y > 0$:

• $p_Y(1) = \varphi(0) = 0.3989$ and $p_Y(e) = \varphi(1)/e = 0.2420/2.7183 = 0.0890$.
• Mode. $\frac{d}{dy}\log p_Y(y) = -\frac{\log y}{y} - \frac1y = 0$ gives $\log y = -1$, so $y = 1/e \approx 0.368$. There $p_Y(1/e) = \varphi(-1)\,e = 0.2420 \cdot 2.7183 = 0.6577$. The term $-\log y$ (the log-Jacobian) moves the peak away from $e^0 = 1$.

(c) One coupling layer in 2-D. The base is $X \sim N(0, I_2)$. Keep $x_1$, and transform $x_2$ with $s(x_1) = x_1/2$ and $t(x_1) = x_1^2$: $y_1 = x_1$, $y_2 = x_2\,e^{x_1/2} + x_1^2$. Evaluate $\log p_Y$ at $y = (1, 2)$. (This is the "coupling layer" map of the second figure.)

1. Invert. $x_1 = y_1 = 1$, and $x_2 = (y_2 - t(1))\,e^{-s(1)} = (2 - 1)\,e^{-0.5} = 0.6065$.
2. Jacobian at $x$. $\partial y_2/\partial x_1 = x_2 e^{s}\,s'(x_1) + t'(x_1) = 1 \cdot 0.5 + 2 = 2.5$ (because $x_2 e^{s} = y_2 - t = 1$). So $Dg = \begin{pmatrix}1 & 0\\ 2.5 & e^{0.5}\end{pmatrix}$, $\det Dg = e^{0.5} = 1.6487$, and $\log|\det Dg| = 0.5 = s(x_1)$. The entry $2.5$ plays no role.
3. Base log-density. $\log p_X(x) = -\log 2\pi - \tfrac12(1^2 + 0.6065^2) = -1.8379 - 0.6839 = -2.5218$.
4. Result. $\log p_Y(1, 2) = -2.5218 - 0.5 = -3.0218$, so $p_Y(1, 2) = 0.0487$.

In closed form, $p_Y(y_1, y_2) = \varphi(y_1)\,\varphi\big((y_2 - y_1^2)e^{-y_1/2}\big)\,e^{-y_1/2}$. It must integrate to 1; the code checks this.

5. In code #

NumPy and SciPy check the hand computations. The notebook has these cells and the exercises.

import numpy as np
from scipy import integrate, stats

# (1) midpoint rule for the integral of x^2 on [0, 1], and Monte Carlo for E[X^2]
for n in [4, 8]:
    m = (np.arange(n) + 0.5) / n                 # midpoints, h = 1/n
    M = np.mean(m**2)                            # h * sum f(m_i)
    print(n, M, 1/3 - M)
# 4 0.328125 0.0052083...
# 8 0.33203125 0.0013020...      (error / 4 when n doubles)
rng = np.random.default_rng(0)
fx = rng.standard_normal(10_000) ** 2
print(fx.mean().round(4), (fx.std(ddof=1) / np.sqrt(len(fx))).round(4))
# 0.9962 0.014                   (true value 1, standard error about 0.0141)

# (2) polar coordinates: integral of exp(-(x^2 + y^2)) over the unit disk
val, _ = integrate.dblquad(lambda y, x: np.exp(-(x*x + y*y)), -1, 1,
                           lambda x: -np.sqrt(1 - x*x), lambda x: np.sqrt(1 - x*x))
print(round(val, 5), round(np.pi * (1 - np.exp(-1)), 5))   # 1.98587 1.98587

# (3) log-normal density phi(log y) / y against scipy
phi = stats.norm.pdf
for y in [1, np.e, np.exp(-1)]:
    print(round(phi(np.log(y)) / y, 4), round(stats.lognorm(s=1).pdf(y), 4))
# 0.3989 0.3989 / 0.089 0.089 / 0.6577 0.6577

# (4) the coupling layer of example (c): Jacobian by central differences, then log p
def g(x):
    return np.array([x[0], x[1] * np.exp(x[0] / 2) + x[0]**2])

x = np.array([1.0, np.exp(-0.5)])
print(g(x))                                     # [1. 2.]
eps = 1e-6
J = np.column_stack([(g(x + eps * e) - g(x - eps * e)) / (2 * eps) for e in np.eye(2)])
print(J.round(4))                               # [[1. 0.] [2.5 1.6487]]: lower triangular
logp = stats.norm.logpdf(x).sum() - np.log(abs(np.linalg.det(J)))
print(round(logp, 4))                           # -3.0218

# the closed-form density of Y integrates to 1
p_y = lambda y2, y1: phi(y1) * phi((y2 - y1**2) * np.exp(-y1 / 2)) * np.exp(-y1 / 2)
total, _ = integrate.dblquad(p_y, -9, 9, lambda a: a**2 - 12 * np.exp(a / 2),
                             lambda a: a**2 + 12 * np.exp(a / 2))
print(round(total, 6))                          # 1.0

In real code, never compute det and then log: for large $d$ the determinant overflows or underflows. Add up the logs of the diagonal of a triangular Jacobian, or use np.linalg.slogdet. Libraries such as torch.distributions.TransformedDistribution do exactly the sum of section 3.7: they invert each transform in reverse order and subtract each log_abs_det_jacobian.

6. Where it shows up #

• Training losses are Monte Carlo estimates. The average loss over a mini-batch estimates the expected loss over the data distribution, with an error that falls like $1/\sqrt{B}$ in the batch size $B$. The ELBO of a VAE and the policy-gradient objective are expectations estimated by sampling; their gradients are the topic of lesson 12.
• Normalizing flows. RealNVP and Glow (coupling layers), MAF and IAF (autoregressive layers) are trained with the exact log-likelihood $\log p_Z(z_0) - \sum_k \log|\det J_k|$. See the generative models track.
• Continuous flows and diffusion likelihoods. Neural ODE flows, flow matching and the probability-flow ODE of diffusion models compute $\log p$ by integrating $-\operatorname{tr}(\partial f/\partial z)$ along the trajectory. The PDE behind it is conservation of probability, the continuity equation $\partial_t p + \nabla \cdot (p f) = 0$.
• Sampling. Inverse-CDF sampling ($X = F^{-1}(U)$ for $U$ uniform); Box–Muller, which uses polar coordinates to turn two uniform numbers into two independent $N(0, 1)$ samples ($X = \sqrt{-2\log U_1}\cos(2\pi U_2)$); and the reparameterization $x = \mu + \sigma \odot \epsilon$ of VAEs.
• Priors and parameterizations. A flat prior on $\sigma$ is not flat on $\log\sigma$: $p(\log\sigma) = p(\sigma)\,\sigma$. So MAP estimates and the Laplace approximation (lesson 11) change when you reparameterize. The maximum likelihood estimate does not, because the likelihood is a function of the parameter, not a density over it.
• ML for science. Models often predict log-transformed targets (log-amplitudes of ground motion, log-precipitation, log-permeability). A model trained on $u = \log y$ has $-\log p_Y(y) = -\log p_U(\log y) + \log y$. So you must add the mean of $\log y$ before you compare its negative log-likelihood with a model trained on $y$. In the same way, a flow on images that rescales pixels from $[0, 256)$ to $[0, 1)$ must add the log-determinant of that scaling before it reports bits per dimension.

7. Common mistakes #

• Wrong direction. $p_Y(y) = p_X(x)/|\det Dg(x)|$, not $p_X(x)\,|\det Dg(x)|$. Check with $Y = 2X$: the density must go down.
• Wrong sign in the log form. Going forward through a layer, you subtract $\log|\det Dg|$. Going backward (from data to base), you add $\log|\det Dg^{-1}|$. They are the same number with opposite signs.
• Forgetting the absolute value. A reflection has $\det = -1$, but densities and volumes are never negative.
• Forgetting the $r$ in polar coordinates. $dx\,dy = r\,dr\,d\theta$, not $dr\,d\theta$.
• Missing branches. For $Y = X^2$ you need both $\pm\sqrt y$; one branch gives half the density.
• Moving the mode or the mean with the map. The mode of $g(X)$ is not $g(\text{mode of } X)$, and $\mathbb E[g(X)] \ne g(\mathbb E[X])$ in general ($\mathbb E[e^X] = e^{1/2}$ for $X \sim N(0, 1)$).
• "A density is at most 1." A density is probability per unit volume. It can be any large number, as long as it integrates to 1.
• Swapping the order of integration without checking. Fubini needs $\int|f| < \infty$ (or $f \ge 0$ for Tonelli); otherwise the two orders can give different answers.
• A full determinant for a stacked flow. Never build the $d \times d$ Jacobian of the whole flow. Add the per-layer log-determinants, and use $\log$ of the diagonal, not $\log(\det)$.
• Dimension change. The density formula needs $g: \R^d \to \R^d$. A map from $\R^d$ into $\R^D$ with $D > d$ has no density on $\R^D$; the volume factor on its image surface is $\sqrt{\det(Dg^\top Dg)}$.

8. Check yourself #

1. $X$ is uniform on $(0, 1)$ and $Y = -\log X$. What is $p_Y$?
   Answer

   $g$ is decreasing; $g^{-1}(y) = e^{-y}$ and $|(g^{-1})'(y)| = e^{-y}$. So $p_Y(y) = 1 \cdot e^{-y}$ for $y > 0$: $Y \sim \mathrm{Exp}(1)$. This is inverse-CDF sampling for the exponential distribution.
2. $x \sim N(0, I_2)$ and $y = Ax$ with $A = \begin{pmatrix}1 & 1\\0 & 2\end{pmatrix}$. Compute $\log p_Y(3, 2)$.
   Answer

   $x = A^{-1}y = (2, 1)$; $\log p_X(x) = -\log 2\pi - 5/2 = -4.3379$; $\log|\det A| = \log 2 = 0.6931$; so $\log p_Y = -5.0310$. Check: $Y \sim N(0, AA^\top)$ with $AA^\top = \begin{pmatrix}2 & 2\\2 & 4\end{pmatrix}$, $\det = 4$, $y^\top(AA^\top)^{-1}y = 5$, and $-\log 2\pi - \frac12\log 4 - \frac52 = -5.0310$.
3. $X \sim N(0, 1)$ and $Y = |X|$. What is $p_Y$?
   Answer

   Two branches, $x = \pm y$, each with $|dx/dy| = 1$: $p_Y(y) = 2\varphi(y)$ for $y > 0$ (the half-normal density). For example $p_Y(1) = 0.4839$.
4. Use polar coordinates to compute $\iint_{x^2 + y^2 \le 4} (x^2 + y^2)\,dx\,dy$.
   Answer

   $\int_0^{2\pi}\!\int_0^2 r^2 \cdot r\,dr\,d\theta = 2\pi \cdot \frac{2^4}{4} = 8\pi \approx 25.13$. Without the factor $r$ you would get $2\pi \cdot \frac83 = 16.76$, which is wrong.
5. A 3-layer flow maps $z_0 \mapsto z_1 \mapsto z_2 \mapsto y$. At one sample, $\log p_Z(z_0) = -3.0$, and the forward log-determinants are $0.2$, $-0.5$ and $1.1$. What is $\log p_Y(y)$?
   Answer

   $-3.0 - (0.2 - 0.5 + 1.1) = -3.8$.
6. Model A is trained on $u = \log y$ and has a test negative log-likelihood (NLL) of $1.20$ nats per point on the $u$ scale. Model B is trained on $y$ and has a test NLL of $1.50$ nats on the $y$ scale. The test mean of $\log y$ is $0.5$. Which model is better?
   Answer

   Convert A to the $y$ scale: $-\log p_Y(y) = -\log p_U(\log y) + \log y$, so its NLL is $1.20 + 0.5 = 1.70 > 1.50$. Model B is better.

Recall cardsAnswer in your head first, then reveal. Add the cards to your review deck to see them again just before you would forget.

1. State the fundamental theorem of calculus (both parts).

   Show answer

2. Write the 1-D substitution rule with an absolute value.

   Show answer

3. Write $\mathbb E[f(X)]$ for a random variable $X$ with density $p$.

   Show answer

4. What is the standard error of a Monte Carlo average of $N$ independent samples?

   Show answer

5. When may you swap the order of a double integral?

   Show answer

6. Why is $dx\,dy = r\,dr\,d\theta$ in polar coordinates?

   Show answer

7. State the change-of-variables formula for integrals.

   Show answer

8. Give the log-density of $Y = g(X)$ for a diffeomorphism $g$.

   Show answer

9. $X \sim N(0, 1)$. What is the density of $Y = X^2$, and why?

   Show answer

10. What is $\log|\det J|$ of the affine coupling layer $y_B = x_B \odot e^{s(x_A)} + t(x_A)$, and why?

    Show answer

Add these cards to my reviewOpen the review →

9. Interview angle #

• "Derive the density of $Y = e^X$ for $X \sim N(\mu, \sigma^2)$." — $p_Y(y) = N(\log y; \mu, \sigma^2)/y$ for $y > 0$; the mode is $e^{\mu - \sigma^2}$ and the mean is $e^{\mu + \sigma^2/2}$.
• "Why must a normalizing flow be invertible, and why do people use coupling layers?" — The likelihood needs $g^{-1}$ and $\log|\det|$. Coupling layers make both cheap ($O(d)$ log-determinant, closed-form inverse), while $s$ and $t$ can still be any networks.
• "Write the log-likelihood of a flow." — $\log p_Z(g^{-1}(y)) + \sum_k \log|\det Dg_k^{-1}|$ along the inverse path (equivalently, minus the forward log-determinants).
• "Is the MAP estimate invariant under reparameterization? The MLE?" — The MLE is; the MAP is not, because the posterior density picks up the Jacobian of the parameter map.
• "Why is Monte Carlo used for high-dimensional integrals?" — Its error $\sigma/\sqrt N$ does not depend on the dimension, while a grid rule needs exponentially many points in $d$.
• "How does a continuous normalizing flow compute $\log p$?" — It integrates $-\operatorname{tr}(\partial f/\partial z)$ along the ODE, with a stochastic (Hutchinson) trace estimate.

10. Practice #

The notebook for this lesson runs in Colab. It has five exercises, each with a check cell and a folded solution:

1. Write riemann_mid and trapezoid, measure their error orders, and write a Monte Carlo estimator with a standard error.
2. Write pushforward_logpdf_1d for maps with one or more branches, and check it against SciPy and against a histogram.
3. Write affine_logpdf for $Y = AX + b$ with slogdet, and compare it with scipy.stats.multivariate_normal.
4. Write coupling_forward and coupling_inverse with their log-determinants, and check that the density of a 2-layer flow integrates to about 1 on a grid.
5. Compute $\iint e^{-(x^2 + y^2)}$ over a disk on a Cartesian grid and with a 1-D integral in polar coordinates.

11. Watch and read #

Integration and the fundamental theorem of calculus | Chapter 8, Essence of calculus3Blue1Brown · Essence of calculus, chapter 8 · 21 min

The integral as a limit of sums of thin rectangles, and why the area function has derivative f.

Why “probability of 0” does not mean “impossible” | Probabilities of probabilities, part 23Blue1Brown · 10 min

What a probability density is, and why it can be larger than 1.

The Jacobian DeterminantKhan Academy (Grant Sanderson) · Multivariable calculus · 9 min

The determinant of the Jacobian as the local area scale factor of a nonlinear map.

Lec 18: Change of variables | MIT 18.02 Multivariable Calculus, Fall 2007Denis Auroux · MIT OpenCourseWare 18.02 · 50 min

Polar coordinates and general substitutions in double integrals, with the Jacobian.

What are Normalizing Flows?Ari Seff · 13 min

The change-of-variables formula for densities, composition of layers, and coupling layers.

• Deisenroth, Faisal and Ong, Mathematics for Machine Learning, free PDF: §6.4 (expectations) and §6.7 (change of variables, theorem 6.16).
• Papamakarios, Nalisnick, Rezende, Mohamed and Lakshminarayanan, Normalizing Flows for Probabilistic Modeling and Inference (JMLR 2021): §2.1 (change of variables and compositions), §3.1 (autoregressive flows and coupling layers), §4 (continuous-time flows).
• Dinh, Sohl-Dickstein and Bengio, Density estimation using Real NVP (ICLR 2017): the affine coupling layer.
• UvA Deep Learning tutorial 11: normalizing flows, sections "Normalizing Flows as generative model" and "Coupling layers".
• Art Owen, Monte Carlo theory, methods and examples, free chapters: chapter 2 (simple Monte Carlo and its error).
• MIT OpenCourseWare, 18.02 Multivariable Calculus (Fall 2007): lecture notes on double integrals, polar coordinates and change of variables.