Marginal Likelihood via the Kalman Filter

Summary

Running the Kalman filter not only estimates the state — it also evaluates the marginal likelihood $p({\mathbf{y}}_{1:T}\mid \mathbf{\theta })$ of the model parameters $\mathbf{\theta }$, for free, as a by-product. The trick is the prediction-error decomposition: the joint density of the data factors into a product of one-step predictive densities, each a Gaussian in the filter’s innovation ${\mathbf{v}}_k$ with covariance ${\mathbf{S}}_k$. This turns an intractable high-dimensional integral over states into a simple recursive sum — which is precisely what lets you put a prior on $\mathbf{\theta }$ and run MCMC, MAP optimization, or EM. It is the engine behind parameter inference in BSTS.

Overview

To learn parameters $\mathbf{\theta }$ (the entries of $\mathbf{A},\mathbf{H},\mathbf{Q},\mathbf{R},{\mathbf{m}}_0,{\mathbf{P}}_0$) we want the marginal posterior $p(\mathbf{\theta }\mid {\mathbf{y}}_{1:T})\propto p({\mathbf{y}}_{1:T}\mid \mathbf{\theta })p(\mathbf{\theta })$. The states must be integrated out:

$$p(\mathbf{\theta }\mid {\mathbf{y}}_{1:T})=\int p({\mathbf{x}}_{0:T},\mathbf{\theta }\mid {\mathbf{y}}_{1:T})\mathrm{d}{\mathbf{x}}_{0:T}.$$

Doing this integral directly is intractable (and grows with $T$). The state-space structure dissolves it.

Main Content

Theorem: Prediction-error decomposition (Särkkä Thm. 12.1, Eq. 12.5)

The marginal likelihood factors into one-step predictive densities:

$$p({\mathbf{y}}_{1:T}\mid \mathbf{\theta })=\prod _{k=1}^{T}p({\mathbf{y}}_k\mid {\mathbf{y}}_{1:k-1},\mathbf{\theta }),\text{\hspace{1em}(12.5)}$$

where each term is obtained by integrating the measurement model against the predicted state:

$$p({\mathbf{y}}_k\mid {\mathbf{y}}_{1:k-1},\mathbf{\theta })=\int p({\mathbf{y}}_k\mid {\mathbf{x}}_k,\mathbf{\theta })p({\mathbf{x}}_k\mid {\mathbf{y}}_{1:k-1},\mathbf{\theta })\mathrm{d}{\mathbf{x}}_k.\text{\hspace{1em}(12.6)}$$

The predictive $p({\mathbf{x}}_k\mid {\mathbf{y}}_{1:k-1},\mathbf{\theta })$ is exactly the Kalman prediction step, and the normalizer $Z_k$ of the update step is this term — so the filter computes it at no extra cost.

Theorem: Energy function for the linear-Gaussian model (Särkkä Thm. 12.3, Eq. 12.38)

Define the energy function (unnormalized negative log-posterior) ${\phi }_T(\mathbf{\theta })=-\log p({\mathbf{y}}_{1:T}\mid \mathbf{\theta })-\log p(\mathbf{\theta })$. For the linear-Gaussian model it is built recursively alongside the Kalman filter:

$${\phi }_k(\mathbf{\theta })={\phi }_{k-1}(\mathbf{\theta })+\frac{1}{2}\log |2\pi {\mathbf{S}}_k(\mathbf{\theta })|+\frac{1}{2}{\mathbf{v}}_k^{\mathsf{T}}(\mathbf{\theta }){\mathbf{S}}_k^{-1}(\mathbf{\theta }){\mathbf{v}}_k(\mathbf{\theta }),\text{\hspace{1em}(12.38)}$$

started from ${\phi }_0(\mathbf{\theta })=-\log p(\mathbf{\theta })$, where ${\mathbf{v}}_k$ (innovation) and ${\mathbf{S}}_k$ (innovation covariance) come straight from the Kalman update step. Each one-step predictive density is the Gaussian

$$p({\mathbf{y}}_k\mid {\mathbf{y}}_{1:k-1},\mathbf{\theta })=\mathcal{N}({\mathbf{y}}_k\mid {\mathbf{H}}_k{\mathbf{m}}_k^{-},{\mathbf{S}}_k),{\mathbf{S}}_k={\mathbf{H}}_k{\mathbf{P}}_k^{-}{\mathbf{H}}_k^{\mathsf{T}}+{\mathbf{R}}_k.\text{\hspace{1em}(12.41)}$$

The marginal parameter posterior is then $p(\mathbf{\theta }\mid {\mathbf{y}}_{1:T})\propto \exp (-{\phi }_T(\mathbf{\theta }))$.

Why this makes inference tractable

• The intractable integral is gone. The full likelihood is a recursive sum of cheap Gaussian terms — one Kalman pass evaluates ${\phi }_T(\mathbf{\theta })$ exactly for any $\mathbf{\theta }$.
• MAP / ML estimates: ${\hat{\mathbf{\theta }}}^{\mathrm{MAP}}=\arg {\min }_{\mathbf{\theta }}{\phi }_T(\mathbf{\theta })$ (Eq. 12.13); the ML estimate is the same with a flat prior $p(\mathbf{\theta })\propto 1$.
• MCMC: Metropolis–Hastings only needs the unnormalized posterior, i.e. ${\phi }_T(\mathbf{\theta })$; the normalizing constant $p({\mathbf{y}}_{1:T})$ is never required. The acceptance ratio uses energy differences directly: $${\alpha }_i=\min \{1,\exp ({\phi }_T({\mathbf{\theta }}^{(i-1)})-{\phi }_T({\mathbf{\theta }}^{\ast })\left)\frac{q({\mathbf{\theta }}^{(i-1)}\mid {\mathbf{\theta }}^{\ast })}{q({\mathbf{\theta }}^{\ast }\mid {\mathbf{\theta }}^{(i-1)})}\right\}.\text{\hspace{1em}(12.17)}$$ So a single Kalman-filter run per proposal gives a full MCMC sampler over model parameters.
• EM: Fisher’s identity expresses $\nabla {\phi }_T$ as an expectation of the complete-data score under the smoothing distribution — the RTS smoother supplies the required sufficient statistics (Särkkä §12.3, Eq. 12.32).
• Laplace approximation: $p(\mathbf{\theta }\mid {\mathbf{y}}_{1:T})\approx \mathcal{N}(\mathbf{\theta }\mid {\hat{\mathbf{\theta }}}^{\mathrm{MAP}},[\mathbf{H}({\hat{\mathbf{\theta }}}^{\mathrm{MAP}}){]}^{-1})$ using the Hessian of ${\phi }_T$ (Eq. 12.15).

Examples

Posterior of the noise variance in the random walk (Särkkä Ex. 12.1)

For the Gaussian random-walk model, fix the process variance and treat the measurement variance $R$ as the unknown $\theta$. With a flat prior $p(R)\propto 1$, running the scalar Kalman filter and accumulating ${\phi }_T(R)$ from Eq. (12.38) yields $p(R\mid y_{1:T})\propto \exp (-{\phi }_T(R))$. Särkkä’s Fig. 12.1 shows the true $R$ sits in the high-density region — but the MAP/ML point estimate is biased low relative to the truth, illustrating why full posterior (MCMC) treatment beats point estimation.

Worked energy step (local level model)

Per time step, with scalar innovation $v_k=y_k-m_k^{-}$ and $S_k=P_k^{-}+R$:

$${\phi }_k={\phi }_{k-1}+\frac{1}{2}\log (2\pi S_k)+\frac{v_k^2}{2S_k}.$$

Summing over $k$ and adding $-\log p(\theta )$ gives the negative log-posterior used inside an MCMC loop.

Connections

• The Kalman Filter — supplies the innovations ${\mathbf{v}}_k$ and covariances ${\mathbf{S}}_k$ that form every likelihood term
• The RTS Smoother — supplies smoothing statistics for EM / Fisher’s-identity gradients
• Linear-Gaussian State-Space Models — the parameters $\mathbf{\theta }=(\mathbf{A},\mathbf{H},\mathbf{Q},\mathbf{R},\dots )$ being inferred
• State-Space Models and the Kalman Filter - Overview — pipeline context

See Also

• Bayesian Structural Time-Series Model — BSTS marginal likelihood / sampler is built on exactly this prediction-error decomposition
• MCMC Inference for CausalImpact — the Gibbs sampler draws variances/coefficients conditional on states, the complement to integrating states out here
• Spike-and-Slab Prior for Covariate Selection — variable selection relies on tractable marginal likelihoods of competing models
• Single Marketing Time Series — ARIMA likelihoods are evaluated the same way via the Kalman filter