Beyond the Quadratic: Bayesian Uncertainty

The hidden assumption

The curvature method does something subtle: it measures the surface only at the peak, then extrapolates. From three numbers — two curvatures and a cross-term — it reconstructs an entire uncertainty picture, as if the summit's immediate steepness continued forever in a perfect quadratic bowl. Equivalently: it assumes the uncertainty is shaped like a single multivariate normal hill, symmetric about the MLE.

For large samples this is usually an excellent bargain — three numbers instead of mapping a whole surface. But the assumption can fail, and it fails most often exactly where GLMs live: parameters with boundaries (probabilities near 0 or 1, rates near 0), small samples, and skewed surfaces where the likelihood falls gently on one side of the peak and steeply on the other.

Interactive: when the parabola lies

A pilot safety study observes 2 adverse events among 20 patients, so $\hat{p} = 0.1$. The teal curve is the actual (normalised) likelihood of $p$ — with a flat prior, also the posterior density. The purple dashed curve is what the curvature method assumes instead: a normal distribution centred at $\hat p$. Move the sliders and watch when the two stories agree and when they don't.

Walking instead of assuming

What is the alternative when the surface isn't quadratic? Measure it — all of it that matters. You have already seen the machinery: on the MCMC page, chains wandered across Arthur's Seat, preferring higher ground but never settling, and the fraction of time they spent in each place converged to the posterior distribution. There the visit-map was a picture; here it becomes the uncertainty estimate itself. Percentiles of the visited values give credible intervals, their spread gives posterior standard deviations — and none of it assumes the surface is a single symmetric MVN-shaped hill around the peak. Skewed, squashed against a boundary, even ridge-shaped: the walker's record captures the joint uncertainty as it actually is.

Interactive: two portraits of the same summit

Now a two-parameter case where the difference matters. A small trial ($n=30$) models treatment response against a standardised biomarker score: $\text{logit}(p) = \beta_0 + \beta_1 x$. Both panels show the same log-likelihood surface. The left panel prices uncertainty the curvature way; on the right, a Metropolis walker builds the posterior the Bayesian way (flat prior). Set it walking, then compare the numbers below.

Highlighted cells mark where the two routes disagree. The posterior columns come from a fine grid evaluation of the posterior (240×240 points) — a “walk everywhere, carefully” version of what the Metropolis walker approximates; both are validated against scripts/py/generate_inference_data.py.

The marginal story for $\beta_1$

Choosing a route

Neither method is simply “better”. The curvature route is instant, deterministic and baked into every regression summary(); for large samples away from boundaries it is essentially exact. The sampling route costs computation and care ( convergence diagnostics, priors) but keeps working when the quadratic story breaks: small $n$, bounded parameters, skewed or multi-modal surfaces — and it extends naturally to quantities the curvature method struggles with, such as uncertainty about predictions passed through non-linear link functions.

A practical habit: run the cheap route always, the expensive route when it matters, and compare. Agreement is reassurance; disagreement is information.

In code

The Python output is a real run on docs/data/inference-logistic.json ($n=30$, seed 42). The R panel shows the standard workflow — rstanarm fits in seconds but is not run here, so no output is reproduced.

# Curvature route: glm() + Wald SEs from the inverse information
fit_ml <- glm(response ~ biomarker, family = binomial(), data = trial)
summary(fit_ml)          # coefficients, curvature-based Std. Errors
confint(fit_ml)          # profile-likelihood intervals: R's middle path —
                         # it re-walks the surface along each axis rather
                         # than trusting the quadratic

# Sampling route: rstanarm (Bayesian, weakly-informative default priors)
library(rstanarm)
fit_bayes <- stan_glm(response ~ biomarker, family = binomial(),
                      data = trial)
posterior_interval(fit_bayes, prob = 0.95)   # credible intervals
draws <- as.matrix(fit_bayes)   # the walker's record: one row per visit

import numpy as np
import statsmodels.api as sm

# Curvature route
fit = sm.Logit(y, sm.add_constant(x)).fit(disp=0)
print(fit.params)        # MLE
print(fit.bse)           # SEs from the inverse information matrix
print(fit.conf_int())    # Wald 95% intervals

# Sampling route (flat-prior grid posterior; PyMC/emcee for real problems)
# -- evaluate the log-likelihood over a fine (b0, b1) grid, exponentiate,
#    normalise, then read quantiles off the marginals.
# Full implementation: scripts/py/generate_inference_data.py

# MLE:                b0 = 0.9059,  b1 = 2.4329
# Wald SEs:           se(b0) = 0.6618, se(b1) = 0.8625
# Wald 95% CI:        b0: (-0.3911, 2.2030)   b1: (0.7425, 4.1233)
#
# Grid posterior (flat prior, 240x240):
#   b0: median 1.0555, sd 0.7473, 95% CrI (-0.2404, 2.7168)
#   b1: median 2.8876, sd 1.0606, 95% CrI ( 1.3720, 5.5291)

Further connections