⚙️ Setup¶
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
rng = np.random.default_rng(43)
plt.rcParams.update({"figure.dpi":110,"font.size":11,"axes.spines.top":False,"axes.spines.right":False})
VIOLET="#7c3aed"; PINK="#db2777"; TEAL="#0d9488"
print("ready")
ready
X = rng.normal(50, 10, size=200_000)
a, b = 2, 5
Y = a*X + b
print(f"X: mean {X.mean():.2f}, sd {X.std():.2f}")
print(f"Y = 2X+5: mean {Y.mean():.2f} (theory {a*50+b}), sd {Y.std():.2f} (theory {abs(a)*10})")
X: mean 49.96, sd 9.97 Y = 2X+5: mean 104.92 (theory 105), sd 19.94 (theory 20)
A linear transform relocates and rescales but does not change the family: a linear function of a normal is still normal. This is exactly the standardization of Chapter 39 run in reverse, and it is why z = (X − mu)/sigma turns any normal into the standard normal.
X = rng.normal(0, 0.5, size=200_000)
Y = np.exp(X) # lognormal
print(f"X is symmetric normal: mean {X.mean():.3f}, skew {stats.skew(X):.3f}")
print(f"Y = exp(X) is right-skewed: mean {Y.mean():.3f}, skew {stats.skew(Y):.3f}")
X is symmetric normal: mean 0.000, skew 0.000 Y = exp(X) is right-skewed: mean 1.134, skew 1.772
fig,(a1,a2)=plt.subplots(1,2,figsize=(8.6,3))
a1.hist(X,bins=60,density=True,color=VIOLET,alpha=0.8); a1.set_title("X ~ Normal (symmetric)")
a2.hist(Y,bins=80,density=True,color=PINK,alpha=0.8); a2.set_title("Y = exp(X): lognormal"); a2.set_xlim(0,6)
plt.tight_layout(); plt.show()
The symmetric normal becomes a skewed lognormal. Nonlinear transforms reshape distributions, and the change-of-variables formula makes this precise: the new density is the old density divided by the magnitude of the transform's derivative, the Jacobian factor.
# verify for Y = exp(X), X ~ N(0,1). The lognormal pdf is the result.
xs = np.linspace(0.05, 6, 300)
analytic = stats.lognorm.pdf(xs, s=1) # scipy lognormal, sigma=1
# build from the formula: f_Y(y) = f_X(ln y) * |d/dy ln y| = phi(ln y)/y
from_formula = stats.norm.pdf(np.log(xs)) / xs
print(f"max |formula - scipy lognormal| = {np.abs(from_formula-analytic).max():.2e}")
max |formula - scipy lognormal| = 1.11e-16
The hand-built density phi(ln y)/y matches SciPy's lognormal exactly. The 1/y is the Jacobian: it is |d/dy ln(y)|, correcting for the fact that the log compresses large values. Every density transform, however exotic, is this same recipe, old density times the inverse-transform's Jacobian.
U = rng.random(200_000)
# generate exponential(rate=1): inverse CDF is -ln(1-U)
X = -np.log(1 - U)
print(f"inverse-transform sample: mean {X.mean():.3f} (exponential mean 1.0)")
print(f"compare numpy exponential: mean {rng.exponential(1.0, 200_000).mean():.3f}")
inverse-transform sample: mean 0.997 (exponential mean 1.0) compare numpy exponential: mean 1.000
fig,ax=plt.subplots(figsize=(7,3))
ax.hist(X,bins=80,density=True,color=VIOLET,alpha=0.75,label="inverse-transform samples")
xs=np.linspace(0,8,200); ax.plot(xs, np.exp(-xs), color=PINK, lw=2.5, label="exponential pdf")
ax.set_title("Uniforms -> exponential via the inverse CDF"); ax.set_xlim(0,8); ax.legend()
plt.tight_layout(); plt.show()
Plain uniform random numbers, pushed through the inverse CDF F⁻¹, become exponential samples that match the target density. This inverse-transform trick is how computers turn a uniform random generator into samples from any distribution you like.
# reparameterization: z = mu + sigma * epsilon, with epsilon ~ N(0,1)
mu, sigma = 3.0, 2.0
eps = rng.normal(0, 1, size=200_000)
z = mu + sigma*eps # a sample from N(mu, sigma) as a transform of fixed noise
print(f"reparameterized sample: mean {z.mean():.3f} (target {mu}), sd {z.std():.3f} (target {sigma})")
print("Because z is a DIFFERENTIABLE function of (mu, sigma), gradients can flow through the sampling step.")
reparameterized sample: mean 2.999 (target 3.0), sd 2.004 (target 2.0) Because z is a DIFFERENTIABLE function of (mu, sigma), gradients can flow through the sampling step.
The reparameterization trick writes a random sample z = mu + sigma·epsilon as a transform of fixed noise epsilon, so the randomness sits in epsilon and the parameters (mu, sigma) enter through a differentiable function, letting backpropagation train variational autoencoders. Normalizing flows stack many invertible transforms (tracking each Jacobian, exactly Demo 3) to morph simple noise into complex data distributions. Change of variables is the mathematics of generative AI.