25. The delta method and Slutsky's theorem

The CLT gives you the limiting distribution of \bar X_n. But you usually want something derived from it — a ratio, a log, an odds ratio, a standard deviation. And you usually have to estimate \sigma rather than know it.

Two tools close both gaps.

Slutsky's theorem

Suppose X_n \xrightarrow{d} X and Y_n \xrightarrow{p} c for a constant c. Then:

X_n + Y_n \xrightarrow{d} X + c, \qquad X_nY_n \xrightarrow{d} cX, \qquad \frac{X_n}{Y_n} \xrightarrow{d} \frac{X}{c} \;\;(c \ne 0)

The requirement that c be a constant is essential — the theorem is false if Y_n converges to a genuine random variable.

Why it matters: replacing \sigma with s

The CLT says

\frac{\bar X_n - \mu}{\sigma/\sqrt n} \xrightarrow{d} N(0,1)

But \sigma is unknown. In practice we use the sample standard deviation s. Is that legitimate?

Write it as a product:

\frac{\bar X_n - \mu}{s/\sqrt n} = \underbrace{\frac{\bar X_n - \mu}{\sigma/\sqrt n}}_{\xrightarrow{d}\; N(0,1)} \times \underbrace{\frac{\sigma}{s}}_{\xrightarrow{p}\; 1}

Since s \xrightarrow{p} \sigma (by the LLN applied to squared deviations), the second factor converges in probability to the constant 1. Slutsky then gives:

\frac{\bar X_n - \mu}{s/\sqrt n} \xrightarrow{d} N(0,1)

That is the justification for every large-sample z-test and confidence interval you will ever compute. Estimating \sigma costs you nothing asymptotically.

For small n the cost is real, and that's precisely the t distribution's correction (§1, and again in §5). Slutsky says the correction vanishes as n \to \infty — which is why t_k \to N(0,1).

The delta method

You know \bar X_n is asymptotically Normal. What about g(\bar X_n)?

Statement. If

\sqrt n\,(T_n - \theta) \xrightarrow{d} N(0, \sigma^2)

and g is differentiable at \theta with g'(\theta) \ne 0, then

\sqrt n\,\big(g(T_n) - g(\theta)\big) \xrightarrow{d} N\!\left(0,\; \sigma^2\big[g'(\theta)\big]^2\right)

The idea is a first-order Taylor expansion. Near \theta:

g(T_n) \approx g(\theta) + g'(\theta)(T_n - \theta)

So g(T_n) - g(\theta) \approx g'(\theta)(T_n - \theta) — a linear function of something already known to be asymptotically Normal. Scaling a Normal by g'(\theta) multiplies its variance by [g'(\theta)]^2.

The remainder is negligible because T_n \to \theta, so the expansion is evaluated over a shrinking neighbourhood where the linear approximation gets arbitrarily good.

In practical form:

g(\bar X_n) \;\approx\; N\!\left(g(\mu),\; \frac{\sigma^2[g'(\mu)]^2}{n}\right)

Read the consequence: a transformation that is steep at \mu amplifies uncertainty; one that is flat damps it. The derivative is the amplification factor.

When g'(\theta) = 0

The method breaks — the linear term vanishes and the whole approximation is the remainder. Then you expand to second order, and the limit is a chi-squared rather than a Normal:

n\big(g(T_n) - g(\theta)\big) \xrightarrow{d} \frac{\sigma^2 g''(\theta)}{2}\chi^2_1

Note the scaling changes from \sqrt n to n: convergence is faster when the function is flat, which makes sense.

Worked example

X_1, \dots, X_n are Bernoulli(p). Find the asymptotic distribution of the log-odds \log\frac{\hat p}{1 - \hat p}.

By the CLT, with \sigma^2 = p(1-p):

\sqrt n\,(\hat p - p) \xrightarrow{d} N\big(0,\; p(1-p)\big)

Take g(x) = \log\frac{x}{1-x} = \log x - \log(1-x). Differentiate:

g'(x) = \frac{1}{x} + \frac{1}{1-x} = \frac{1}{x(1-x)}

Apply the delta method:

\sqrt n\left(g(\hat p) - g(p)\right) \xrightarrow{d} N\!\left(0,\; p(1-p)\cdot\frac{1}{p^2(1-p)^2}\right) = N\!\left(0, \frac{1}{p(1-p)}\right)

So the log-odds has asymptotic standard error

\text{SE} = \frac{1}{\sqrt{n\,p(1-p)}}

This is exactly the standard error reported for coefficients in logistic regression (§6) — the delta method is where it comes from.

Notice something useful. For \hat p itself, the variance p(1-p)/n depends heavily on p. Near p = 0.5 it's largest; near 0 or 1 it collapses, and the Normal approximation gets poor because the distribution is squashed against a boundary.

On the log-odds scale there is no boundary — the transformation maps (0,1) to all of \mathbb{R} — and the Normal approximation is far better in the tails. This is why odds ratios are analysed on the log scale and only converted back at the end. Choosing a transformation to make the Normal approximation behave is called a variance-stabilising transformation, and it's a standard move.

Doing it in Python

Slutsky in action — swapping \sigma for s costs nothing at large n:

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

rng = np.random.default_rng(0)
trials, mu, sigma = 40_000, 5.0, 2.0

print(f"{'n':>6} {'KS p (known sigma)':>20} {'KS p (estimated s)':>20}")
for n in (5, 10, 30, 100, 1000):
    x = rng.normal(mu, sigma, size=(trials, n))
    known = (x.mean(axis=1) - mu) / (sigma / np.sqrt(n))
    est = (x.mean(axis=1) - mu) / (x.std(axis=1, ddof=1) / np.sqrt(n))
    print(f"{n:>6} {kstest(known, 'norm').pvalue:>20.4f} {kstest(est, 'norm').pvalue:>20.4f}")

print("\nWith sigma known, standard Normal at every n.")
print("With s estimated, it's t-distributed -- but converges to Normal as n grows.")

The delta method for the log-odds, checked by simulation:

import numpy as np

rng = np.random.default_rng(1)
n, p, trials = 500, 0.3, 100_000

phat = rng.binomial(n, p, size=trials) / n
phat = np.clip(phat, 1e-9, 1 - 1e-9)          # guard the log at the boundary
logodds = np.log(phat / (1 - phat))

theory_se = 1 / np.sqrt(n * p * (1 - p))
print(f"log-odds simulated sd: {logodds.std():.3f}")
print(f"delta-method SE      : {theory_se:.3f}")
print()
print("centre simulated:", round(logodds.mean(), 5))
print("g(p) = log(p/(1-p)):", round(np.log(p / (1 - p)), 5))

Where the delta method breaks — g'(\theta) = 0:

import numpy as np
from scipy.stats import skew, kurtosis

rng = np.random.default_rng(2)
n, trials = 2_000, 100_000

# X ~ N(0, 1); estimate mu = 0 and transform with g(x) = x^2, so g'(0) = 0.
xbar = rng.standard_normal((trials, n)).mean(axis=1)
g = xbar**2

print("g'(0) = 0, so the delta method's Normal limit does not apply.")
print("  skewness of n*g:", round(float(skew(n * g)), 3), " (a Normal would be 0)")
print("  excess kurtosis:", round(float(kurtosis(n * g)), 3), " (a Normal would be 0)")
print("\nchi-squared(1) has skewness 2.83 and excess kurtosis 12 -- that's the limit here.")

Your turn

1. \bar X_n estimates \mu = 4 with \sigma = 2, n = 100. Find the approximate distribution of \sqrt{\bar X_n}.

2. Why does replacing \sigma by s not change the asymptotic distribution?

3. \hat p = 0.5 from n = 100. Find the SE of \hat p and of the log-odds.

Solutions

1. g(x) = \sqrt x, so g'(x) = \frac{1}{2\sqrt x} and g'(4) = \frac{1}{4}.

Centre: g(\mu) = \sqrt 4 = 2.

Variance: $\frac{\sigma^2[g'(\mu)]^2}{n} = \frac{4 \times (1/4)^2}{100} = \frac{0.25}{100} = 0.0025$.

\sqrt{\bar X_n} \;\approx\; N(2,\; 0.0025), \qquad \text{SE} = 0.05

Compare: \bar X_n itself has SE = 2/\sqrt{100} = 0.2. The square root has damped the uncertainty by the factor g'(4) = 0.25, because it's a flattening transformation at that point.

2. Because s \xrightarrow{p} \sigma, so \sigma/s \xrightarrow{p} 1 — a constant. Slutsky's theorem then says multiplying a sequence converging in distribution by a sequence converging in probability to a constant just scales the limit, here by 1:

\frac{\bar X_n - \mu}{s/\sqrt n} = \frac{\bar X_n - \mu}{\sigma/\sqrt n} \cdot \frac{\sigma}{s} \xrightarrow{d} N(0,1) \times 1

The intuition: with enough data, s is so close to \sigma that the difference is swamped by the sampling variability of \bar X_n itself.

For small n it is not negligible, the ratio is genuinely random, and the result is a t distribution with heavier tails — which is exactly why the t exists.

3. For the proportion:

\text{SE}(\hat p) = \sqrt{\frac{\hat p(1-\hat p)}{n}} = \sqrt{\frac{0.25}{100}} = 0.05

For the log-odds, using the delta-method result:

\text{SE} = \frac{1}{\sqrt{n\hat p(1-\hat p)}} = \frac{1}{\sqrt{100 \times 0.25}} = \frac{1}{5} = 0.2

The log-odds SE is four times larger — which is just $g'(0.5) = \frac{1}{0.5 \times 0.5} = 4$, the amplification factor. Different scales carry different uncertainties, and the delta method converts between them.

Check yourself in code

Verify the delta method for the log-odds transformation: simulate \hat p, transform it, and compare the empirical standard deviation against 1/\sqrt{np(1-p)}.

Print exactly this:

simulated sd 0.098
delta-method SE 0.098
centre matches: True

Use default_rng(1), n = 500, p = 0.3, and 100000 trials. Print both standard deviations to 3 decimal places. Report centre matches: True if the simulated mean log-odds is within 0.01 of \log\frac{p}{1-p}.

Why 3 decimals and not more? A standard deviation estimated from 100,000 draws carries a relative error of roughly 1/\sqrt{2 \times 10^5} \approx 0.2\%, which lands on the third decimal of 0.098. Printing a fourth digit would be quoting noise as if it were signal — and, in practice, that digit also shifts between NumPy versions, because the algorithm behind binomial changes with n. Report only the digits your simulation actually supports.

import numpy as np

rng = np.random.default_rng(1)
n, p, trials = 500, 0.3, 100_000

phat = rng.binomial(n, p, size=trials) / n
phat = np.clip(phat, 1e-9, 1 - 1e-9)
logodds = np.log(phat / (1 - phat))

print(f"simulated sd {logodds.std():.3f}")

# Print the delta-method SE 1/sqrt(n*p*(1-p)), then check the simulated
# centre against log(p/(1-p)) to within 0.01.

import numpy as np

rng = np.random.default_rng(1)
n, p, trials = 500, 0.3, 100_000

phat = rng.binomial(n, p, size=trials) / n
phat = np.clip(phat, 1e-9, 1 - 1e-9)
logodds = np.log(phat / (1 - phat))

print(f"simulated sd {logodds.std():.3f}")

se = 1 / np.sqrt(n * p * (1 - p))
print(f"delta-method SE {se:.3f}")

target = np.log(p / (1 - p))
print("centre matches:", bool(abs(logodds.mean() - target) < 0.01))

---

Slutsky's theorem lets you swap an unknown constant for a consistent estimate of it without changing the limiting distribution — which is why estimating \sigma is free at large n. The delta method pushes asymptotic normality through a smooth function, with the derivative acting as the amplification factor on the standard error.

That closes §3. Next: with the limit theory in place, we can finally build estimators and say how good they are.