Fast Bayesian inference of the multivariate Ornstein-Uhlenbeck process

TL;DR

This paper introduces an efficient Bayesian approach for estimating parameters of the multivariate Ornstein-Uhlenbeck process from discrete data, enabling detailed analysis of systems like colloidal particles in optical traps.

Contribution

It presents an $O(N)$ Bayesian method using exact likelihoods for parameter estimation and model comparison in multivariate Ornstein-Uhlenbeck processes, including applications to physical systems.

Findings

Efficient $O(N)$ Bayesian estimation of drift and diffusion matrices.

Explicit maximum a posteriori estimates with standard errors.

Bayesian model comparison between Kramers and Smoluchowski limits.

Abstract

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.

Tables

Table 1. Table 1: Bayesian MAP estimates and standard errors of the parameters of the Brownian harmonic oscillator. There is excellent agreement between the two Bayesian methods and with the parameter values used to generate the paths.

$m\left(ng\right)$
Simulation	Bayes I	Bayes II
$1$	$0.994\pm 0.006$	$0.993\pm 0.008$
$1.5$	$1.521\pm 0.014$	$1.517\pm 0.012$
$2$	$2.107\pm 0.015$	$2.092\pm 0.016$

Equations

$$d{x}_i=-{\lambda }_{ij}{x}_{j}dt+\sigma {}_{ij}d{W}_j,$$

$$\mathcal{L}(\mathbf{x})=\frac{\partial }{\partial x_i}{\lambda }_{ij}{x}_j+\frac{1}{2}\frac{{\partial }^2}{\partial x_i\partial x_j}(\mathbf{\sigma }{\mathbf{\sigma }}^{\text{tr}}{)}_{ij}.$$

$${\mathbf{x}}^{\prime },t^{\prime }|\mathbf{x},t\sim \mathcal{N}\left(\mathbf{\mu },\mathbf{\Sigma }\right),$$

$$\mathbf{\mu }=\mathbf{\Lambda }\mathbf{x},\mathbf{\Sigma }=\mathbf{c}-\mathbf{\Lambda }\mathbf{c}{\mathbf{\Lambda }}^{\text{tr}},\mathbf{\Lambda }=e^{-\mathbf{\lambda }|\Delta t|}$$

$$\mathbf{x}\sim \mathcal{N}\left(0,\mathbf{c}\right).$$

$$\mathbf{\lambda c}+(\mathbf{\lambda c}{)}^{\text{tr}}=\mathbf{\sigma }{\mathbf{\sigma }}^{\text{tr}}.$$

$$\mathbf{c}={\left({\mathbf{\lambda }}^{\mathrm{ir}}\right)}^{-1}\left(\mathbf{\sigma }{\mathbf{\sigma }}^{\text{tr}}\right),$$

$$\mathbf{C}(t-t^{\prime })\equiv ⟨\mathbf{x}(t){\mathbf{x}}^{\text{tr}}(t^{\prime })⟩=e^{-\mathbf{\lambda }|\Delta t|}\mathbf{c},$$

$$\mathbf{C}(\Omega )=(\text{textminus}i\Omega 1+\mathbf{\lambda }{)}^{-1}\left(2\mathbf{D}\right)(i\Omega 1+{\mathbf{\lambda }}^{\text{tr}}{)}^{-1},$$

$$P\left(\mathbf{X}|\mathbf{\theta }\right)=\prod _{n=1}^{N-1}{P}_{1|1}\left({\mathbf{x}}_{n+1}|{\mathbf{x}}_n,\mathbf{\theta }\right)P_1\left({\mathbf{x}}_1|\mathbf{\theta }\right).$$

$$P\left(\mathbf{\theta }|\mathbf{X}\right)=\frac{P\left(\mathbf{X}|\mathbf{\theta }\right)P\left(\mathbf{\theta }\right)}{P\left(\mathbf{X}\right)}.$$

$$\ln$$

$${\mathbf{a}}_0={\mathbf{\Sigma }}^{-1}+\mathbf{\Lambda }{}^{\text{tr}}{\mathbf{\Sigma }}^{-1}\mathbf{\Lambda },{\mathbf{a}}_1=-{\mathbf{\Sigma }}^{-1}\mathbf{\Lambda },{\mathbf{a}}_0^{\prime }={\mathbf{\Sigma }}^{-1}+{\mathbf{c}}^{-1},{\mathbf{a}}_0^{\prime \prime }={\mathbf{\Sigma }}^{-1}.$$

$$\ln P\left(\bm{\theta}|\bm{X}\right)=-\frac{\bm{\Sigma}^{-1}}{2}:\Big{[}\left(\bm{\Lambda}-\bm{T}_{2}\bm{T}_{3}^{-1}\right)\bm{T}_{3}\left(\bm{\Lambda}-\bm{T}_{2}\bm{T}_{3}^{-1}\right)^{\text{tr}}+\left(\bm{T}_{1}-\bm{T}_{2}\bm{T}_{3}^{-1}\bm{T}_{2}^{\text{tr}}\right)\Big{]}-\tfrac{N-1}{2}\ln|2\pi\boldsymbol{\Sigma}|-\tfrac{\bm{c}^{-1}}{2}:\boldsymbol{T}_{4}-\tfrac{1}{2}\ln|2\pi\boldsymbol{c}|.$$

$${\mathbf{T}}_1$$

$${\mathbf{T}}_3$$

$${\mathbf{\Lambda }}^{\ast }$$

$${\mathbf{\Sigma }}^{\ast }$$

$${\mathbf{c}}^{\ast }={\mathbf{\Sigma }}^{\ast }{\left[1-({\mathbf{\Lambda }}^{\ast }\mathbf{ϵ}{)}^{\text{tr}}({\mathbf{\Lambda }}^{\ast }\mathbf{ϵ}{)}^{\text{tr}}\right]}^{-1}.$$

$$\ln P\left(\mathbf{c}|\mathbf{X}\right)$$

$${\mathbf{c}}^{\ast }=\frac{1}{N}\sum _{n=1}^N{\mathbf{x}}_n{\mathbf{x}}_n^{\text{tr}}=\frac{1}{N}{\mathbf{T}}_3,$$

$$P({\mathcal{M}}_{\alpha }|\mathbf{X})\propto P(\mathbf{X}|{\mathcal{M}}_{\alpha })P({\mathcal{M}}_{\alpha }),$$

$$P(\mathbf{X}|{\mathcal{M}}_{\alpha })=\int P(\mathbf{X}|\mathbf{\theta },{\mathcal{M}}_{\alpha })P(\mathbf{\theta }|{\mathcal{M}}_{\alpha })d\mathbf{\theta }.$$

$$P(\mathbf{X}|{\mathcal{M}}_{\alpha })\simeq P(\mathbf{X}|{\mathbf{\theta }}^{\ast },{\mathcal{M}}_{\alpha })P({\mathbf{\theta }}^{\ast }|{\mathcal{M}}_{\alpha })\Delta \mathbf{\theta }.$$

$${\mathbf{x}}_{n+1}=\mathbf{\Lambda }{\mathbf{x}}_n+\sqrt{\mathbf{\Sigma }}{\mathbf{\xi }}_n,$$

$$m\dot{v}+\gamma v+\nabla U(x)=\xi .$$

$$dv=\left(-{\omega }_0^{2}x-v/\tau \right)dt+\sigma d{W}_v,dx=vdt.$$

$$d\left(\begin{array}[]{c}x\\ v\end{array}\right)=-\bm{\lambda}\left(\begin{array}[]{c}x\\ v\end{array}\right)dt+\boldsymbol{\sigma}\left(\begin{array}[]{c}dW_{x}\\ dW_{v}\end{array}\right),$$

$$\mathbf{\lambda }=\left(\begin{array}{cc}0& -1\\ {\omega }_0^2& 1/\tau \end{array}\right),$$

$$\mathbf{\sigma }=\left(\begin{array}{cc}0& 0\\ 0& \sqrt{2D}/\tau \end{array}\right).$$

$$\mathbf{c}=\left(\begin{array}{cc}{k}_{B}T/k& 0\\ 0& k_{B}T/m\end{array}\right).$$

Full text

DOI: 10.1103/PhysRevE.98.012136

I Introduction

The multivariate Ornstein-Uhlenbeck process is widely used in many branches of science and engineering to describe the regression of a system to its stationary mean. Its importance arises from the fact that it is the only continuous stochastic process that is simultaneously stationary, Gaussian and Markovian doob1942brownian ; doob1953stochastic . Therefore, the process is fully characterized by its stationary and conditional distributions, each of which is a multivariate Gaussian, with mean and variance determined by the regression rates and diffusion coefficients of the process van1992stochastic ; gardiner1985handbook . This simplicity implies that the likelihood associated with discretely observed sample paths can be calculated explicitly and exactly in terms of the process parameters and that the posterior distributions of the latter can be obtained, without approximation, from discrete observations.

Here we present the details of such an exact Bayesian estimation method for the parameters of the $M$-dimensional Ornstein-Uhlenbeck process. We show that the likelihood can be expressed in terms of four matrices whose elements are the self- and mutual-correlations of the stochastic variables at equal times and one observation time apart. We derive the maximum a posteriori (MAP) estimates and the error bars of the parameters in terms of these four sufficient statistics. We address the problem of model selection, i. e., of choosing the Ornstein-Uhlenbeck model with the least number of parameters that best explains the data, within the framework of a Bayesian model comparison.

The principal computational cost of our method is in evaluating the four sufficient statistic matrices. For a time series of $N$ discrete observations, only $O(N)$ operations are required to compute the necessary self- and mutual-correlations. Further, additional $N^{\prime}$ observations can be incorporated incrementally and the sufficient statistics recomputed at a cost of $O(N^{\prime})$, making our algorithm suitable for “online” estimation. Together with the consistency and optimality of Bayesian estimators, this yields a fast and accurate parameter estimation method for the multivariate Ornstein-Uhlenbeck process.

We apply our general results to the Brownian harmonic oscillator, a bivariate Ornstein-Uhlenbeck process that is ubiquitous in the physical sciences. This provides a Bayesian method for jointly estimating the mass, friction, and spring constant of the oscillator and new ways of estimating its correlation function and power spectral density. The results are pertinent for the analysis of time-series data of inertial motion of colloidal particles in optical traps and for selecting between the underdamped and overdamped models of oscillator motion, when the parameters are not known *apriori. *

We make a few remarks on how this work relates to the preceding literature. Bayesian inference for general Gaussian processes, of which the Ornstein-Uhlenbeck process is a special case, has been an area of extensive research in the past few decades mackay2003information ; rasmussen2006gaussian ; lindgren2011explicit . For $N$ observations of the process, inference requires the inversion of an $N\times N$ symmetric covariance matrix whose computational cost, using a direct method, is $O(N^{3})$. This makes inference superlinear in $N$ and for large $N$ approximate methods become necessary smola2000sparse ; williams2001using ; csato2002sparse ; quinonero2005unifying . Imposing the Markov property on a general Gaussian process reduces the covariance matrix to symmetric tridiagonal form and, imposing stationarity further simplifies it to symmetric, tridiagonal and Toeplitz form. In this limit, the inverse can be obtained analytically in terms of sufficient statistics, and the only cost is in computing the latter, which is proportional to the number of data points $N$. This special feature of the stationary Gauss-Markov processes, which allows for fast, yet exact, inference, appears not to have been exploited earlier.

The remainder of the paper is organized as follows. In section II, we briefly review key properties of the Ornstein-Uhlenbeck process, and we make use of them in section III to obtain MAP estimates for the parameters and the model odds. In section IV, we present an exact path-sampling algorithm, which is used in V to generate sample paths of the Brownian harmonic oscillator and to validate the results of section III. We conclude in section VI with a discussion on further applications.

II Multivariate Ornstein-Uhlenbeck Process

The multivariate Ornstein-Uhlenbeck is defined by the It$\bar{\text{o}}$ stochastic differential equation gardiner1985handbook

[TABLE]

where $-\lambda_{ij}$ is a stable matrix of mean regression rates, $\sigma{}_{ij}$ is the volatility matrix, $W_{i}(t)$ are Wiener processes, and $i,j=1,\ldots,M$. We denote $(x_{1},\ldots,x_{M})^{\text{tr}}$ by the vector $\boldsymbol{x}$ and $\lambda_{ij}$ by the matrix $\boldsymbol{\lambda}$, with similar bold-face notation for other vectors and matrices, when convenient. The volatility matrix is related to $\boldsymbol{D}$, a symmetric positive-semi-definite matrix of diffusion coefficients, as $\boldsymbol{\sigma}\boldsymbol{\sigma}^{\text{tr}}\equiv 2\boldsymbol{D}$.

The probability density of a displacement from $\boldsymbol{x}$ at time $t$ to $\boldsymbol{x}^{\prime}$ at time $t^{\prime}$, P_{1|1}$$(\boldsymbol{x}^{\prime},t^{\prime}|\boldsymbol{x},t), obeys the Fokker-Planck equation $\partial_{t}P_{1|1}=\mathcal{L}P_{1|1}$, where the Fokker-Planck operator is

[TABLE]

The solution is a multivariate normal distribution

[TABLE]

where

[TABLE]

and $\Delta t=t^{\prime}-t$. This solution is exact and holds for arbitrary values of $\Delta t$. The stationary distribution $P_{1}(\boldsymbol{x})$ obeys the steady-state Fokker-Planck equation $\mathcal{L}P_{1}=0$, and the solution is, again, a normal distribution,

[TABLE]

Then, $\boldsymbol{c}=\langle\boldsymbol{x}\boldsymbol{x}^{\text{tr}}\rangle$ can be identified as the matrix of covariances in the stationary state, and $\mathcal{L}P_{1}=0$ implies that the matrices $\boldsymbol{\lambda}$, $\boldsymbol{c}$, and $\boldsymbol{\sigma}$ are not all independent but are related by the stationarity condition

[TABLE]

This is a Lyapunov matrix equation for $\boldsymbol{c},$ given $\boldsymbol{\lambda}$ and $\boldsymbol{\sigma}$. Solutions are considerably simplified when the Fokker-Planck operator obeys detailed balance $\mathcal{L}(\boldsymbol{x)}P_{1}(\boldsymbol{x})=P_{1}(\boldsymbol{\epsilon x})\mathcal{L}^{\dagger}(\boldsymbol{\epsilon x})$, where $\mathcal{L^{\dagger}}$ is the adjoint Fokker-Planck operator, \text{\boldsymbol{\epsilon} is a diagonal matrix of the parities }\epsilon_{i}=\pm 1 of $x_{i}$ under time reversal, and the stationary distribution is time-reversal invariant, $P(\boldsymbol{x})=P(\boldsymbol{\epsilon x})$. This implies Onsager-Casimir symmetry $\boldsymbol{\epsilon}(\boldsymbol{\lambda c})=(\boldsymbol{\lambda c})^{\text{tr}}\boldsymbol{\epsilon}$ for the regression matrix and $\boldsymbol{\epsilon c}=\boldsymbol{c\epsilon}$ for the covariance matrix. The matrix of covariances is then determined as

[TABLE]

where $\lambda_{ij}^{\mathrm{ir}}=(\lambda_{ij}+\epsilon_{i}\epsilon_{j}\lambda_{ij})$ .

The Gauss-Markov property of the Ornstein-Uhlenbeck process ensures that the correlation function

[TABLE]

decays exponentially and that its Fourier transform, the power spectral density

[TABLE]

is a multivariate Lorentzian in the angular frequency $\Omega$ van1992stochastic .

In what follows, we shall take $\boldsymbol{\Lambda}$, the matrix exponential of the mean regressions rates and $\boldsymbol{c}$, the covariance matrix, to be the independent parameters. Estimates of the parameters in the diffusion matrix can then be obtained from the estimates of $\boldsymbol{\Lambda}$ and $\boldsymbol{c}$ through the stationarity condition. Thus, there are $M^{2}+M$ independent parameters for a $M$-variate Ornstein-Uhlenbeck process. For notational brevity, the set of all unknown parameters is collected in $\boldsymbol{\theta}=(\boldsymbol{\Lambda},\boldsymbol{c})$.

III Bayesian inference

*Parameter Estimation: *Consider now the discrete time series \boldsymbol{X}=\big{\{}\boldsymbol{x}_{1},\boldsymbol{x}_{2},...,\boldsymbol{x}_{N}\big{\}}, consisting of $N$ observations of the sample path $\boldsymbol{x}(t)$ at the discrete times $t=n\Delta t$ with $n=1,\ldots,N.$ Each observation $\boldsymbol{x}_{n}$ is an $M$-dimensional vector corresponding to the number of components of the multivariate Ornstein-Uhlenbeck process. From the Markov property of the process, the probability of the path, given the parameters in $\boldsymbol{\theta}$, is

[TABLE]

The probability $P\left(\boldsymbol{\theta}|\boldsymbol{X}\right)$ of the parameters, given the sample path, is given by Bayes theorem to be

[TABLE]

The denominator $P\left(\bm{X}\right)$ is a normalization independent of the parameters, and thus it can be ignored in parameter estimation. Using informative uniform priors for $P\left(\boldsymbol{\theta}\right)$, the logarithm of the posterior probability, after using the explicit forms of $P_{1|1}$ and $P_{1}$, is in matrix form

[TABLE]

where

[TABLE]

We note that the above matrix of covariances is real, symmetric, tridiagonal, and almost-Toeplitz (only the first and last elements of the diagonal differ from the remaining elements). Symmetry follows from the Gaussian property rasmussen2006gaussian , tridiagonality from the Markov property kavcic2000matrices , and the Toeplitz character from stationarity van1992stochastic . Exploiting these properties, with some elementary manipulations detailed in Appendix A, the posterior probability can be written in terms of the four matrix sufficient statistics as

[TABLE]

The sufficient statistics matrices are

[TABLE]

We note that these are matrices of covariances at equal times and one observation time apart. Therefore, they can be evaluated at a cost proportional to $2N-1$. It is this property that is results in our method being “fast”, that is, requiring only $O(N)$ computational cost.

The MAP estimates for mean regression rates and the covariances turn out to be

[TABLE]

The standard errors are obtained from the Hessian matrix $\bm{A}$ of the posterior probability evaluated at the maximum jeffreys1998theory ; jaynes2003probability ; sivia2006data . Its explicit form is provided in the Appendix B. The four $M\times M$ sufficient statistic matrices $\boldsymbol{T}_{i}$ rather than the considerably larger $M\times N$ times series $\boldsymbol{X}$ contains all the information relevant for inference. Their use reduces both the storage and computational cost of the inference algorithm.

The estimate for the mean regression rate is recognized to be the multivariate generalization of our earlier result for the univariate Ornstein-Uhlenbeck process bera2017fast . An explicit estimate for $\boldsymbol{c}$ can be obtained from that of $\boldsymbol{\Sigma}$ when the Fokker-Planck operator obeys detailed balance. Then Onsager-Casimir symmetry $\boldsymbol{\epsilon}(\boldsymbol{\lambda c})=(\boldsymbol{\lambda c})^{\text{tr}}\boldsymbol{\epsilon}$ implies $\boldsymbol{\Lambda}\boldsymbol{c}=\boldsymbol{c}\boldsymbol{\epsilon}\boldsymbol{\Lambda}^{\text{tr}}\boldsymbol{\epsilon}$ and, from the definition of $\boldsymbol{\Sigma}$, it follows that the MAP estimate for $\boldsymbol{c}$ is

[TABLE]

The MAP estimate for the matrix of diffusion coefficients then follows from the stationarity condition. This is the multivariate generalization of our earlier result for the diffusion coefficient for the univariate Ornstein-Uhlenbeck process. We refer to this method as “Bayes I”, following the terminology of bera2017fast .

In the absence of detailed balance, such explicit expressions can no longer be found and linear systems have to be solved to obtain the MAP estimate of $\boldsymbol{c}$ from that of $\boldsymbol{\Sigma}$ and, then, to relate it to the matrix of diffusion coefficients. We shall pursue this elsewhere.

An alternative Bayesian procedure for directly estimating $\boldsymbol{c}$ is arrived at by interpreting the time series $\boldsymbol{X}$ as independent repeated samples from the stationary distribution $\mathcal{N}\left(\boldsymbol{0},\boldsymbol{c}\right)$. Using non-informative priors, the expression of the logarithm of the posterior probability in this approach is given as

[TABLE]

from which the MAP estimate

[TABLE]

follows straightforwardly. The Bayesian inference method described above, using the stationary distribution of the Ornstein-Uhlenbeck process, is referred to as “Bayes II”.

Bayes I exploits both the Gaussian and Markovian character of the process, while Bayes II exploits only its Gaussian character. Therefore, numerical agreement between the above two methods of estimating the covariance matrix provides a stringent test of the stationary, Gaussian, and Markovian characters of the process. Since the Ornstein-Uhlenbeck process is the only continuous process with all of these properties, agreement affirms it as the data generating model. The preceding equations (23-29) are the main results of this paper.

*Model comparison: *Thus far we assumed that the data generating model was given and that only the parameters of the model needed to be estimated. In certain circumstances, though, the model itself may be uncertain and it becomes necessary to estimate the probability of different models $\mathcal{M}_{\alpha}$ jeffreys1998theory ; kashyap1977bayesian ; mackay1992bayesian ; zellner1984 ; gregory2005bayesian . The probability of a model, given the data, is

[TABLE]

where the first term on the right is the “evidence” of the model and the second term is the prior probability of the model. We shall assume all models to be, a priori, equally likely. The evidence is the normalizing constant in Eq.(10), given as an integral over the space of parameters $\boldsymbol{\theta}$ contained in the drift and diffusion matrices:

[TABLE]

For unimodal posterior distributions, the height at the MAP value $\boldsymbol{\theta}^{\ast}$ times the width $\Delta\boldsymbol{\theta}$ of the distribution is, often, a very good approximation for the evidence,

[TABLE]

The first term is the best fit likelihood while the second term, the product of the prior for the MAP estimate and the standard error of this estimate is called the Ockham factor. Thus models that achieve a compromise between the degree of fit to the data and the number of parameters required for the fit are ones that are favoured by the Bayesian model selection procedure. This avoids the over-fitting that would occur if the degree of fit was made the sole criterion for model selection, and it encodes the commonsense “principle of parsimony”, attributed to William of Ockham, which states that between two models that fit the data equally well, the simpler one is to be preferred. Both the evidence and the Ockham factor can be obtained straightforwardly for the Ornstein-Uhlenbeck models, and we shall make use of this below for model selection within a family of Ornstein-Uhlenbeck models.

To summarize, we derive two $O(N)$ methods to compute the parameters of a multivariate Ornstein-Uhlenbeck process in terms of sufficient statistics. The cost of computing these statistics for an additional $N^{\prime}$ data points is only $O(N^{\prime})$. Thus previous computations can be reused with the arrival of fresh data, making the method suitable for “online” estimation. The fact that new data can be incrementally added to the sufficient statistics is also useful for comparisons within a family of $M$-variate Ornstein-Uhlenbeck models. For large $M$, the determinant of the Hessian matrix needed for model comparison can be computed efficiently using iterative Krylov subspace methods, thereby allowing for a full “online” inference of both models and their parameters.

IV Path sampling

The solution of the Fokker-Planck equation, Eq.(3), provides a method for sampling paths of the multivariate Ornstein-Uhlenbeck process exactly. Given an initial state $\boldsymbol{x}$ at time $t$ and final state $\boldsymbol{x}^{\prime}$ at time $t^{\prime}$ , the quantity $\boldsymbol{x}^{\prime}-\boldsymbol{\Lambda}\boldsymbol{x}$ is normally distributed with mean zero and variance $\boldsymbol{\Sigma}$, a property that was first recognized by Uhlenbeck and Ornstein. Therefore, a sequence of states at times $t=n\Delta t$, where $n$ is a positive integer, forming a discrete sampling of a path can be obtained from the following iteration uhlenbeck1930theory

[TABLE]

where $\sqrt{\boldsymbol{\Sigma}}$ is a matrix square-root of $\boldsymbol{\Sigma}$, and $\boldsymbol{\xi}_{n}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{1})$ is an $M$-dimensional uncorrelated normal variate with zero mean and unit variance. The exponential of the regression matrix and the square-root of the variance matrix are the two key quantities in the iteration. They can be obtained analytically in low dimensional problems, but for high-dimensional problems they will, in general, have to be obtained numerically parker2012sampling . The sampling interval must satisfy $\lambda_{\mathrm{max}}\Delta t\ll 1$ such that the shortest time scale in the dynamics, corresponding to the inverse of the largest eigenvalue $\lambda_{\mathrm{max}}$ of the regression matrix, is resolved in the samples. In the following section, we use the above method to sample paths of the Brownian harmonic oscillator, which is equivalent to a bivariate Ornstein-Uhlenbeck process, to test the accuracy of our Bayesian estimates.

V Brownian harmonic oscillator

*Inference problem: *We now apply the results above to the physically important case of a massive Brownian particle confined in a harmonic potential described by the Langevin equation

[TABLE]

Here the pair $(x,v)$ describes the state of the particle in its phase space of position and velocity, while $m$ and $\gamma$ are the particle mass and friction coefficient respectively. The potential $U=\frac{1}{2}kx^{2}$ is harmonic with a stiffness $k$. $\xi\left(t\right)$is a zero-mean Gaussian white noise with variance $\langle\xi\left(t\right)\xi\left(t^{\prime}\right)\rangle=2k_{B}T\gamma\delta\left(t-t^{\prime}\right)$ that satisfies the fluctuation-dissipation relation chandrasekhar1943stochastic . The inference problem is to jointly estimate the triplet of parameters $(m,\gamma,k)$ from discrete observations of the position and velocity and to estimate the correlation functions and the spectral densities from these observations.

*Bivariate Ornstein-Uhlenbeck process: *The Langevin equation can be recast as a bivariate Ornstein-Uhlenbeck process in phase space for the pair ($x,v$) as

[TABLE]

Here $\omega_{0}^{2}=k/m$ is the natural frequency of the undamped harmonic oscillator, $\tau=m/\gamma$ is the characteristic time scale associated with the thermalization of the momentum due to viscous dissipation, and $\sigma=\sqrt{2D}$, where the diffusion coefficient $D$ of the particle is defined, as usual, by the Einstein relation $D=k_{B}T\gamma^{-1}$ einstein1905theory ; kubo1966fluctuation . Thus, we have constructed three independent parameters from the four dependent parameters. The resulting bivariate system can be written in matrix form as

[TABLE]

where the mean regression matrix is

[TABLE]

and the volatility matrix is

[TABLE]

The structure of the volatility matrix ensures that the positional Wiener process $dW_{x}$ does not enter the dynamics.

At thermal equilibrium, the joint distribution of position and velocity factorize into the Gibbs distribution for the position and the Maxwell-Boltzmann distribution for the velocity to give a diagonal covariance matrix

[TABLE]

It is easily verified that the stationarity condition, Eq. (6), is satisfied by the above matrices. The condition of micro-reversibility translates, here, into Onsager-Casimir symmetry, $\boldsymbol{\lambda}^{\mathrm{ir}}\boldsymbol{c}=\boldsymbol{\sigma}\boldsymbol{\sigma}^{\text{tr}}$, where $\lambda_{ij}^{\mathrm{ir}}=(\lambda_{ij}+\epsilon_{i}\epsilon_{j}\lambda_{ij})$ is the irreversible part of the drift coefficient for variables that are, respectively, even or odd under time reversal. Then, the only non-zero entry of $\boldsymbol{\lambda}^{\mathrm{ir}}$ is $\lambda_{22}^{\mathrm{ir}}=\tau^{-1}$, and it is trivial to verify the Onsager-Casimir symmetry.

*Path sampling: *To sample paths of the Brownian harmonic oscillator exactly, it is necessary to obtain the exponential of the regression matrix $\boldsymbol{\lambda}$ and the square-root of the variance matrix $\boldsymbol{\Sigma}$. From the Cayley-Hamilton theorem, the former is easily found to be

[TABLE]

while the latter is obtained from a Cholesky factorization of $\boldsymbol{\Sigma}$ into a lower triangular matrix and its transpose,

[TABLE]

See appendix C for the derivation of the coefficients.

*Parameter estimation: *We now take the time series \boldsymbol{X}=\big{\{}(x_{1},v_{1})^{\text{tr}},\ldots,(x_{N},v_{N})^{\text{tr}}\big{\}} of discrete observations of positions and velocities obtained from the exact path sampling and construct from it the sufficient statistics

[TABLE]

From these, we compute the MAP estimate for the regression matrix,

[TABLE]

and the covariance matrix,

[TABLE]

Using the above, $\boldsymbol{\lambda}^{*}$ yields the MAP estimates for the natural frequency and the relaxation time scale

[TABLE]

while $\boldsymbol{c}^{*}$ yields the MAP estimate of the spring constant and the mass, in units of $k_{B}T$

[TABLE]

The friction constant is estimated by eliminating the mass between two of the previous ratios,

[TABLE]

The three preceding equations provide the Bayes I map estimates of the oscillator parameters.

Next we use Eq.(39) for the covariance matrix in the bivariate analog of Eq.(28) for the logarithm of the posterior probability. The MAP estimates and the error bars for mass $m$ and spring constant $k$ are then

[TABLE]

These correspond to Bayes II estimates of the oscillator parameters.

In Table 1, we provide the MAP estimates and corresponding error bars for three different time-series data obtained from the simulation of the underdamped Brownian harmonic oscillator. This clearly shows that the MAP estimates of both Bayesian methods and values of the simulation parameters are in excellent agreement. The error bars have been calculated from the Hessian matrix. The Bayesian standard error in estimating the parameters is less than $1\%$ for each case.

In Fig.(1) we show the posterior distribution of the parameters as corner plots foreman2016corner , with panel (a) corresponding to Bayes I and panel (b) to Bayes II. The distributions have been shifted to the MAP estimates, which, therefore are always at the origin marked by the black dot, and scaled by the variances. The Bayesian credible intervals corresponding to $70\%$, $95\%$, and $99\%$ of the posterior probability are shown in shades of gray.

*Correlation functions and spectral densities: *The MAP estimates of the parameters provide a novel way of estimating the correlation functions and spectral densities. Their expressions can be obtained from Eq.(7) and Eq.(8) using the explicit forms of the covariance and the regression matrices. The autocorrelations are

[TABLE]

and the spectral densities are

[TABLE]

These expressions are evaluated at the Bayesian MAP estimates for the parameters. In Fig.(2), the result is compared with the autocorrelation function and the discrete Fourier transform of the time series. There is excellent agreement between the two autocorrelations. We emphasize that no numerical function fitting is required to obtain this estimate. The spectral density is, arguably, even more impressive as it interpolates through the noisy discrete Fourier transform in a sensible manner. Our estimation of the spectral density shares a methodological similarity with Burg’s classic maximum-entropy method burg1967maximum . The principal differences are (i) that we choose the Ornstein-Uhlenbeck process as the data generating model, rather than the discrete autoregressive process assumed by Burg and (ii) that we work directly on the space of trajectories rather than in the space of correlations. From a Bayesian perspective, our estimation procedure encodes the prior information that the data are generated by an Ornstein-Uhlenbeck process and, therefore, will outperform all methods berg2004power ; tassieri2012microrheology that do not incorporate this prior information, such as the one in Bera:16 .

*Model comparison: *To observe the inertial motion of a colloidal particle in an optical trap, it is necessary to observe the motion at time intervals much smaller than the momentum relaxation time, that is, to ensure $\Delta t\ll\tau$. In the opposite limit, of $\Delta t\gg\tau$, the inertia of the colloid is no longer relevant, and only purely diffusive motion can be observed. These correspond, respectively, to the Kramers and Smoluchowski limits of the Brownian harmonic oscillator. In experiment, it is often not a priori clear if the observational time scale places the system in the Kramers, Smoluchowski, or crossover regime. A Bayesian model comparison provides a principled way to answer this question, as we show below.

We consider two Ornstein-Uhlenbeck data generating models

[TABLE]

corresponding, respectively, to motion on the Kramers** (underdamped) and Smoluchowski (overd**amped) time scales. The latter is the formal $m\rightarrow 0$ limit of Eq.(34) and can be obtained systematically by adiabatically eliminating the velocity as a fast variable gardiner1984adiabatic . The overdamped oscillator is an univariate Ornstein-Uhlenbeck model with two parameters, $\gamma$ and $k$, and Bayesian parameter estimation for it was presented in bera2017fast . The two-parameter overdamped model is simpler than the three-parameter underdamped model. The question we try to answer is which of these models provides the best explanation of the data for the least number of parameters.

To do so, we compare the posterior probability of the model, given the data, by approximating the model evidence in terms of the maximum likelihood and the Ockham factors, assuming equal prior probabilities for the models as in section III. The evidence for both models is straightforwardly computed from the general expression in the Appendix B. For the overdamped model the sufficient statistics are the scalars

[TABLE]

as was first pointed out in bera2017fast .

In Fig.(3), we plot the logarithm of the Bayesian evidence $\ln P(\boldsymbol{X}|\mathcal{M}_{\alpha})$ in natural log units (nits) for overdamped and underdamped Brownian motion. The panel on the left contains the plot of the logarithm of the evidence for underdamped data while the right panel is for overdamped data. The comparison of $\ln P(\boldsymbol{X}|\mathcal{M}_{\alpha})$ for the overdamped and underdamped models clearly shows that the evidence is greater in each case for the true model of the data.

VI Summary

In summary, we have presented two Bayesian methods for inferring the parameters of a multivariate Ornstein-Uhlenbeck process, given discrete observations of the sample paths. An exact path sampling procedure has been presented and utilized to validate the Bayesian methods for the Brownian harmonic oscillator. The problem of Bayesian model comparison has been addressed and applied to select between Kramers and Smoluchowski limits of the Brownian harmonic oscillator. Future work will address the problem of parameter inference and model selection when either or both of the Markov and Gaussian properties of the process are relaxed.

Acknowledgements.

We thank the two anonymous referees for their constructive criticism, which led to an improvement in the presentation of our results. RA thanks the Isaac Newton Trust for an Early Career Support grant. RS is funded by a Royal Society-SERB Newton International Fellowship.

Appendix A MAP estimates

In this appendix, we make explicit the derivation of Eq.(24) from Eq.(23). To this end, we consider the quadratic form $\boldsymbol{\Delta}_{n}^{\text{tr}}\bm{\Sigma^{-1}\Delta_{n}}=\left(\bm{x_{n+1}}-\bm{\Lambda}\bm{x_{n}}\right)^{\text{tr}}\bm{\Sigma^{-1}}\left(\bm{x_{n+1}}-\bm{\Lambda}\bm{x_{n}}\right)$ that occurs in the expression of the logarithm of the posterior probability in Eq.(23). Expanding terms and completing the summation we have the following identity

[TABLE]

Here, we have used the definitions of sufficient statistics in Eq.(25). With this identification, the expression of the logarithm of the posterior probability in terms of the sufficient statistics is given in Eq.(24). Inspecting the expression, we recognize that the posterior is normal in $\boldsymbol{\Lambda}$ which directly yields its MAP estimate $\boldsymbol{\Lambda}^{*}=\bm{T}_{2}\bm{T}_{3}^{-1}$. Taking the derivative of the above equation with respect to $\boldsymbol{\Sigma}^{-1}$ and using the matrix identity $\partial\ln\left(\det\bm{A}\right)/\partial\bm{A}=\left(\bm{A}^{\text{tr}}\right)^{-1}$, we obtain

[TABLE]

At the maximum, we obtain the MAP estimate as $\bm{\Sigma^{\ast}}=\dfrac{1}{N}\left(\bm{T}_{1}-\boldsymbol{T}_{2}\boldsymbol{T}_{3}^{-1}\boldsymbol{T}_{2}^{\text{tr}}\right).$ The estimates of the parameters of the underdamped oscillator follow from $\bm{\Lambda^{\ast}}$ and $\bm{\Sigma^{\ast}}$ after algebraic manipulations, see Appendix of rasmussen2006gaussian for further details.

Appendix B Standard errors and evidence

We now obtain an explicit expression for the standard errors and evidence from the derivatives of the logarithm of the posterior distribution. The second partial derivatives of the logarithm of the posterior distribution with respect to $\boldsymbol{\Sigma}^{-1}$ and $\boldsymbol{\Lambda}$, at the maximum, are

[TABLE]

The mixed partial derivative at the maximum is

[TABLE]

The Hessian matrix, $\boldsymbol{A}=-\nabla\nabla\ln P\left(\bm{\theta}|\bm{X}\right)$, at the maximum is then

[TABLE]

The above expression of the Hessian matrix has been used to obtain the standard errors of the MAP estimates and the evidence of a given model.

The expression of the Bayesian evidence of a model $\mathcal{M}_{\alpha}$ with several parameters is given in terms of the likelihood and the Hessian matrix as mackay2003information

[TABLE]

This expression has been used to obtain the logarithm of the evidence for a $M$-dimensional multivariate Ornstein-Uhlenbeck model $\mathcal{M}_{\alpha}$. Using explicit forms of the likelihood and the Hessian matrix, at the maximum, the expression of the logarithm of the evidence for a model $\mathcal{M}_{\alpha}$ is given as

[TABLE]

Here $\boldsymbol{A}_{ij}$ are the elements of the Hessian matrix.

Appendix C Path sampling

In this appendix, we provide explicit expressions for $e^{-\boldsymbol{\lambda}\Delta t}$ and $\sqrt{\boldsymbol{\Sigma}}$, needed for the path sampling of a bivariate Ornstein-Uhlenbeck process in section V. From the Cayley-Hamilton theorem, $e^{-\boldsymbol{\lambda}\Delta t}$ is obtained in terms of the identity matrix and $\boldsymbol{\lambda}$ as shown in Eq.(40). The coefficients are

[TABLE]

Here $\omega=\sqrt{\omega_{0}^{2}-1/(2\tau)^{2}}$ is the frequency of the damped oscillator. $\sqrt{\boldsymbol{\Sigma}}$ is obtained from its Cholesky factorization into a lower triangular matrix and its transpose, as shown in Eq.(41). The elements of the Cholesky factor are

[TABLE]

We use the above to elucidate the results in Eq.(33) in order to obtain exactly sampled trajectories of the Brownian harmonic oscillator. In Fig.(4) we show a typical sample path of the positions and velocities and their corresponding histograms. The histograms of the positions and the velocities clearly show that the distributions are normal.