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

The section forks here — choose a route

With standard errors in hand, the curvature route continues into hypothesis testing. Or cross to the sampling route, which prices uncertainty by walking the surface instead of assuming its shape. Neither route requires the other first.