Transformations and the Jacobian method · Probability & Statistics

You know the distribution of X. You need the distribution of Y = g(X).

This comes up constantly — you measure a radius but need the area, you model log-returns but need returns, you generate uniforms but need normals. There are two reliable methods, and one very common mistake.

The mistake first

You cannot simply substitute. If f_X is the density of X and Y = g(X), it is not true that f_Y(y) = f_X(g^{-1}(y)).

The reason is that a density is probability per unit length, and g stretches or compresses the axis. If g doubles the spacing between points, the same probability now spreads over twice the width, so the density must halve.

Getting that stretch factor right is the entire content of this lesson.

Method 1: the CDF method (always works)

Go through the CDF, where no such correction is needed — because a CDF is a probability, not a density.

$F_Y(y) = P(Y \le y) = P(g(X) \le y)$

Solve the inequality for X, evaluate using F_X, then differentiate to get the density.

This is the method to reach for when in doubt. It handles non-monotone functions, discrete pieces, and boundary cases that trip up the formula below.

Worked through: Y = X^2 with X \sim N(0,1)

$F_Y(y) = P(X^2 \le y) = P(-\sqrt y \le X \le \sqrt y) = \Phi(\sqrt y) - \Phi(-\sqrt y)$

for y \ge 0. By the symmetry of the standard Normal, \Phi(-a) = 1 - \Phi(a):

$F_Y(y) = 2\Phi(\sqrt y) - 1$

Differentiate, using the chain rule (\frac{d}{dy}\sqrt y = \frac{1}{2\sqrt y}):

$f_Y(y) = 2\varphi(\sqrt y)\cdot\frac{1}{2\sqrt y} = \frac{\varphi(\sqrt y)}{\sqrt y} = \frac{1}{\sqrt{2\pi y}}e^{-y/2}$

That is the chi-squared with 1 degree of freedom — precisely as §1 claimed when it defined \chi^2_k as a sum of k squared standard Normals.

Note that g(x) = x^2 is not monotone on \mathbb{R}: both +\sqrt y and -\sqrt y map to y. The CDF method handled that automatically, which is exactly why it's the safe default.

Method 2: the change-of-variables formula (monotone g)

When g is monotone and differentiable, the CDF method's result can be packaged into a formula:

$\boxed{\;f_Y(y) = f_X\big(g^{-1}(y)\big)\left|\frac{d}{dy}g^{-1}(y)\right|\;}$

The derivative term is the stretch correction. The absolute value handles decreasing g, where the derivative is negative but a density can't be.

Writing x = g^{-1}(y), it's often clearer as

$f_Y(y) = f_X(x)\left|\frac{dx}{dy}\right|$

which reads as: probability mass in dx equals probability mass in dy, so f_Y(y)|dy| = f_X(x)|dx|.

The multivariate case: the Jacobian

Now (X_1, X_2) \to (Y_1, Y_2). The stretch factor becomes the Jacobian determinant — the local volume-scaling factor of the transformation:

$J = \det\begin{pmatrix} \dfrac{\partial x_1}{\partial y_1} & \dfrac{\partial x_1}{\partial y_2} \\[2ex] \dfrac{\partial x_2}{\partial y_1} & \dfrac{\partial x_2}{\partial y_2} \end{pmatrix}$

$f_{Y_1,Y_2}(y_1, y_2) = f_{X_1,X_2}(x_1, x_2)\,\big|J\big|$

The recipe:

Solve for the old variables in terms of the new: x_1 = \dots, x_2 = \dots
Build the matrix of partial derivatives \partial x_i / \partial y_j
Take its determinant, take the absolute value
Substitute into f_X and multiply

Watch the direction. The Jacobian must be of the old variables with respect to the new ones. Getting it upside down inverts your answer — and since J_{\text{wrong}} = 1/J_{\text{right}}, the error is easy to make and easy to miss.

Finally, work out the new support: where the transformed variables can actually live. This is the step most often forgotten, and it's usually where the interesting constraints appear.

Worked example

X \sim \text{Uniform}(0,1). Find the density of Y = -\ln X.

g(x) = -\ln x is decreasing on (0,1), so the formula applies.

Invert: y = -\ln x \implies x = e^{-y}.

Range: as x runs over (0,1), y runs over (0, \infty).

Derivative: \frac{dx}{dy} = -e^{-y}, so \left|\frac{dx}{dy}\right| = e^{-y}.

Assemble, remembering f_X = 1 on (0,1):

$f_Y(y) = 1 \cdot e^{-y} = e^{-y}, \qquad y > 0$

That's the Exponential with rate 1.

This is not a curiosity — it's the inverse transform method, and it's how exponential random numbers are actually generated: draw a uniform, take -\ln U. We'll build that out properly in §11.

Sanity check the density: \int_0^\infty e^{-y}dy = 1. ✓ And a rough intuition check: X near 0 maps to large y, and X near 0 is where the transformation stretches hardest, which is why the density decays there.

Doing it in Python

Verify a derived density by simulating the transformation and comparing against the formula:

import numpy as np

rng = np.random.default_rng(0)
U = rng.uniform(0, 1, 500_000)
Y = -np.log(U)                      # claim: Exponential(1)

print("simulated mean:", round(Y.mean(), 4), " theory 1.0")
print("simulated var :", round(Y.var(), 4), "  theory 1.0")

# Compare the empirical density with e^{-y} on a few bins
counts, edges = np.histogram(Y, bins=np.linspace(0, 5, 11), density=True)
centres = (edges[:-1] + edges[1:]) / 2
print(f"\n{'y':>6} {'empirical':>11} {'e^-y':>9}")
for c, d in zip(centres, counts):
    print(f"{c:>6.2f} {d:>11.4f} {np.exp(-c):>9.4f}")

The X^2 result, checked against SciPy's chi-squared:

import numpy as np
from scipy.stats import chi2, norm

rng = np.random.default_rng(1)
Z = rng.standard_normal(500_000)
Y = Z**2

print("simulated mean:", round(Y.mean(), 4), " chi2(1) mean:", chi2.mean(1))
print("simulated var :", round(Y.var(), 4), " chi2(1) var :", chi2.var(1))

for q in (0.25, 0.5, 0.9, 0.99):
    print(f"quantile {q}: simulated {np.quantile(Y, q):8.4f}   chi2(1) {chi2.ppf(q, 1):8.4f}")

And a direct demonstration that naive substitution is wrong — the stretch factor is not optional:

import numpy as np
from scipy.integrate import quad

# X ~ Uniform(0,1), Y = X**2. Correct density: f_Y(y) = 1/(2*sqrt(y)).
correct = lambda y: 1 / (2 * np.sqrt(y))
naive = lambda y: 1.0                       # "just substitute" -> still uniform

print("integral of correct density:", round(quad(correct, 1e-12, 1)[0], 6))
print("integral of naive density  :", round(quad(naive, 0, 1)[0], 6))

rng = np.random.default_rng(2)
Y = rng.uniform(0, 1, 200_000) ** 2
print("\nsimulated P(Y < 0.25) =", round((Y < 0.25).mean(), 4))
print("correct   P(Y < 0.25) =", round(quad(correct, 1e-12, 0.25)[0], 4))
print("naive     P(Y < 0.25) =", round(quad(naive, 0, 0.25)[0], 4), " <- wrong")

Both densities integrate to 1, so that check alone won't catch the error — but the naive version gets the actual probability badly wrong.

Your turn

1. X \sim \text{Uniform}(0,1). Find the density of Y = X^2.

2. X \sim N(\mu, \sigma^2). Show Y = \frac{X - \mu}{\sigma} is N(0,1).

3. X \sim \text{Exponential}(\lambda). Find the density of Y = X^{1/2}.

Solutions

1. g(x) = x^2 is increasing on (0,1), so the formula applies.

Invert: x = \sqrt y, with y \in (0,1). Derivative: \frac{dx}{dy} = \frac{1}{2\sqrt y}.

$f_Y(y) = 1 \cdot \frac{1}{2\sqrt y} = \frac{1}{2\sqrt y}, \qquad 0 < y < 1$

Check: \int_0^1 \frac{1}{2\sqrt y}dy = [\sqrt y]_0^1 = 1. ✓

Note the density is unbounded as y \to 0 — it blows up. That's legal (§1: only the area is capped), and it reflects squaring compressing values near 0 into a tiny interval, concentrating the probability there.

2. Linear, increasing since \sigma > 0.

Invert: x = \mu + \sigma y. Derivative: \frac{dx}{dy} = \sigma.

$f_Y(y) = f_X(\mu + \sigma y)\cdot\sigma = \frac{1}{\sigma\sqrt{2\pi}}\exp\!\left(-\frac{(\mu + \sigma y - \mu)^2}{2\sigma^2}\right)\cdot \sigma$

The \sigma from the Jacobian cancels the \sigma in the normalising constant, and the exponent simplifies to -\frac{\sigma^2y^2}{2\sigma^2} = -\frac{y^2}{2}:

$f_Y(y) = \frac{1}{\sqrt{2\pi}}e^{-y^2/2}$

Standard Normal. ✓ This is the standardisation from §1, now actually proved — and note it only works because the Jacobian factor was included.

3. g(x) = \sqrt x is increasing on (0, \infty).

Invert: x = y^2, with y > 0. Derivative: \frac{dx}{dy} = 2y.

$f_Y(y) = \lambda e^{-\lambda y^2}\cdot 2y = 2\lambda y\,e^{-\lambda y^2}, \qquad y > 0$

This is the Rayleigh distribution. It's what you get for the distance from the origin to a point whose two coordinates are independent Normals — the distribution of wind speeds, and of the magnitude of 2-D random noise.

Check yourself in code

Confirm that Y = -\ln U turns a uniform into an Exponential(1), by comparing the simulated mean, variance, and a tail probability against theory.

Print exactly this:

mean 1.0
var 1.0
P(Y > 2) 0.1358
theory   0.1353

Round the mean and variance to 1 decimal place, and both probabilities to 4. Use numpy.random.default_rng(0) and 500000 draws — the simulated tail probability lands near, but not exactly on, the theoretical e^{-2}; that gap is sampling noise, and §3 will tell you how big to expect it to be.

import numpy as np

rng = np.random.default_rng(0)
U = rng.uniform(0, 1, 500_000)
Y = -np.log(U)

print("mean", round(Y.mean(), 1))

# Print the variance to 1 dp, then the simulated P(Y > 2) and the
# theoretical value e^{-2}, both to 4 dp.

import numpy as np

rng = np.random.default_rng(0)
U = rng.uniform(0, 1, 500_000)
Y = -np.log(U)

print("mean", round(Y.mean(), 1))
print("var", round(Y.var(), 1))
print("P(Y > 2)", round((Y > 2).mean(), 4))
print("theory  ", round(float(np.exp(-2)), 4))

A density can't just be substituted through a function — it has to be scaled by how much the transformation stretches the axis. In one dimension that's |dx/dy|; in several it's the absolute Jacobian determinant, taken of the old variables with respect to the new. And when the function isn't monotone, fall back to the CDF method, which never needs the correction at all.

Next: what happens when you sort a sample and ask about its smallest, largest, or middle value.

18. Transformations and the Jacobian method