Least Squares from Gaussian Likelihood

Derivation

Least Squares from Gaussian Likelihood

D-120 Home PU-106 Threads chance Depends on The Method of Maximum Likelihood, The Central Limit Theorem

Statement

Given N independent measurements y_i at known inputs x_i, each scattered about a model prediction f(x_i; θ) by a Gaussian error of standard deviation σ_i, the parameters θ that maximise the joint likelihood are exactly those that minimise the weighted sum of squared residuals χ² = Σ_i (y_i − f(x_i; θ))² / σ_i². Maximum likelihood and least squares coincide.

Why it matters

Least squares is usually taught as an ad hoc recipe: "minimise the squared residuals". This derivation shows it is not arbitrary. It is the precise consequence of one physical assumption — that measurement noise is Gaussian and independent — and that assumption is itself justified by the central limit theorem, since real noise is a sum of many small independent perturbations.

The result also tells you exactly when least squares is the right thing to do, and when it is not. The moment the errors are non-Gaussian, correlated, or heteroscedastic in a way you ignore, the equivalence breaks and blind least squares becomes statistically suboptimal or biased. Knowing where the equivalence comes from is knowing where it ends.

Assumptions

Errors are Gaussian.If dropped, the log-likelihood is no longer a sum of squared residuals; the optimal estimator changes shape entirely (e.g. Laplacian errors give a sum of absolute residuals, the L¹ / median estimator). Errors are statistically independent across data points.If dropped, the joint likelihood no longer factorises into a product, so the log-likelihood picks up cross terms and χ² must use the full inverse covariance matrix rather than a plain diagonal sum. The standard deviations σ_i are known and fixed (not parameters to be fitted).If dropped, the normalisation prefactor 1/(σ_i√2π) is no longer constant in the fit, so the −Σ ln σ_i term must be retained and minimising χ² alone is wrong. Errors act only on the responses y_i; the inputs x_i are exact.If dropped, vertical residuals are the wrong metric; the correct problem is errors-in-variables / total least squares, and ordinary least squares becomes biased (attenuation). The model f(x; θ) is the true generating function for some θ.If dropped (model misspecification), the estimator still returns a "best fit", but likelihood maximisation no longer targets a true parameter and interval interpretation fails.

Derivation

p(y_i | θ) = (1 / σ_i√2π) · exp[ −(y_i − f(x_i; θ))² / 2σ_i² ]

The probability density of one measurement given the model is a Gaussian centred on the prediction with width σ_i. This is the Gaussian-error assumption written out. A

L(θ) = ∏_i=1^N p(y_i | θ) = ∏_i=1^N (1 / σ_i√2π) · exp[ −(y_i − f(x_i; θ))² / 2σ_i² ]

Independence lets the joint likelihood factorise into a product of the single-point densities. This is the only place independence is used. A

ln L(θ) = −½ Σ_i (y_i − f(x_i; θ))² / σ_i² − Σ_i ln(σ_i√2π)

Take the logarithm. It is strictly increasing, so it has the same argmax as L; the product becomes a sum and the exponentials linearise. A

ln L(θ) = −½ χ²(θ) + C, χ²(θ) ≡ Σ_i (y_i − f(x_i; θ))² / σ_i², C = −Σ_i ln(σ_i√2π)

Separate the θ-dependent part from the constant. Because the σ_i are fixed, C does not depend on θ and cannot affect the location of the maximum. B

argmax_θ ln L(θ) = argmax_θ [ −½ χ²(θ) ] = argmin_θ χ²(θ)

Adding the constant C and multiplying by the negative constant −½ leaves the location of the extremum unchanged but flips maximisation into minimisation. This is the equivalence. B

∂χ²/∂θ_k = −2 Σ_i (y_i − f(x_i; θ)) / σ_i² · ∂f(x_i; θ)/∂θ_k = 0

The interior minimum satisfies stationarity in every parameter. For a linear model f = Σ_k θ_k g_k(x) this yields the linear normal equations with a unique solution; for nonlinear f it is solved iteratively (Gauss–Newton). C

Result

θ̂_MLE = argmax_θ L(θ) = argmin_θ χ²(θ), χ²(θ) = Σ_i=1^N (y_i − f(x_i; θ))² / σ_i²

Reading. Fitting by weighted least squares is not a convenient approximation to the "real" statistical answer — under independent Gaussian errors it is the maximum-likelihood answer, exactly. Each residual is weighted by 1/σ_i², so precise points (small σ) pull the fit hard and noisy points count for little.

Units check. Each term (y_i − f)²/σ_i² is [y]² / [y]² = dimensionless, so χ² is a pure number, as any quantity destined to be compared with N degrees of freedom must be. The likelihood L has units of 1/[y]^N (a product of densities), but only its argmax is used, so its dimension is irrelevant to the result.

Limiting cases

Equal errors, σ_i = σ: the weights drop out and χ² ∝ Σ(y_i − f)² — plain ordinary least squares.
Constant model, f = μ: the estimator is the inverse-variance weighted mean μ̂ = Σ(y_i/σ_i²) / Σ(1/σ_i²); with equal σ it collapses to the arithmetic mean.
Linear model, f = a + bx: stationarity gives closed-form normal equations for a and b — the familiar straight-line least-squares formulae.
Single data point, N = 1: χ² = 0 is achievable by any model passing through it; parameters are undetermined, mirroring an unbounded likelihood ridge.

Breaks when

Errors are heavy-tailed or contain outliers. For non-Gaussian noise ln L is no longer −½χ²; least squares over-weights outliers (a squared penalty), so the MLE is a different, robust estimator (Huber, L¹, or a Student-t likelihood). Least squares still runs, but it is no longer optimal and is easily wrecked by a single bad point.
Errors are correlated. With a non-diagonal covariance matrix C the joint density does not factorise; the correct statistic is χ² = r^TC⁻¹r with off-diagonal terms. Summing plain squared residuals then mis-weights the data and gives wrong parameters and wrong error bars.
The inputs x_i also carry error. Vertical-residual least squares is biased (regression attenuation); the likelihood must marginalise over the true x, giving total least squares / errors-in-variables regression.
σ_i is itself unknown and estimated jointly. The −Σlnσ_i term is no longer constant, so minimising χ² alone is incomplete and can drive σ → 0 spuriously; the full log-likelihood must be optimised.

Failure modes

Dropping the 1/σ_i² weights when the errors genuinely differ — running ordinary least squares on heteroscedastic data, letting noisy points bias the fit.
Believing least squares requires Gaussian errors to be defined. It does not; it is just that its optimality as the MLE requires them. Gauss–Markov guarantees it is the best linear unbiased estimator regardless.
Discarding the −Σlnσ_i term while treating σ as a fitted parameter — valid only when σ is fixed and known.
Confusing the minimised value χ²_min with the goodness-of-fit statistic. The equivalence is about the argmin; interpreting χ²_min/(N−p) needs the errors to be correctly scaled.
Minimising Σ|r_i| (absolute deviations) and calling it least squares — that is the MLE for Laplacian, not Gaussian, errors.

Discussion

Historically this equivalence runs the other way from how we teach it. Gauss defined the normal distribution partly by demanding that the arithmetic mean be the most probable estimate of a location parameter; the Gaussian is essentially the unique error law for which maximum likelihood reproduces least squares. So "Gaussian errors" and "least squares is optimal" are two faces of the same statement, not an accident that they happen to line up.

The central limit theorem supplies the physical reason to expect Gaussian errors in the first place: whenever the total measurement error is a sum of many small, independent contributions — thermal jitter, readout noise, mechanical vibration, quantisation — its distribution tends to Gaussian regardless of the individual sources. This is why least squares is the default across the physical sciences, and why it fails precisely in the regimes where one contribution dominates or the tails are fat (a single cosmic-ray hit, a stuck bit) — those are exactly the non-Gaussian, outlier-prone cases where the CLT justification does not apply.

The result also seeds the whole framework of parameter uncertainties. Near the minimum, χ²(θ) ≈ χ²_min + ½(θ − θ̂)^TH(θ − θ̂) with H the Hessian of χ²; because ln L = −½χ² + const, this is a Gaussian likelihood in θ with covariance 2H⁻¹, so the parameter covariance is the inverse of the curvature of χ². The familiar "Δχ² = 1 for a 1σ interval" rule is just this quadratic expansion. Adding a Gaussian prior on θ turns the same calculation into ridge (Tikhonov) regression, and with a flat prior the maximum-a-posteriori estimate coincides with the MLE — a clean bridge to the Bayesian view.

Common misconceptions. Least squares does not assume the data are Gaussian — it assumes the errors are; y_i can span any range as long as the residuals scatter normally. And the equivalence says nothing about whether the model is correct: minimising χ² always returns a best fit, but only a valid model plus correct σ_i make χ²_min/(N−p) ≈ 1 meaningful.

Worked examples

Example 1 — Inverse-variance weighted mean. Three independent measurements of one voltage, with different precisions: y = 5.10 ± 0.10 V, 5.30 ± 0.20 V, 5.00 ± 0.05 V. Combine them by maximum likelihood (constant model f = μ).

μ̂ = Σ_i(y_i/σ_i²) / Σ_i(1/σ_i²), σ_μ = 1/√(Σ_i 1/σ_i²)

Setting ∂χ²/∂μ = 0 for f = μ gives the weighted mean; curvature of χ² gives its uncertainty. B

w_i = 1/σ_i² = 100, 25, 400 V⁻²; Σw_i = 525 V⁻²

Compute the weights; the σ = 0.05 V point dominates. A

Σw_iy_i = 100(5.10) + 25(5.30) + 400(5.00) = 2642.5 V⁻¹

Numerator of the weighted mean. A

μ̂ = 2642.5 / 525 = 5.033 V; σ_μ = 1/√525 = 0.0436 V

Divide, then take the inverse square root of the summed weights. A

μ̂ = 5.033 ± 0.044 V

Reading. The combined precision (±0.044 V) beats even the best single measurement (±0.05 V), and the answer sits close to that best point, as inverse-variance weighting demands.

Example 2 — Straight-line fit with equal errors. Data (x, y): (1, 2.1), (2, 3.9), (3, 6.1), (4, 8.0), each with σ = 0.2. Fit f = a + bx by least squares (= MLE here) and report χ².

b̂ = (NΣx_iy_i − Σx_iΣy_i) / (NΣx_i² − (Σx_i)²), â = (Σy_i − b̂Σx_i)/N

Normal equations from ∂χ²/∂a = ∂χ²/∂b = 0 with equal σ (weights cancel). B

N = 4, Σx = 10, Σy = 20.1, Σx² = 30, Σxy = 60.2

Accumulate the sums from the four points. A

b̂ = (4·60.2 − 10·20.1)/(4·30 − 10²) = 39.8/20 = 1.99

Insert numbers into the slope formula. A

â = (20.1 − 1.99·10)/4 = 0.2/4 = 0.05

Back-substitute the slope into the intercept formula. A

r_i = 0.06, −0.13, 0.08, −0.01; χ² = Σr_i²/σ² = 0.027/0.04 = 0.68

Residuals against f = 0.05 + 1.99x, weighted by 1/σ². B

y = 0.05 + 1.99 x, χ² = 0.68 (ν = N−2 = 2)

Reading. χ²/ν = 0.34 < 1: the fit is excellent — arguably the quoted σ = 0.2 slightly over-states the scatter. The near-zero intercept and unit-ish slope match a clean linear relation.

Problems

(A) Combine two measurements. A length is measured as 10.2 ± 0.3 cm and 10.8 ± 0.4 cm. Give the maximum-likelihood combined value and its uncertainty.
Solution
Weights w₁ = 1/0.3² = 11.11, w₂ = 1/0.4² = 6.25 cm⁻²; Σw = 17.36. Σwy = 11.11(10.2) + 6.25(10.8) = 113.33 + 67.50 = 180.83. μ̂ = 180.83/17.36 = 10.42 cm. σ_μ = 1/√17.36 = 0.24 cm. Result: 10.42 ± 0.24 cm.
(A) Mean from χ². For a constant model f = μ with equal errors σ, show that minimising χ² gives the arithmetic mean.
Solution
χ² = Σ_i(y_i − μ)²/σ². Then dχ²/dμ = −(2/σ²)Σ_i(y_i − μ) = 0 ⇒ Σ_iy_i = Nμ ⇒ μ̂ = (1/N)Σ_iy_i = ȳ. The second derivative 2N/σ² > 0 confirms a minimum. The MLE of a common mean under equal Gaussian errors is the sample mean.
(B) Derive the weighted-mean uncertainty. With μ̂ = Σw_iy_i/Σw_i and w_i = 1/σ_i², show σ_μ² = 1/Σw_i by error propagation.
Solution
μ̂ is linear in the independent y_i, so σ_μ² = Σ_i(∂μ̂/∂y_i)²σ_i² with ∂μ̂/∂y_i = w_i/Σw. Hence σ_μ² = Σ_i w_i²σ_i²/(Σw)². Since w_iσ_i² = 1, the numerator is Σw_i, giving σ_μ² = Σw_i/(Σw_i)² = 1/Σw_i = 1/Σ(1/σ_i²). Equivalently, 1/σ_μ² is the curvature ∂²(½χ²)/∂μ².
(C) σ as a parameter — the N vs N−1 bias. For N points with a common unknown σ and known mean, write the full log-likelihood, find the MLE σ̂², and show it is biased.
Solution
−ln L = Σ_i[ r_i²/(2σ²) + lnσ ] + const, with r_i the residuals. Then ∂(−ln L)/∂σ = −(1/σ³)Σr_i² + N/σ = 0 ⇒ σ̂² = (1/N)Σ_ir_i². When the residuals are taken about a fitted mean, one degree of freedom is used up, so E[Σr_i²] = (N−1)σ² and E[σ̂²] = (N−1)σ²/N < σ². The MLE divides by N and is biased low; the unbiased estimator divides by N−1. This is why the −Σlnσ term cannot be dropped once σ is fitted.
(C) Correlated pair. Two measurements y₁ = 1.0, y₂ = 1.4 of a common value μ, with σ₁ = σ₂ = 0.2 and correlation ρ = 0.5. Using χ² = r^TC⁻¹r, find μ̂ and σ_μ, and comment on the effect of ρ.
Solution
C = 0.04·[[1, 0.5],[0.5, 1]], so C⁻¹ = (1/0.04)(1/0.75)[[1, −0.5],[−0.5, 1]] = 33.33·[[1, −0.5],[−0.5, 1]]. With design vector 1 = (1,1)^T: 1^TC⁻¹1 = 33.33(1 − 0.5 − 0.5 + 1) = 33.33; 1^TC⁻¹y = 33.33·0.5(y₁+y₂) = 33.33(1.2) = 40.0. Then μ̂ = (1^TC⁻¹y)/(1^TC⁻¹1) = 40.0/33.33 = 1.20, and σ_μ² = 1/(1^TC⁻¹1) = 1/33.33 = 0.030, so μ̂ = 1.20 ± 0.17. The estimate is still the plain average (equal σ), but the variance 0.030 exceeds the uncorrelated value σ²/2 = 0.020: positive correlation means the two points partly share their error, so averaging removes less of it. Ignoring ρ would understate the uncertainty.