27. The method of moments

The last lesson judged estimators. This one constructs them — with the oldest and simplest recipe in statistics, due to Karl Pearson (1894).

The idea

A distribution's theoretical moments are functions of its parameters:

\mu_1' = E[X] = g_1(\theta_1, \dots, \theta_k), \qquad \mu_2' = E[X^2] = g_2(\theta_1, \dots, \theta_k), \;\dots

The corresponding sample moments are computable from data:

m_1' = \bar X = \frac{1}{n}\sum X_i, \qquad m_2' = \frac{1}{n}\sum X_i^2, \;\dots

By the Law of Large Numbers, m_k' \xrightarrow{p} \mu_k' — sample moments converge to theoretical ones.

So set them equal and solve for the parameters.

\mu_k'(\theta) = m_k' \quad \text{for } k = 1, \dots, \text{(as many as there are parameters)}

One equation per unknown parameter, starting from the lowest moment.

The recipe

1. Count the parameters — say k of them.
2. Write the first k theoretical moments as functions of the parameters.
3. Set each equal to the corresponding sample moment.
4. Solve the system.

That's it. No calculus, no optimisation.

Worked example 1: the Normal

Estimate \mu and \sigma^2 from a Normal sample.

Two parameters, so two equations.

E[X] = \mu \;\Longrightarrow\; \hat\mu = \bar X

E[X^2] = \sigma^2 + \mu^2 \;\Longrightarrow\; \hat\sigma^2 + \hat\mu^2 = \frac{1}{n}\sum X_i^2

Substitute \hat\mu = \bar X:

\hat\sigma^2 = \frac{1}{n}\sum X_i^2 - \bar X^2 = \frac{1}{n}\sum(X_i - \bar X)^2

Note this is the divide-by-n version — the biased one from last lesson. The method of moments does not care about bias; it just matches moments.

Worked example 2: the Uniform

X_1, \dots, X_n \sim \text{Uniform}(0, \theta). Estimate \theta.

One parameter, so one equation:

E[X] = \frac{\theta}{2} = \bar X \;\Longrightarrow\; \hat\theta = 2\bar X

Which is exactly the estimator we compared last lesson — and found to be more than twice as bad as the sample maximum on MSE.

This is the method's central weakness. It has no idea that X_{(n)} exists, because the mean discards the information carried by the largest observation.

Worse, the estimate can be logically impossible. Suppose you observe 1, 2, 12. Then \bar X = 5 and \hat\theta = 10 — but you saw a 12, so \theta \ge 12. The method of moments has produced an estimate ruled out by the data itself.

Worked example 3: the Gamma

Estimate the shape k and rate \lambda of a Gamma.

E[X] = \frac{k}{\lambda}, \qquad \operatorname{Var}(X) = \frac{k}{\lambda^2}

Using the variance is equivalent to using the second moment and often tidier. Setting them equal to the sample versions:

\frac{\hat k}{\hat\lambda} = \bar X, \qquad \frac{\hat k}{\hat\lambda^2} = \hat\sigma^2

Divide the second into the first:

\hat\lambda = \frac{\bar X}{\hat\sigma^2}, \qquad \hat k = \bar X \hat\lambda = \frac{\bar X^2}{\hat\sigma^2}

Clean closed forms. The maximum likelihood estimates for a Gamma, by contrast, require solving an equation involving the digamma function numerically — which is exactly when the method of moments earns its keep.

Properties

Consistent. Sample moments converge to theoretical moments (LLN), so if the solving step is a continuous function of them, the estimates converge too.

Asymptotically Normal, usually — by the CLT plus the delta method (§3).

Often biased, as the Normal variance case showed.

Not efficient. It generally has larger variance than maximum likelihood, because it uses only the first few moments and discards the rest of the data's information.

Can fail outright. Estimates may fall outside the parameter space (a negative variance, a \hat\theta below an observed value), and for distributions without finite moments — Cauchy again — the method has nothing to work with at all.

Its virtue is simplicity. Closed forms, no optimisation, no convergence failures. In practice it's often used to generate starting values for the iterative algorithms that maximise a likelihood.

Doing it in Python

The Gamma example, where the closed form is genuinely useful:

import numpy as np
from scipy.stats import gamma

rng = np.random.default_rng(0)
true_k, true_lambda = 3.0, 2.0
x = rng.gamma(shape=true_k, scale=1/true_lambda, size=5_000)

xbar = x.mean()
var = x.var()                       # divide by n, matching the method

lam_hat = xbar / var
k_hat = xbar**2 / var

print(f"true    k = {true_k},  lambda = {true_lambda}")
print(f"MoM     k = {k_hat:.4f}, lambda = {lam_hat:.4f}")

# Compare with maximum likelihood (numerical, from SciPy)
k_mle, loc, scale_mle = gamma.fit(x, floc=0)
print(f"MLE     k = {k_mle:.4f}, lambda = {1/scale_mle:.4f}")

The uniform failure, made concrete:

import numpy as np

# A small sample where the method of moments contradicts the data
sample = np.array([1.0, 2.0, 12.0])
mom = 2 * sample.mean()

print("sample      :", sample)
print("MoM estimate:", mom)
print("largest observed value:", sample.max())
print("impossible? ", bool(mom < sample.max()),
      "-- theta cannot be below a value we actually saw")

How often does that happen? More than you'd guess:

import numpy as np

rng = np.random.default_rng(1)
theta = 10.0
print(f"{'n':>6} {'P(MoM estimate < observed max)':>32}")
for n in (3, 5, 10, 30, 100):
    s = rng.uniform(0, theta, size=(200_000, n))
    impossible = (2 * s.mean(axis=1) < s.max(axis=1)).mean()
    print(f"{n:>6} {impossible:>32.4f}")

And a head-to-head against MLE on the Normal variance:

import numpy as np

rng = np.random.default_rng(2)
mu, sigma, n, trials = 0.0, 2.0, 10, 200_000
x = rng.normal(mu, sigma, size=(trials, n))

mom = x.var(axis=1, ddof=0)          # method of moments: divide by n
unb = x.var(axis=1, ddof=1)          # unbiased: divide by n-1

for name, est in [("MoM (/n)", mom), ("unbiased (/n-1)", unb)]:
    bias = est.mean() - sigma**2
    print(f"{name:>18}  bias {bias:>8.4f}  MSE {est.var() + bias**2:>8.4f}")
print("\nMoM is biased low -- but note its MSE is actually the smaller of the two.")

Your turn

1. X \sim \text{Exponential}(\lambda). Find the method of moments estimator of \lambda.

2. X \sim \text{Bernoulli}(p). Find \hat p.

3. X \sim \text{Uniform}(a, b) with both endpoints unknown. Find \hat a and \hat b.

Solutions

1. One parameter, one equation. Since E[X] = 1/\lambda:

\frac{1}{\hat\lambda} = \bar X \implies \hat\lambda = \frac{1}{\bar X}

Sensible: a rate is the reciprocal of a mean waiting time. (This one happens to coincide with the MLE, though it is biased — E[1/\bar X] \ne 1/E[\bar X], by Jensen again.)

2. E[X] = p, so

\hat p = \bar X = \frac{\text{number of successes}}{n}

The obvious estimator, and here it's also the MLE and unbiased.

3. Two parameters, two equations:

E[X] = \frac{a+b}{2} = \bar X, \qquad \operatorname{Var}(X) = \frac{(b-a)^2}{12} = \hat\sigma^2

From the second, b - a = \sqrt{12\hat\sigma^2} = 2\sqrt{3}\,\hat\sigma. With a + b = 2\bar X, solve the pair:

\hat a = \bar X - \sqrt3\,\hat\sigma, \qquad \hat b = \bar X + \sqrt3\,\hat\sigma

The same objection applies: these are built only from the mean and variance, so nothing stops \hat a from exceeding the smallest observation. The sensible estimators here are X_{(1)} and X_{(n)}, which the method of moments cannot find.

Check yourself in code

Fit a Gamma by the method of moments and confirm the closed-form estimates land close to the truth.

Print exactly this:

k hat 3.0182
lambda hat 1.9861
k within 5%: True
lambda within 5%: True

Use default_rng(0), true k = 3, \lambda = 2, and 5000 samples. Round both estimates to 4 decimal places. Compute the variance with ddof=0, matching the method.

import numpy as np

rng = np.random.default_rng(0)
true_k, true_lambda = 3.0, 2.0
x = rng.gamma(shape=true_k, scale=1 / true_lambda, size=5_000)

xbar, var = x.mean(), x.var()

lam_hat = xbar / var
k_hat = xbar**2 / var

print("k hat", round(k_hat, 4))
# Print lambda hat, then whether each estimate is within 5% of the truth.

import numpy as np

rng = np.random.default_rng(0)
true_k, true_lambda = 3.0, 2.0
x = rng.gamma(shape=true_k, scale=1 / true_lambda, size=5_000)

xbar, var = x.mean(), x.var()

lam_hat = xbar / var
k_hat = xbar**2 / var

print("k hat", round(k_hat, 4))
print("lambda hat", round(lam_hat, 4))
print("k within 5%:", bool(abs(k_hat - true_k) / true_k < 0.05))
print("lambda within 5%:", bool(abs(lam_hat - true_lambda) / true_lambda < 0.05))

---

Match theoretical moments to sample moments and solve. It's consistent, usually asymptotically Normal, and often the only method with a closed form — but it is frequently biased, statistically inefficient, and can produce estimates the data plainly contradicts.

Next: the method that fixes those problems by asking a better question — which parameter value makes the data we actually saw most probable?