Hierarchical Inference in Gaussian Models

Summary

The parameters $\Theta =\{m,k,A,b,\Lambda \}$ of a Gaussian model shape its posterior, but the “right” ones are usually unknown. Hierarchical inference places uncertainty on these hyperparameters and uses the marginal likelihood (evidence / type-II likelihood) to infer them. For the scalar mean and scale of a Gaussian, a conjugate prior (Gauss–Gamma; multivariate Gauss-inverse-Wishart) yields closed-form posterior updates and a Student-$t$ predictive. This is what lets a PN solver calibrate its uncertainty at runtime — cheaply estimating the one remaining scale so that its error bars are trustworthy.

Overview

A recurring PN thesis (Computation as Probabilistic Inference) is that classical methods arise from a family of Gaussian priors sharing one mean (the classical estimate) but differing in an overall uncertainty scale. The posterior standard deviation contracts at a rate tied to the classical convergence rate, differing only by that constant scale. If — and only if — that scale is calibrated well does the posterior variance become a trustworthy notion of numerical error. Hierarchical inference supplies the machinery to estimate it with minimal computational overhead, often during the solver’s run. Without it, PN uncertainty is arbitrary; with it, the error bars mean something.

Main Content

Hyperparameters and the evidence

In a Gaussian model $p(f\mid \Theta )=\mathcal{N}(f;m,k)$, $p(y\mid f,\Theta )=\mathcal{N}(y;Af+b,\Lambda )$, the numbers $\Theta$ that specify the model — mean $m$, kernel $k$ (e.g. output scale $\theta$, length-scale $\lambda$), noise $\Lambda$ — are hyperparameters. The distinction between a variable (of interest) and a parameter (nuisance) is conceptual, not formal; so if $\Theta$ is uncertain, assign it a distribution too. Integrating out $f$ gives the marginal likelihood (a.k.a. evidence or type-II likelihood) for $\Theta$:

$$p(y\mid \Theta )=\int p(f\mid \Theta )p(y\mid f,\Theta )\mathrm{d}f=\mathcal{N}(y;Am+b,k+\Lambda ).$$

Because this “models the model,” it is called hierarchical inference.

Hierarchical inference generally raises cost and yields a non-Gaussian (nonlinear) optimisation/inference problem. But for special cases — notably the scale and mean of a Gaussian — analytic forms exist via conjugate priors: priors whose posterior has the same functional form (only updated hyperparameters).

Scalar case: Gauss–Gamma conjugate prior

Gauss–Gamma prior and posterior

Let $y=[y_i{]}_{i=1,\dots ,N}$ be i.i.d. from $\mathcal{N}(\alpha ,{\beta }^2)$ with unknown mean $\alpha$ and variance ${\beta }^2$, so $p(y\mid \alpha ,\beta )=\mathcal{N}(y;\alpha 1,{\beta }^{2}I)$. The conjugate Gauss–Gamma prior (hyperparameters ${\theta }_0=({\mu }_0,{\lambda }_0,a_0,b_0)$) puts a Gaussian on $\alpha$ (variance $\propto {\beta }^2$) and a Gamma on ${\beta }^{-2}$:

$$\eta (\alpha ,\beta \mid {\mu }_0,{\lambda }_0,a_0,b_0)=\mathcal{N}(\alpha ;{\mu }_0,\frac{{\beta }^2}{{\lambda }_0})\mathcal{G}({\beta }^{-2};a_0,b_0),\mathcal{G}(z;a,b)=\frac{b^a{z}^{a-1}}{\Gamma (a)}{e}^{-bz}.$$

Multiplying by the likelihood and simplifying via the sufficient statistics — sample mean $\bar{\alpha }=\frac{1}{N}{\sum }_i{y}_i$ (Eq. 6.5) and sample variance ${\bar{\beta }}^2=\frac{1}{N}{\sum }_i(y_i-\bar{\alpha }{)}^2$ (Eq. 6.6) — the posterior is again Gauss–Gamma, $\eta (\alpha ,\beta \mid {\mu }_N,{\lambda }_N,a_N,b_N)$, with

$${\mu }_N=\frac{{\lambda }_0{\mu }_0+N\bar{\alpha }}{{\lambda }_0+N},{\lambda }_N={\lambda }_0+N,a_N=a_0+\frac{N}{2},b_N=b_0+\frac{1}{2}(N{\bar{\beta }}^2+\frac{{\lambda }_{0}N}{{\lambda }_0+N}(\bar{\alpha }-{\mu }_0{)}^2).\text{\hspace{1em}(6.8)}$$

Interpretation: ${\lambda }_N$ and $a_N$ act as counters for how many observations have been collected (with ${\lambda }_0,2a_0$ the “pseudo-observations” encoded in the prior); ${\mu }_N,b_N$ are sufficient statistics for the population mean and variance. As $N\to \infty$, ${\mu }_N\to \bar{\alpha }$ (sample mean) and $\mathbb{E}({\beta }^2)=b_N/a_N\to$ the sample variance plus the correction $\frac{{\lambda }_{0}N}{{\lambda }_0+N}(\bar{\alpha }-{\mu }_0{)}^2$, which corrects for the downward bias of ${\bar{\beta }}^2$.

Student- $t$ predictive

The predictive over the next sample $y_{N+1}$, marginalising the posterior over $(\alpha ,\beta )$, is a Student-$t$ distribution:

$$p(y_{N+1}\mid y)=\int \int p(y_{N+1}\mid \alpha ,\beta )p(\alpha ,\beta \mid y)\mathrm{d}\alpha \mathrm{d}\beta =\mathrm{St}(y_{N+1};{\mu }_N,\frac{b_N}{a_N},2a_N).$$

(Eq. 6.9.) Heavier tails than a Gaussian reflect the residual uncertainty about the variance ${\beta }^2$.

Multivariate case

For $Y=[y_1,\dots ,y_N]\in {\mathbb{R}}^{M\times N}$ i.i.d. from $\mathcal{N}(\alpha ,B)$ with unknown mean vector $\alpha$ and covariance $B$, the conjugate prior is the Gauss-inverse-Wishart, $p(\alpha ,B)=\mathcal{N}(\alpha ;{\mu }_0,\frac{1}{{\lambda }_0}B)\mathcal{W}(B^{-1};W_0,{\nu }_0)$. The posterior is again Gauss-inverse-Wishart with ${\mu }_N,{\lambda }_N,{\nu }_N={\nu }_0+N$ and $W_N=W_0+N\bar{B}+\frac{{\lambda }_{0}N}{{\lambda }_0+N}(\bar{\alpha }-{\mu }_0)(\bar{\alpha }-{\mu }_0{)}^{\top }$, and the predictive is a multivariate Student-$t$ (Eq. 6.12).

Hierarchical inference inside filters (the PN payoff)

Runtime scale calibration during filtering

For a linear-Gaussian state-space model, the marginal likelihood factorises over the Markov chain (prediction-error decomposition):

$$p(y\mid \theta )=\prod _{i=1}^{N}p(y_i\mid y_{1:i-1},\theta )=\prod _{i=1}^N\mathcal{N}(y_i;H{m}_i^{-},{\theta }^{2}H{\tilde{P}}_i^{-}{H}^{\top }),$$

where the output scale $\theta$ (from the SDE’s $L=\theta \tilde{L}$, see Gauss-Markov Processes and SDEs) is factored out using a unit-scale predictive covariance ${\tilde{P}}_i^{-}$ (defined via ${\tilde{Q}}_{t_i}={\int }_0^h{e}^{F\tau }\tilde{L}{\tilde{L}}^{\top }{e}^{F^{\top }\tau }\mathrm{d}\tau$). With an inverse-Gamma prior $p({\theta }^{-2})=\mathcal{G}({\theta }^{-2};{\alpha }_0,{\beta }_0)$, the posterior on $\theta$ updates recursively during the filtering pass:

$$p({\theta }^{-2}\mid y_{1:N})=\mathcal{G}({\theta }^{-2};{\alpha }_0+\frac{N}{2},{\beta }_0+\frac{1}{2}\sum _{i=1}^N\frac{(y_i-H{m}_i^{-}{)}^2}{H{\tilde{P}}_i^{-}{H}^{\top }}).\text{\hspace{1em}(6.14)}$$

In Filter (Alg. 5.1) this is realised by using $\tilde{Q}$ instead of $Q$ and adding two scalar updates per step (initialised $\alpha \leftarrow {\alpha }_0,\beta \leftarrow {\beta }_0$). The corresponding Student-$t$ marginal on the state can be computed locally each step.

This is exactly why “runtime hyperparameter calibration matters for PN uncertainty”: the output scale $\theta$ — the single free constant that turns a contraction rate into an absolute error bar — is estimated on the fly with negligible ($\mathcal{O}(1)$ per step) overhead, keeping the solver’s uncertainty well-calibrated.

Empirical Bayes vs full Bayes

• Empirical Bayes (type-II ML): pick a point estimate $\hat{\theta }=\arg {\max }_{\theta }p(y\mid \theta )$ (maximise the evidence) and plug it in. Cheap; ignores hyperparameter uncertainty.
• Full Bayes: integrate over $\theta$ (the predictive Student-$t$ above), propagating hyperparameter uncertainty into the final answer. More faithful; here tractable thanks to conjugacy. One could iterate hierarchical layers ad infinitum, but this is not done in practice to keep cost finite.

Examples

Scaled-observation variant (Exercise 6.1)

If instead of $y_i$ one observes individually scaled samples $a_i=y_i{s}_i$ with known scales $s_i$, the same Gauss–Gamma machinery applies with modified sufficient statistics $\bar{\alpha }=\frac{1}{N}{\sum }_i{a}_i/s_i$ and ${\bar{\beta }}^2=\frac{1}{N}{\sum }_i(a_i/s_i-\bar{\alpha }{)}^2$. This heteroscedastic variant is exactly what appears when calibrating linear-solver uncertainty (§III).

Calibrating an ODE filter's error bars

An ODE filter with an IWP prior has an unknown diffusion scale $\theta$ governing how fast its uncertainty grows between steps. Running the recursive Gamma update (Eq. 6.14) alongside the Kalman filter yields a data-driven $\theta$, so the solver’s reported credible intervals contract at the right absolute magnitude — the difference between a decorative error bar and a usable one.

Connections

• Infers the parameters of the Gaussians defined in Gaussian Distributions and Algebra, Gaussian Process Regression, and Gauss-Markov Processes and SDEs.
• The recursive-scale-in-filtering scheme feeds directly into ODE Filters and Smoothers and Theory of ODE Filters and Smoothers.
• Calibration of the remaining scale is applied in Uncertainty Calibration for Linear Solvers and Convergence and Priors in Bayesian Quadrature.

See Also

• Gaussian Process Regression — the model whose kernel hyperparameters are inferred here.
• Bayesian Filtering and Smoothing — where runtime scale calibration is embedded (Eq. 6.14).
• Uncertainty Calibration for Linear Solvers — the same idea for linear-algebra solvers.
• Computation as Probabilistic Inference — why calibrated scale makes PN error bars meaningful.