Fast Bayesian inference of the multivariate Ornstein-Uhlenbeck process

The multivariate Ornstein-Uhlenbeck process is used in many branches of science and engineering to describe the regression of a system to its stationary mean. Here we present an $O(N)$ Bayesian method to estimate the drift and diffusion matrices of the process from $N$ discrete observations of a sample path. We use exact likelihoods, expressed in terms of four sufficient statistic matrices, to derive explicit maximum a posteriori parameter estimates and their standard errors. We apply the method to the Brownian harmonic oscillator, a bivariate Ornstein-Uhlenbeck process, to jointly estimate its mass, damping, and stiffness and to provide Bayesian estimates of the correlation functions and power spectral densities. We present a Bayesian model comparison procedure, embodying Ockham’s razor, to guide a data-driven choice between the Kramers and Smoluchowski limits of the oscillator. These provide novel methods of analyzing the inertial motion of colloidal particles in optical traps.

💡 Research Summary

The paper introduces a fast, exact Bayesian framework for inferring the drift (λ) and diffusion (σ, or equivalently D) matrices of a multivariate Ornstein‑Uhlenbeck (OU) process from a discrete time series of length N. The OU process is the unique continuous stochastic process that is simultaneously Gaussian, stationary, and Markovian, which means its transition density and stationary distribution are both multivariate normal and can be written in closed form.

Key methodological contributions:

1. Sufficient‑statistics reduction – By exploiting the Gaussian‑Markov structure, the full log‑posterior can be expressed solely in terms of four M × M matrices (T₁–T₄) that are simple averages of products of successive observations. T₁ and T₂ involve one‑step‑ahead covariances, T₃ is the empirical covariance of the current state, and T₄ is the outer product of the initial observation. Computing these matrices requires only a single pass over the data, i.e., O(N) operations, and they can be updated incrementally for online estimation.

2. Closed‑form MAP estimates – The maximum‑a‑posteriori (MAP) estimates for the exponential of the drift matrix (Λ = e^{‑λΔt}) and the stationary covariance Σ are obtained analytically:
   Λ* = T₂ T₃⁻¹, Σ* = (1/N)