Standard Errors from Curvature

From one curvature to a matrix of them

On the previous page, one parameter meant one curvature: a single number describing how quickly the log-likelihood falls away from the peak. Now consider the model from the 2D optimisation page — a straight line through data, $y = \beta_0 + \beta_1 x$ — fitted to our synthetic resting heart rate sample (heart rate in bpm against weekly exercise hours, $n = 60$).

The log-likelihood is now a surface over the $(\beta_0, \beta_1)$ plane, and a surface can fall away at different rates in different directions. It can also fall away diagonally in a coordinated way, which turns out to encode how the two estimates co-vary. To capture all of this we need the curvature in every direction — the Hessian matrix of second derivatives:

$$H = \begin{pmatrix} \dfrac{\partial^2 \ell}{\partial \beta_0^2} & \dfrac{\partial^2 \ell}{\partial \beta_0 \partial \beta_1} \\[1em] \dfrac{\partial^2 \ell}{\partial \beta_1 \partial \beta_0} & \dfrac{\partial^2 \ell}{\partial \beta_1^2} \end{pmatrix}\Bigg|_{\beta = \hat\beta}$$

Interactive: probing the summit, barefoot and blind

JonStats introduces the Hessian with a memorable image: you are standing at the summit, barefoot and blindfolded, and your only instrument is your feet — you can step a little way in a chosen direction and feel how far the ground drops. Three probes are enough to map the local shape completely. Try them:

The pipeline: Hessian → information → variance → standard errors

Once the Hessian is measured at the peak, standard errors follow in three mechanical steps. This is the multi-parameter version of $\text{SE} = 1/\sqrt{-\ell''}$, and it is the computation running behind every summary() table you have seen in the tutorials (numbers below: our heart rate model at $n=60$):

Interactive: the confidence ellipse

The variance–covariance matrix has a direct picture: an ellipse around the summit containing the parameter values compatible with the data. Its width in each direction comes from the diagonal of $\Sigma$; its tilt comes from the covariance. Change the sample size and watch the ellipse shrink — standard errors fall as $1/\sqrt{n}$. Then scatter simulated draws from $\text{MVN}(\hat\beta, \Sigma)$ — the sampling distribution that the curvature approximation implies.

The optimiser already knew

Here is the elegant part. On the 4D+ page, Newton–Raphson used the Hessian to decide how far to step: sharp curvature → the peak is close → small steps; gentle curvature → a long plateau → stride out. The same matrix, evaluated at the summit it helped find, prices the uncertainty of the estimates. Fitting and inference are two uses of one object: the curvature that guides the climb also calibrates the confidence.

And the purple ellipse stamped without explanation on the Multi-Optima page's representativeness comparison? It was this confidence ellipse all along.

In code: extract the Hessian, invert, done

Both snippets fit the resting heart rate line by hand-rolled maximum likelihood, then read the standard errors straight off the inverted (negated) Hessian. Outputs are from real runs on docs/data/inference-rest-hr.json ($n=60$) — validation scripts: scripts/R/validate-inference.R, scripts/py/generate_inference_data.py.

llNormal <- function(pars, y, X) {
  beta   <- pars[1:ncol(X)]
  sigma2 <- exp(pars[ncol(X) + 1])        # keep sigma^2 positive
  -1/2 * sum(log(2 * pi * sigma2) + (y - (X %*% beta))^2 / sigma2)
}

fit <- optim(par = c(0, 0, 0), fn = llNormal, method = "BFGS",
             control = list(fnscale = -1),   # maximise, not minimise
             hessian = TRUE,                 # <-- ask for the curvature
             y = y, X = cbind(1, x))

vcov_ml <- solve(-fit$hessian)   # negate, then invert
sqrt(diag(vcov_ml))              # standard errors

# Simulate the implied sampling distribution (JonStats-style):
draws <- MASS::mvrnorm(n = 10000, mu = fit$par[1:2],
                       Sigma = vcov_ml[1:2, 1:2])

# point estimates:      71.7714  -1.2710   (log sigma^2: 3.6995)
# solve(-fit$hessian):
#           [,1]    [,2]
# [1,]    2.9675 -0.4415
# [2,]   -0.4415  0.0850
# standard errors:       1.7226   0.2915
#
# cross-check with lm(y ~ x): coef 71.7714, -1.2711; SEs 1.7521, 0.2965
# (lm's SEs are slightly larger: lm divides the residual variance by
#  n - 2 where maximum likelihood divides by n)

import numpy as np
import statsmodels.api as sm
from scipy import optimize

X = sm.add_constant(x)

# The convenience route: statsmodels computes vcov from the curvature
ols = sm.OLS(y, X).fit()
print(ols.params, ols.bse)         # estimates and standard errors
print(ols.cov_params())            # the full variance-covariance matrix

# The by-hand route: minimise the negative log-likelihood ...
def negll(pars):
    b0, b1, eta = pars
    s2 = np.exp(eta)
    mu = b0 + b1 * x
    return 0.5 * np.sum(np.log(2 * np.pi * s2) + (y - mu)**2 / s2)

fit = optimize.minimize(negll, x0=np.zeros(3), method="BFGS")

# ... then invert the Hessian of the negative log-likelihood
# (already the information matrix, so no negation needed here)
H = numeric_hessian(negll, fit.x)      # central finite differences
vcov = np.linalg.inv(H)
print(np.sqrt(np.diag(vcov)))

# OLS  coefficients:  b0 = 71.7714, b1 = -1.2711
# OLS  std errors:    se(b0) = 1.7521, se(b1) = 0.2965
# ML   std errors from inverse Hessian:
#                     se(b0) = 1.7226, se(b1) = 0.2915
# vcov (ML):
#   [[ 2.9675, -0.4415],
#    [-0.4415,  0.0850]]

The small print