<figure><img src="image_1.png" /></figure>

HAL
open science

# Fast analytical method to integrate multivariate Gaussians over hidden variables

François C Simon, Justin Cardona, Sean Molesky, Lucien E Weiss

<figure><img src="image_2.png" /></figure>

## To cite this version:

François C Simon, Justin Cardona, Sean Molesky, Lucien E Weiss. Fast analytical method to integrate multivariate Gaussians over hidden variables. 2025. hal-04692487v2

HAL Id: hal-04692487

https://hal.science/hal-04692487v2

Preprint submitted on 24 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Fast analytical method to integrate multivariate Gaussians over hidden variables

Francois Simon¹*, Justin Cardona¹, Sean Molesky¹, Lucien E Weiss¹*

¹Department of Engineering Physics, Polytechnique Montréal, Montréal, Québec, Canada;
*Corresponding authors

## Abstract

Many biological and physical datasets can be modeled as Gaussian processes, whose parameters can be estimated using maximum likelihood estimation. The likelihood can be expressed as a product of conditional probabilities, where each term provides information about the model's behavior and dynamics. Here, we present a versatile framework for retrieving the physical properties of these models using products of univariate Gaussians that depend on a set of observed and hidden variables. To do so, our model uses a recurrence formula to integrate over the hidden variables, enabling fast and accurate computing.

This framework is particularly effective for characterizing recurrent Gaussian processes. In the last section, we show the use of our framework via three microscopy applications: determining the type of motion of a particle, modeling time-dependent changes in FRET efficiency, and detecting confinement between two particles.

## Introduction

Gaussian distributions are frequently encountered in biology and physics, for instance, when measuring a multifactorial trait in a population such as human height, estimating measurement errors, or thermal diffusion (1–3). The ubiquity of Gaussians is a direct consequence of the central limit theorem, which states that a sum of $n$ independent random variables tends towards a Gaussian distribution when $n$ is large, regardless of the distribution characteristics of the individual variables. In many complex applications, the data can be efficiently modeled by a Gaussian processes, including image generation (4), weather forecasting (5, 6), resource management (7, 8), and risk estimation (9).

Fitting a model to a dataset consists of adjusting the model parameters to maximize a goodness-of-fit metric. Bayesian inference and maximum likelihood estimation are powerful tools for performing the fitting step as they provide the most precise unbiased parameters (10). They both require computing the likelihood of the model for a given data set $X$. In general, multivariate-Gaussian distributed data is composed of $n$ observations and $n \cdot (n+1)/2$ parameters, which becomes intractable when $n$ is large. To address this, various strategies have been developed, including: assuming independence between variables (11), using low-rank approximations like principal component analysis (8, 12), or employing structured covariance matrices (13, 14).

When analyzing temporal data, such as particle trajectories, structured covariances are especially relevant as sequential data can often be modeled as a product of Gaussian conditional probabilities that depend on a linear combination of observed and hidden variables (13–15). One way to analyze temporal data in a tractable probabilistic manner with a limited number of parameters is to consider a linear state-space model (16) where hidden variables $x_i$ are updated at each time step according to $x_{i+1} = F \cdot x_i + w_i$ with $F$ called the state transition matrix and $w_i$ called the process noise and where the observed variables $y_i = H \cdot x_i + v_i$, with $H$ the observation matrix and $v_i$ the measurement noise. Such models are ubiquitous in physics as they
can model Newtons laws and stochastic processes (Langevin equations). The most widespread method to study state-space models is Kalman filtering which iteratively describes the probability distribution of the hidden states based on a factor called the Kalman Gain (16, 17). This recurrent technique is accurate and fast, with computation times that scale linearly with the number of time steps. However, the matrix multiplications and inversions required for Kalman filtering induce computation times which scale with the cube of the number of variables hindering its scalability to models with high numbers of observed and hidden variables.

An alternative way to analytically compute the likelihood of a recurrent Gaussian process is to express it as a product of Gaussian conditional probabilities, which depend on observed and hidden variables (11, 15, 18-21). For example, consider a diffusive particle that has consecutive positions $R = (r_1, r_2, ..., r_n)^t$ that are not directly visible to the observer. During tracking, the observed positions $C = (c_1, c_2, ..., c_n)^t$ are the real positions obscured by a Gaussian noise arising from photon statistics and camera noise. The joint probability of $C$ and $R$ is a multivariate Gaussian that can be expressed as a product of conditional probabilities

$$ 
\begin{aligned}
P(C, R) &= P(r_1) P(c_1|r_1) P(r_2|r_1) P(c_2|r_2) \dots P(r_n|r_{n-1}) P(c_n|r_n) \\
&= f_0(r_1) \left( \prod_{i=1}^{n-1} f_e(r_i - c_i) f_d(r_{i+1} - r_i) \right) f_e(r_n - c_n)
\end{aligned} 
\quad (1) $$

where $f_e$ and $f_d$ are two Gaussian functions with mean 0 and standard deviations $e$ the localization error and $d$ the diffusion length respectively, and $f_0$ is the Gaussian function corresponding to $P(r_1)$.

From this joint probability density function, the probability of an observed track $P(C)$ can be retrieved by integrating over the possible real positions. While this integration step is often done using sampling techniques (22), it can also be performed analytically in a far more precise and time-efficient manner using either a multivariate Gaussian (13, 23) or a recurrence formula such as in Relich et al., Simon et al., Simon et al.

Here, we generalize the recurrent formula used in (15) to build a framework that uses conditional probabilities and latent variables to efficiently model recurrent Gaussian processes. This framework, which we call a Conditional Gaussian Process (CGP) model, allows users to build custom physics-informed models for complex datasets. Considering conditional probabilities that define the relationships between observed and hidden variables enables sparse and physically relevant models that are easier to interpret with lower computation times.

We show how to deploy our method and demonstrate its capability to fit and generate new data, extract meaningful model parameters, and build statistical tests for detecting significant differences in behaviors. Specifically, we apply our approach to three examples from single-molecule microscopy: tracked particles undergoing anomalous diffusion, FRET-efficiency measurements, and interactions between multiple mobile species such as two proteins or a protein and a vesicle.

# Main

## 1 Using hidden variables to model Gaussian processes

### 1 General expression of the probability of a multivariate Gaussian with hidden states

When studying a random process, it is often difficult to determine the relationships between all the observed variables. In many cases, introducing hidden variables enables the expression of the problem as a product of conditional probabilities where each conditional probability only depends on a linear combination of a small number of observed and hidden variables as exemplified by the above noisy diffusion model (13, 21).

Consider a Gaussian process that can be described as the product of conditional probabilities that depend on linear combinations of observed data $X = (x_1, x_2, ..., x_n)^t$, and hidden variables $Y = (y_1, y_2, ..., y_m)^t$. The likelihood of the observed variables can be retrieved by integrating the joint probability of $X$ and $Y$ over
the possible hidden variables Y as shown in the following formula:

$$L(\theta | X) = P(X | \theta) = \int_{Y} \prod_{i=1}^{n+m} f_i \left( \sum_{j=1}^{n} a_{i,j} x_j + \sum_{j=1}^{m} a_{i,n+j} y_j \right) dY \quad (2)$$

with Gaussian functions $f_i$ of mean $\mu_i$ and standard deviation $\sigma_i$. Here, the joint probability density function is expressed as a product of univariate Gaussians. The joint probability density function can also be expressed as a multivariate Gaussian by performing a change of variables. With A the matrix with elements $a_{i,j}$ and D the diagonal matrix of diagonal elements ($\sigma_1^2, \sigma_2^2, \dots, \sigma_{n+m}^2$), the joint probability is a multivariate Gaussian of covariance matrix $\Sigma = A^{-1}D(A^{-1})^t$ and bias vector $\mu = A^{-1}(\mu_1, \mu_2, \dots, \mu_{n+m})^t$.

Expressing the joint probability as a product of univariate Gaussians is particularly useful when A is sparse, as it encodes a dense structured covariance matrix. However, this expression is not useless when A possesses more non-null coefficients than the number of parameters of $\Sigma_X$.

From $A^{-1}$ and D, we can compute $\Sigma$, which can be separated by observed, X, and hidden variables, Y, to identify the covariance matrix $\Sigma_X$ and the mean of the observed data $\mu_X$.

$$\Sigma = \begin{pmatrix} \Sigma_X & \Sigma_{XY} \\ \Sigma_{YX} & \Sigma_Y \end{pmatrix}, \quad \mu = \begin{pmatrix} \mu_X \\ \mu_Y \end{pmatrix}, \qquad (3)$$

Once $\Sigma_X$ is known, we must be inverted to compute the likelihood. The computation time of this inversion scales with the cube of the number of observations in the sequence.

Using conditional probabilities, we can also express the probability density of the hidden variables given the observed variables $P(Y|X = x)$ as a multivariate Gaussian with covariance matrix $\Sigma_Y - \Sigma_{YX}\Sigma_X^{-1}\Sigma_{XY}$ and mean $\mu_Y + \Sigma_{YX}\Sigma_X^{-1}(x - \mu_X)$ (24). This conditional probability density is particularly useful for inferring the hidden states. In the context of tracking, this enables denoising the track positions to estimate the velocity vector or the confinement area.

## 2 Analyzing sequential data with hidden states

Recurrent models are efficient in many applications as they use a fixed number of parameters to study sequences of any length. We can modify eq. 1 to describe a recurrent probabilistic model, as shown in the tracking example eq. . In a recurrent framework, at each time step $k$, we consider $n$ observed variables, $m$ hidden variables, and $n+m$ Gaussian conditional probabilities. These Gaussians specify the relationships between the current observed variables $x_k = (x_{k,1}, x_{k,2}, \dots, x_{k,n})$ and between the current and next hidden variables, $y_k = (y_{k,1}, y_{k,2}, \dots, y_{k,m})$ and $y_{k+1} = (y_{k+1,1}, y_{k+1,2}, \dots, y_{k+1,m})$, respectively. These relations can be represented by a bipartite graph, Fig 1.

For a recurrent process of $N$ steps, Eq. 1 can be written as:

$$L(\theta | X) = \int_{Y} I(y_1) \prod_{k=1}^{N} \prod_{i=1}^{n+m} f_i \left( \sum_{j=1}^{n} a_{i,j} x_{k,j} + \sum_{j=1}^{m} a_{i,n+j} y_{k,j} + \sum_{j=1}^{m} a_{i,n+m+j} y_{k+1,j} \right) dY \quad (4)$$

where $I(y_1)$ is an initialization term corresponding to a product of Gaussians with the initial hidden variables, specifically, $I(y_1) = \prod_{i=1}^{m} g_i (\sum_{j=1}^{m} b_{i,j} y_{1,j})$.

In this process, the covariance matrix $\Sigma_X$ has $N(N+1)/2$ parameters that are structured based on a model with a maximum of $(n+m) \cdot (n+2m) + m^2$ parameters. Using the Gaussian elimination principle, the number of useful parameters is further reduced to $(n+m)(n+m+1)/2 + 2m^2$. In our recurrence framework, A is banded with a maximum bandwidth of $n+2m$. Inverting A then scales linearly with the number of time steps $N$ with a maximum complexity of $N(n+m)(n+2m)^2$.

## 3 Computing the likelihood with a recurrence formula

As seen in the previous section, computing the likelihood of a model with the multivariate Gaussian approach is computationally demanding. As in the general case, it requires a matrix inversion that scales as
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Illustrated example of Recurrence Gaussian network with observed and hidden variables. Each Gaussian depends on a few variables specified by the branches of the graph. Colors represent the same variable or Gaussian at different time points.</figcaption></figure>

$n^3$, where $n$ is the number of observed variables which quickly becomes prohibitively expensive for long sequences.

An alternative to the multivariate Gaussian approach is to iteratively compute the probability of the data from the product of the univariate Gaussians by integrating over the hidden variables. Based on the property that the product of two Gaussians, $f(x)$ and $g(x)$, with means $\mu_f$, $\mu_g$ and standard deviations $\sigma_f$, $\sigma_g$ can be rearranged into a Gaussian $\eta(x)$ times a scaling factor $\phi$ which is also a Gaussian of $\mu_f$ and $\mu_g$: $f(x) \cdot g(x) = \phi \cdot \eta(x)$ where $\phi$ and $\eta$ are the following Gaussian distributions:

$$ \phi = \frac{\exp\left(-\frac{(\mu_f - \mu_g)^2}{2(\sigma_f^2 + \sigma_g^2)}\right)}{\sqrt{2\pi}(\sigma_f^2 + \sigma_g^2)} \quad \text{and} \quad \eta(x) \text{ a Gaussian of mean } \mu_{fg} = \frac{\sigma_f^2\mu_g + \sigma_g^2\mu_f}{\sigma_f^2 + \sigma_g^2} \quad \text{and variance } \mu_{fg} = \frac{\sigma_f^2\sigma_g^2}{\sigma_f^2 + \sigma_g^2}. \quad (5) $$

Applied to the integration problem, this enables the integration to be computed analytically. For a given recurrent multivariate Gaussian, the most efficient conditional probability expression possesses the fewest dependencies on the hidden states.

When all the Gaussians are fully connected to the current and next hidden variables in the graph from Fig 1, the Gaussian-elimination principle can be applied to reduce the number of connections:

$$ L(\theta | X) = \int_Y I'(y_1) \prod_{k=1}^{N} \prod_{i=1}^{n+m} f_i \left( \sum_{j=1}^{n} a_{i,j} x_{k,j} + \sum_{j=i}^{m} a_{i,n+j} y_{k,j} + \sum_{j=\max(1, i-m)}^{m} a_{i,n+m+j} y_{k+1,j} \right) dY \quad (6) $$
with $I'(y_1) = \prod_{i=1}^{m} g(\sum_{j=i}^{m} a_{0,j} y_{1,j})$.

At a given time step $k$, there are a maximum of $n+2m$ Gaussians that depend on a given hidden variable $y_{k,j}$, where each Gaussian has a maximum of $2m+1$ terms (after calculating the contribution of the observed data, $\sum_{j=1}^{n} a_{i,j} x_{k,j}$ which can be included in the bias). We then have to apply equation $3n+2m-1$ times. Integration over $y_{k,1}$ results in the removal of the last Gaussian that depends on $y_{k,1}$. This leaves $n+m+1$ Gaussians that depend on $y_{k,2}$. Repeating this process for all hidden variables results in a maximal time complexity of $Nm(n+2m-1)(2m+1)$. However, the time complexity can be strongly reduced if we limit the number of hidden variables affecting each Gaussian. For example, in the method applications described in Section 2, our Gaussians depend on a maximum of 3 hidden variables. Limiting the number of hidden variables affecting each Gaussian reduces the computational time to be proportional to the square of the number of hidden variables instead of the cube.

## 2 Application

Our method can be applied to any random phenomena that can be approximated by Gaussian processes. In this section, we show several examples of applications in the field of single-particle measurements.
# 1 Determining and characterizing the type of motion

To create a model that characterizes the type of motion of single-particle tracks, we must account for diffusion, localization error, and either directed motion or confinement (15, 19, 20, 25).

By assuming that the motion is independent in the x and y axes, we can provide a good description of an anomalous track with the following recurrent Gaussian model: At each new time step i, we consider 3 new univariate Gaussians that depend on the observed position $o_i$ and on two hidden variables: the real particle position r (without localization error) and the anomalous value u. The first Gaussian represents the localization error, the second Gaussian represents the evolution of the particle's real positions, and the third Gaussian represents the evolution of the anomalous value.

$$ \begin{array}{c|cccccc} & o_i & r_i & u_i & r_{i+1} & u_{i+1} & \text{bias} \\ \hline \text{Gaussian 1} & a_{11} & a_{12} & 0 & 0 & 0 & b_1 \\ \text{Gaussian 2} & 0 & a_{22} & a_{23} & a_{24} & 0 & b_2 \\ \text{Gaussian 3} & 0 & 0 & a_{33} & 0 & a_{35} & b_3 \end{array} \qquad (7) $$

A major advantage of this model is that it can fit both directed motion and confined motion depending on its parameter values (15). To verify the capability of our model to fit Brownian, directed, and confined motion behaviors, we compared a simulated set of tracks using a predefined model to tracks generated with our fitted model trained on the original dataset (Fig 2). For Brownian, directed, and confined tracks, we found that the CGP model generates qualitatively similar tracks both visually and in terms of Mean Squared Displacements (MSD) (Fig 2). Notably, there is a discrepancy in the velocity of individually directed tracks with fixed velocities and random orientations. This is because Gaussian processes cannot account for a ring-shaped probability density function and thus cannot perfectly match a dataset with fixed velocities with variable orientations. Nevertheless, the ensemble MSDs of the original tracks and the tracks generated by the CGP model are very similar (Fig 2c).

To show that our method can also be used to determine the model parameters, we varied the diffusion length and the confinement factor or the velocity. This produced accurate estimations for a wide range of biologically relevant parameters (Fig 2, right).

# 2 Modeling FRET dynamics

FRET (Förster Resonance Energy Transfer) microscopy enables the study of nanoscale interactions between or within molecules (26–28). In FRET, a donor fluorophore is excited, and a photon is either emitted by the same fluorophore, or the energy is transferred to a nearby acceptor fluorophore. This process is highly distance-dependent (Eq. 8) and thus can provide detailed spatial information retrieved using the FRET efficiency and calibrations (28, 29).

In FRET, the efficiency follows:

$$ e = \frac{F_A}{F_D + F_A} = \frac{1}{1 + (R/R_0)^6} \qquad (8) $$

where e is the FRET efficiency, $F_D$ the donor fluorescence, $F_A$ the acceptor fluorescence, R the donor/acceptor distance, and $R_0$ the Förster distance.

When an observed FRET efficiency measurement follows a Gaussian-like distribution (28, 30), this may arise from either random fluctuations in photon emission or physical changes in the efficiency, e.g. distance or orientation changes between the donor and acceptor. While efficiency changes are often interpreted as a change of state, there are cases where it makes more sense to consider distance changes occurring within a single state, such as in intrinsically disordered proteins (31, 32). In these cases, we propose to model the evolution of the FRET efficiency as a diffusive process confined in a fixed potential well characterized by a confinement factor.

To quantify the dynamics of FRET efficiency changes within states using a CGP model, we need to approximate the relationship between the signal and the FRET efficiency as linear. A simple model considers 3 variables at each time point i: the observed donor and acceptor signals ($F_i^d$ and $F_i^a$) and the hidden FRET efficiency $e_i$. These variables are linked by 3 Gaussians: The first Gaussian represents the diffusion
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Results of the fitting of different types of motions by our CGP model (diffusive, confined and directed). The models were trained on data sets composed of 1,000 tracks with 40 time steps, 0.02 µm localization error, and frame lengths of 0.02 s. Left: original tracks (simulations) and tracks generated by the model presented in Equ. 7 after fitting to the data sets. Middle: Comparison of the ensemble MSD curves for the original tracks with the tracks generated with our model. Right: parameter estimate when varying the diffusion length, the velocity or the confinement factor. **a:** Brownian tracks with a diffusion length per step $d = \sqrt{2D\Delta t} = 0.1\mu m$. **b:** Confined diffusive tracks with a diffusion length $d = 0.1\mu m$ and a confinement factor per step of 0.3. **c:** Directed tracks with a constant velocity of 0.02 µm per step and Brownian changes of orientations with an angular diffusion coefficient of 0.5 Rad²/s (diffusion length $d = 0\mu m$).</figcaption></figure>

of the FRET efficiency $e_i$ in a potential well towards its next value $e_{i+1}$. The second and the third Gaussians represent the probability of the donor and acceptor signals, linking these signals and the hidden FRET efficiency $e_i$. The joint probability density function is then the product of conditional probabilities $P(e_1) \prod_{i=1}^{N} P(e_{i+1}|e_i) P(F_i^D|n_i, e_i) P(F_i^A|e_i)$.

$$ \begin{array}{c} \begin{array}{c|ccccc} & F_i^D & F_i^A & e_i & e_{i+1} & \text{bias} \\ \hline \text{Gaussian 1} & 0 & 0 & a_{13} & a_{14} & b_1 \\ \text{Gaussian 2} & a_{21} & 0 & a_{23} & 0 & b_2 \\ \text{Gaussian 3} & 0 & a_{32} & a_{33} & 0 & b_3 \end{array} \end{array} \tag{9} $$

To verify the capacity of our model to estimate relevant parameters, we simulated FRET signal traces of particles for which the distance between the donor and acceptor follows Brownian motion confined in a potential well around an equilibrium distance (Fig. 3a). Such a potential well can be characterized by a confinement factor $I$ representing the movement towards the equilibrium distance (15) (Fig. 3b). Modifying the value of the confinement factor affects the FRET donor and acceptor signals as illustrated in Fig. 3c as long as the equilibrium value is not too far from the Förster distance (linear regime around the Förster distance). Of course, if the equilibrium distance between the donor and acceptor is outside the dynamic range
around the Förster radius 3-8 nm, the confinement factor cannot be estimated as the FRET efficiency will be constant. To evaluate our model's ability to predict the confinement factor, we simulated noisy FRET donor and acceptor signals of particles with a range of confinement factors *l* (See Methods for simulation details). We find that our model accurately estimates the confinement factor for a wide range of values (Fig. 3d). Next, we tested our FRET model on experimental data using the publicly available FRET data from *Asher et al.*, 2021 (Fig. 3e). By comparing the likelihood of our model to a model without correlations between donor and acceptor signals, our model detects significant changes in FRET efficiency (p-value = 0.0079, see Methods and Fig. S1a for more details). As a negative control, we applied the same test on the FRET signals obtained from a fixed condition where we do not expect movement, and therefore, the observed changes in FRET efficiency should only be due to stochastic noise. As expected, we rejected the alternative hypothesis of donor/acceptor distance changes (p-value = 0.36, Fig. S1b).

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

**Figure 3.** **a:** Potential well of a donor/acceptor distance centered around an equilibrium distance of 6.5 nm. **b:** Stepwise displacement resulting from the spring potential energy with *h* the equilibrium distance and *x*ᵢ the current position. **c:** Examples of simulated FRET signals with different confinement factors. **d:** Confinement factor predicted by our model as a function of the true confinement factor of our simulated tracks (1.000 tracks of 100 time steps). **e:** FRET signal from (30) that we tested with our statistical tool.

### 3 Unraveling multi-particle interactions

Modeling the interactions between multiple mobile species, for example, proteins diffusing on mobile vesicles, poses a challenge due to the many possible observable behaviors; however, our framework can be easily adapted to characterize a wide range of multibody problems. To verify the capability of our model to consider complex interactions between multiple bodies, we simulated tracks of a protein and a vesicle where the protein is trapped. The vesicle moves in a slow, directed motion, and the protein shows confined diffusion around the vesicle center. Using this model, we tested our capacity to detect correlated motion and quantify the parameters that characterize the interaction.

We can model such a process as a product of five Gaussians that depend on 2 observed variables (*o*¹, *o*²) and 3 hidden variables (*r*¹, *r*², *u*) with the vesicle positions *r*¹, estimated by *o*¹, the confined particle positions *r*², estimated by *o*² and the drift vector of the vesicle *u*:

<table><thead><tr><th></th><th>o<sub>i</sub><sup>1</sup></th><th>o<sub>i</sub><sup>2</sup></th><th>r<sub>i</sub><sup>1</sup></th><th>r<sub>i</sub><sup>2</sup></th><th>u<sub>i</sub><sup>2</sup></th><th>r<sub>i+1</sub><sup>1</sup></th><th>r<sub>i+1</sub><sup>2</sup></th><th>u<sub>i+1</sub><sup>2</sup></th><th>bias</th></tr></thead><tbody><tr><td>Gaussian 1</td><td>0</td><td>0</td><td>0</td><td>0</td><td>a<sub>15</sub></td><td>0</td><td>0</td><td>a<sub>18</sub></td><td>b<sub>1</sub></td></tr><tr><td>Gaussian 2</td><td>0</td><td>0</td><td>a<sub>23</sub></td><td>0</td><td>a<sub>25</sub></td><td>a<sub>26</sub></td><td>0</td><td>0</td><td>b<sub>2</sub></td></tr><tr><td>Gaussian 3</td><td>0</td><td>0</td><td>a<sub>33</sub></td><td>a<sub>34</sub></td><td>0</td><td>0</td><td>a<sub>27</sub></td><td>0</td><td>b<sub>3</sub></td></tr><tr><td>Gaussian 4</td><td>a<sub>41</sub></td><td>0</td><td>a<sub>43</sub></td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>b<sub>4</sub></td></tr><tr><td>Gaussian 5</td><td>0</td><td>a<sub>52</sub></td><td>0</td><td>a<sub>54</sub></td><td>0</td><td>0</td><td>0</td><td>0</td><td>b<sub>5</sub></td></tr></tbody></table>
Here, the first Gaussian represents the evolution of the drift vector, the second represents the evolution of the vesicle position, the third represents the evolution of the confined particle position, and the final two Gaussians represent the localization errors of the estimated positions.

We first simulated tracks with known parameters (Fig 4a, left) and fit our model (Eq. 10) to these tracks. The first important question that our model can address is to test whether or not we can statistically show two tracks to be confined to each other. To do so, we can consider a null hypothesis without confinement between the two particles ($a_{34} = 0$), which is a specific case of Eq. 10, that represents the alternative hypothesis $H_1$. Then, we compared the distribution of the likelihood ratio statistics $\Lambda = -2(\log(L_0) - \log(L_1))$ of the original tracks and the tracks generated by our model under the null hypothesis (Fig 4b). There, we can see that 97.9 % of the original tracks have a $\Lambda > 0$ while only 4.2 % of the tracks generated from the $H_0$ hypothesis have a $\Lambda$ above 0. This proves that the test has both high sensitivity and high specificity. We then varied the confinement factor to verify the capability of our model to estimate parameters. Fig. 4c shows very accurate predictions at low confinement factors and a small bias at higher confinement factor values.

<figure><img src="image_11.png" /><figcaption>**Figure 4.** **a:** 30 time-point tracks of the vesicle (red) and the confined particles (blue) for the original tracks that we trained our model with, for tracks generated by our model assuming the null hypothesis $H_0$ (no correlation: $a_{34} = 0$) or the alternative hypothesis $H_1$. Original tracks: the vesicle is in linear motion with fixed velocity (0.05µm/step of 0.02 s), a random orientation, and changes of orientations with a rotational diffusion coefficient of 0.1 Rad²/s and a localization error of 0.02µm. The diffusive particle is confined around the vesicle with a (continuous) confinement factor of 0.4 per second, which correspond to a stepwise confinement factor of 0.33, has a diffusion coefficient of 0.25 µm/s and a localization error of 0.02µm. The diameter of the confinement is then 0.45µm (with a 95% interval). **b:** Histograms of the likelihood ratio statistics from the original tracks and from the tracks generated by our model assuming the null hypothesis. **c:** Curve of the estimated confinement factor depending on the true confinement factor (per step) for similar tracks than in **a**.</figcaption></figure>

### 3 Discussion

The presented CGP framework uses a fast and adaptable analytical recurrence formula to compute the likelihood of sequential datasets. We showed that CGP models are particularly useful for analyzing phenomena whose probability can be expressed as a product of univariate Gaussians that depend on a linear combination of a limited number of observed and hidden variables. Using conditional probabilities improves interpretability and enables sparse expressions of multivariate Gaussians with dense covariance matrices. This sparsity enables the recurrent integration process to have a time complexity that increases linearly with the number of time points in the sequence, a major improvement compared to classical Gaussian processes whose time-complexity scales with the cube of the sequence length. Like Kalman filtering, CGP models have the advantage of requiring very little memory compared to the matrix approach (16).

While the computation time scales with the cube of the number of observed variables in Kalman filtering, CGP models offer two levers to reduce its time complexity: 1) Limiting the number of hidden variables considered by each univariate Gaussian, which naturally occurs when breaking down complex physical models. 2) Using specific patterns of relationships between hidden variables coupled with graph theory to find optimal integration orders. This enables model complexities that scale with the square or even linearly with the number of variables (see the Methods Section for an example).
Another major advantage of CGP models is that it does not require a specific order of integration. Therefore, the integration can be performed in parallel for every pair time steps ($k = 2i$), dividing the number of Gaussians by two. For a sequence of $N$ time steps, repeating this process (for $k = 4i, k = 8i, ...$) results in a runtime proportional to $\log(N)$ (33). In comparison, in classical RNNs or in Kalman filtering, computing the current hidden states requires the previous hidden states, preventing parallelization along sequences. As a result, CGP models can offer major speedup with runtimes proportional to $\log(N)$ instead of $N$.

Related to the last point, CGP models do not depend on the directionality of the sequences. This property is particularly relevant when aiming to tackle multidimensional data such as spatial or spatiotemporal maps. In these configurations, the speed-up of CGP models is even more drastic, as the runtime decreases from $N_1 \cdot N_2 \cdot ...$ to $\max(\log(N_1), \log(N_2), ...)$.

This model can also be used for inference. First, it enables to infer the probability density function of the hidden variables during the sequence. This can be particularly useful for denoising data such as seismic data (9) and improving the mapping of contaminants (8, 34). In tracking, this is also useful for refining the particle position (21, 35) and improving spatiotemporal resolution (36). Second, we can also predict the probability distribution of future (or past) observed and hidden variables. Predictions of the next states are particularly relevant in many applications such as weather forecasting (5, 6) or resource management (7).

In our framework, both interpolation and extrapolation can be performed quite easily by integrating over all the hidden variables except the hidden variable of interest. Forecasting applications require considering the next steps in the joint probability density function. This can be done by taking the exact same relations between hidden and measured variables and considering the measured variables as additional hidden variables.

In this implementation, the method cannot handle abrupt random state transitions; however, modeling such transitions is important for many processes. In tracking for instance, we often have transitions between immobile and diffusive states (11, 21, 37, 38). The main challenge when considering random state transitions is that the number of state sequences to consider increases exponentially with the sequence length. Monte Carlo sampling could provide a satisfactory solution to limit the number of sequences to consider as it scales with the number of samples instead of the power of the track length (19, 20). In ExTrack, a tool that considers diffusive and bound states, we constrained the number of sequences by merging sequences that produce similar internal continuous states at a given time point during the integration process (21). In principle, this approach can be extended to manage transitions between states with more complex behaviors made possible by the CGP model. Such an approach might provide better trade-offs between accuracy and speed than sampling methods.

We believe this architecture would be well-suited to be the basic block of a new type of physics-informed neural network that considers multiple Gaussian processes expressed from conditional probabilities. We expect the resulting network to be more reliable than networks that are not probabilistic and faster than those that rely on slow and inaccurate sampling (34) or variational Bayesian inference (6, 39). In the field of single-particle tracking, machine learning has proved useful, but surprisingly, non-probabilistic approaches are preferred to probabilistic approaches (40) even if the latter seem particularly adapted to tracking where random processes play a major role (15, 19, 20). This preference might be explained by a higher exposure of the classical networks (CNN, RNN, LSTM, transformers) that are reputed to provide easy-to-deploy solutions for any type of data without prior knowledge. However, they have several disadvantages that make them little adapted for scientific application. First, they require large datasets that can be very challenging or impossible to collect. Second, they are supervised models, meaning that they require a training step on data that are labeled with the parameters that we aim to estimate. In tracking or in other highly stochastic problems, these parameters are difficult to assess precisely from real data. The classical workaround to obtain reliable labels is to train the models with simulated data (40). However, these approaches have difficulties to generalize to unseen data (15, 41). In (15) for instance, we show that a state-of-the art RNN model designed to characterize anomalous diffusion struggles to detect simple directed motion with localization error. Another well-known disadvantage of the classical networks is that they are difficult if not impossible to interpret beyond their outputs, preventing us to gain knowledge of the data and their underlying mechanisms. The capacity of the CGP architecture to include prior knowledge makes it more difficult to deploy
than the classical machine learning approaches but it ensures better interpretability, better generalizability and it works for small datasets that classical methods cannot analyze.

In conclusion, CGP models are easy-to-deploy, scalable, and interpretable frameworks that are well-suited to study random processes that can be approximated as a product of conditional Gaussian probabilities. They can be used to characterize the properties of the studied signal trajectories, to build a statistical test to compare alternative hypotheses, and be used for inference. As a generative model, CGP models can also be used to simulate signals. Finally, they could prove useful to build more complex physics-informed neural networks.

# 4 Methods

## 1 Creating a recurrent Gaussian process

To create a valid model, one needs to respect some conditions. An important prerequisite is that the number of variables must equal the number of univariate Gaussians. This has two consequences: First, at each time step, the number of new univariate Gaussians must equal the number of new variables (observed + hidden). Second, with *N* the number of time steps, *n* the number of observed variables and *m* the number of hidden variables, the recurrence form of our model describes *N*(*n* + *m*) univariate Gaussians that depend on *N*(*n* + *m*) + *m* variables. This is why in eq. 3, we need a term *I* that initializes the distribution of the hidden variables and obtains *N*(*n* + *m*) + *m* univariate Gaussians. Another important condition is that *A* is not a singular matrix. In practice, this means that every hidden variable needs at least one coefficient that cannot be canceled. In general, as long as the relationships are inferred from a probabilistically valid conditional model all these conditions should be respected. NB: as exemplified in eq. 11, if we want to consider a time independent variable, it is important to inform the coefficients of that variables for the next step instead of the current step. Indeed, informing the coefficient of that variable for the current step will be problematic as *I* and the Gaussians from the first iteration of the recurrence process will cancel one Gaussian and result in a singular matrix.

## 2 FRET simulations

We simulated the continuous Brownian motion of the donor/acceptor distance with 100 sub-steps per sampled time step $\Delta t = 0.02$ s. At each substep, the distance *x* evolves according to Brownian motion ($D = 0.000025\mu m/s$) and is attracted towards the center of the potential well *h* by $l(h - x)/100$. At each acquisition time point (stroboscopic imaging), we considered that a donor fluorophore emits a random number of photons *n* that follow a Poisson distribution of mean 200 photons. The number of photons absorbed by the acceptor *k* then follows a binomial distribution of parameters *n* and *e* the efficiency. The signal captured by the camera is then $a(n-k)+\eta_1$ for the donor signal and $b k+\eta_2$ for the acceptor signal, with *a* and *b* factors that account for the filters (resp. 0.7 and 0.75 in our case). $\eta_1$ and $\eta_2$ are random variables that follow a Poisson distribution of mean 10 to simulate the readout noise of the camera.

$$ \begin{array}{c|ccccc} & F_i^D & F_i^A & e_i & e_{i+1} & \text{bias} \\ \text{Gaussian 1} & 0 & 0 & 0 & a_{14} & b_1 \\ \text{Gaussian 2} & a_{21} & 0 & a_{23} & 0 & b_2 \\ \text{Gaussian 3} & 0 & a_{32} & a_{33} & 0 & b_3 \end{array} \qquad (11) $$

## 3 Determining the confinement factor

We determine the confinement factor by fitting our FRET model to a dataset using the fact that the Gaussian 1 of 9 is a Gaussian of $e_{i+1} - (1 - l)e_i + lh$ and standard deviation *d*. Here *d* represents the diffusion of the efficiency and *h* the equilibrium position (bias).
## 4 Model significance by likelihood ratio test

We performed statistical tests by considering a null hypothesis $H_0$ which is a specific case of a more complex alternative hypothesis $H_1$. We then compute the likelihood ratio statistics, $\Lambda$, where $\Lambda = -2(\log(L_0) - \log(L_1))$ with $L_0$ and $L_1$ the likelihoods of the null and alternative hypotheses respectively. To obtain the p-value from $\Lambda$, we need to know the distribution of $\Lambda_0$ the statistics under the null hypothesis. To do so, we generate 100.000 FRET signals of 69 time points under the null hypothesis with its most likely parameters (bootstrapping) and compute $\Lambda_0$ distribution. We can then retrieve the p-value as, by definition, it equals $P(\Lambda_0 > \Lambda_{obs})$.

### Track simulations

While measurements are made at discrete time intervals, samples evolve continuously. To simulate approximately continuous tracks, each measurement time step was subdivided into 100 substeps for simulating movement (15). We applied the Brownian diffusion and anomalous movements to each of these substeps. For confined motion, the latter step moves the particle toward the center of the potential well center proportionally to the distance multiplied by a scaled confinement factor, $l/20$. In simulations with a moving potential well, the center moves according to the well's diffusion coefficient, $q$. For directed motion, we applied a shift of constant velocity at each time step after to the diffusion step. The orientation of the directed motion could vary according to a rotational diffusion coefficient. The particle's position after the 20 sub-steps was set to the particle's real position, $r_i$, and the localization error was added to create the observed position $c_i$. This process is equivalent for undivided frames for directed motion, as the diffusion and the directed motion steps are independent. However, for the confined model, where the diffusion of the particle influences the anomalous step, thus there is a dependence between the diffusion and confinement steps.

### Model parameter estimations

In order to estimate model parameters, we typically need to set constraints on the values taken by the $a_{i,j}$ elements. In case of anomalous motion for instance, the values depend on 4 actual parameters: the localization error $\sigma$, the diffusion length $d$, the anomalous parameter $l$ (either the velocity or the confinement factor) and the diffusion length of the change of the anomalous parameter $q$ (15). With constraints, Eq. 7 becomes:

$$ \begin{array}{c|cccccc} & o_i & r_i & u_i & r_{i+1} & u_{i+1} & \text{bias} \\ \hline \text{Gaussian 1} & 1/\sigma & -1/\sigma & 0 & 0 & 0 & 0 \\ \text{Gaussian 2} & 0 & (1-l)/d & l/d & -1/d & 0 & 0 \\ \text{Gaussian 3} & 0 & 0 & 1/q & 0 & -1/q & 0 \end{array} \qquad (12) $$

for confined particles and:

$$ \begin{array}{c|cccccc} & o_i & r_i & u_i & r_{i+1} & u_{i+1} & \text{bias} \\ \hline \text{Gaussian 1} & 1/\sigma & -1/\sigma & 0 & 0 & 0 & 0 \\ \text{Gaussian 2} & 0 & 1/d & 1/d & 1/d & 0 & 0 \\ \text{Gaussian 3} & 0 & 0 & 1/q & 0 & -1/q & 0 \end{array} \qquad (13) $$

for directed motion (NB: in this directed motion model, the scale of the velocity is stored in the initialization Gaussians ($u_0$ follows a Gaussian distribution of mean 0 and standard deviation $l$). Eq. 10 also becomes:

$$ \begin{array}{c|cccccccc} & o_i^1 & o_i^2 & r_i^1 & r_i^2 & u_i^2 & r_{i+1}^1 & r_{i+1}^2 & u_{i+1}^2 & \text{bias} \\ \hline \text{Gaussian 1} & 0 & 0 & 0 & 0 & 1/q_2 & 0 & 0 & -1/q_2 & b_1 \\ \text{Gaussian 2} & 0 & 0 & 1/d_1 & 0 & 1/d_1 & -1/d_1 & 0 & 0 & b_2 \\ \text{Gaussian 3} & 0 & 0 & (1-l)/d_2 & l/d_2 & 0 & 0 & -1/d_2 & 0 & b_3 \\ \text{Gaussian 4} & 1/\sigma & 0 & -1/\sigma & 0 & 0 & 0 & 0 & 0 & b_4 \\ \text{Gaussian 5} & 0 & 1/\sigma & 0 & -1/\sigma & 0 & 0 & 0 & 0 & b_5 \end{array} \qquad (14) $$

For more details, see the function 'constraint_function' in the file Tracking.py of the Github page https://github.com/FrancoisSimon/CGP.
## Data Generation from CGP models

To generate data from CGP models we reconstructed the covariance matrix and the bias vector of the multivariate Gaussian expression from the product of univariate Gaussian expression fitted by our model (see section 2). We then used the python function numpy.random.multivariate_normal to generate the data.

## Example of model that scale linearly with the number of hidden variables

We can imagine many models with computation times that scale linearly with the number of hidden variables. An important prerequisite to create a relevant hidden variable is that it must be connected to the observed data. Here, we describe a toy example of such a type of model:

<table><thead><tr><th></th><th>o<sub>i</sub></th><th>r<sup>1</sup><sub>i</sub></th><th>r<sup>2</sup><sub>i</sub></th><th>r<sup>3</sup><sub>i</sub></th><th>r<sup>4</sup><sub>i</sub></th><th>r<sup>1</sup><sub>i+1</sub></th><th>r<sup>2</sup><sub>i+1</sub></th><th>r<sup>3</sup><sub>i+1</sub></th><th>r<sup>4</sup><sub>i+1</sub></th><th>bias</th></tr></thead><tbody><tr><td>Gaussian 1</td><td>a<sub>11</sub></td><td>0</td><td>0</td><td>0</td><td>a<sub>15</sub></td><td>0</td><td>0</td><td>0</td><td>0</td><td>b<sub>4</sub></td></tr><tr><td>Gaussian 2</td><td>0</td><td>a<sub>22</sub></td><td>0</td><td>0</td><td>a<sub>25</sub></td><td>a<sub>26</sub></td><td>0</td><td>0</td><td>0</td><td>b<sub>5</sub></td></tr><tr><td>Gaussian 3</td><td>0</td><td>0</td><td>a<sub>33</sub></td><td>0</td><td>a<sub>35</sub></td><td>0</td><td>a<sub>37</sub></td><td>0</td><td>0</td><td>b<sub>2</sub></td></tr><tr><td>Gaussian 4</td><td>0</td><td>0</td><td>0</td><td>a<sub>44</sub></td><td>a<sub>45</sub></td><td>0</td><td>0</td><td>a<sub>48</sub></td><td>0</td><td>b<sub>3</sub></td></tr><tr><td>Gaussian 5</td><td>0</td><td>0</td><td>0</td><td>0</td><td>a<sub>55</sub></td><td>0</td><td>0</td><td>0</td><td>a<sub>59</sub></td><td>b<sub>1</sub></td></tr></tbody></table>

(15)

In the context of tracking, the Gaussian 1 could describe the discrepancy between an observed particle position and its real position, while the other Gaussians can describe a range of processes that affect the particle true position (e.g. different types of confined, directed or oscillatory motion).

# 5 Supplemental figures
a

<figure><img src="image_12.png" /></figure>

b

<figure><img src="image_13.png" /></figure>

**Supplementary Figure 1.** Distribution of the Λ statistic of data that follow the null hypothesis called Λ₀ obtained from sampling 100.000 data sequences of 69 time steps that follow the null hypothesis with its most likely parameters. Vertical line: statistics Λ of the observed FRET signal. The p-value is the area under the curve. **a:** Distribution obtained after fitting of the models to the observed FRET signal shown in Fig. 3. **b:** Distribution obtained after fitting of the models to the observed FRET signal shown in the sub-plot. Signal retrieved from the bottom left plot of Extended Data Fig. 3 from Asher et al. (2021), time points 3.8 s to 6.12 s.
# Code and data availability

Our framework CGP for recurrent Gaussian processes is available at https://github.com/FrancoisSimon/CGP.

# Competing interests

The authors declare no competing interests.

# References

1. Shinichi Nakagawa and Holger Schielzeth. Repeatability for gaussian and non-gaussian data: a practical guide for biologists. *Biological Reviews*, 85(4):935-956, 2010.

2. Ronald Forrest Fox. Gaussian stochastic processes in physics. *Physics Reports*, 48(3):179-283, 1978.

3. Albert Einstein. Über die von der molekularkinetischen theorie der wärme geforderte bewegung von in ruhenden flüssigkeiten suspendierten teilchen. *Annalen der physik*, 4, 1905.

4. Jascha Sohl-Dickstein, Eric Weiss, Niru Maheswaranathan, and Surya Ganguli. Deep unsupervised learning using nonequilibrium thermodynamics. In *International conference on machine learning*, pages 2256–2265. PMLR, 2015.

5. Niya Chen, Zheng Qian, Ian T Nabney, and Xiaofeng Meng. Wind power forecasts using gaussian processes and numerical weather prediction. *IEEE Transactions on Power Systems*, 29(2):656–665, 2013.

6. James Donnelly, Alireza Daneshkhah, and Soroush Abolfathi. Forecasting global climate drivers using gaussian processes and convolutional autoencoders. *Engineering Applications of Artificial Intelligence*, 128:107536, 2024.

7. Aaron Zeng, Hodde Ho, and Yao Yu. Prediction of building electricity usage using gaussian process regression. *Journal of Building Engineering*, 28:101054, 2020.

8. Dmitrii Shadrin, Artjom Nikitin, Polina Tregubova, Vera Terekhova, Raghavendra Jana, Sergey Matveev, and Maria Pukalchik. An automated approach to groundwater quality monitoring—geospatial mapping based on combined application of gaussian process regression and bayesian information criterion. *Water*, 13 (4):400, 2021.

9. Ricard Durall, Ammar Ghanim, Mario Ruben Fernandez, Norman Ettrich, and Janis Keuper. Deep diffusion models for seismic processing. *Computers & Geosciences*, 177:105377, 2023.

10. Hossein Pishro-Nik. Introduction to probability, statistics, and random processes. 2016.

11. Fredrik Persson, Martin Lindén, Cecilia Unoson, and Johan Elf. Extracting intracellular diffusive states and transition rates from single-molecule tracking data. *Nature methods*, 10(3):265, 2013.

12. Neil Lawrence. Gaussian process latent variable models for visualisation of high dimensional data. *Advances in neural information processing systems*, 16, 2003.

13. Peter K Relich, Mark J Olah, Patrick J Cutler, and Keith A Lidke. Estimation of the diffusion constant from intermittent trajectories with variable position uncertainties. *Physical Review E*, 93(4):042401, 2016.

14. Jack M Wang, David J Fleet, and Aaron Hertzmann. Gaussian process dynamical models for human motion. *IEEE transactions on pattern analysis and machine intelligence*, 30(2):283–298, 2007.

15. Francois Simon, Guillaume Ramadier, Ines Fonquernie, Janka Zsok, Sergiy Patskovsky, Michel Meunier, Caroline Boudoux, Elisa Dultz, and Lucien E Weiss. Detecting directed motion and confinement in single-particle trajectories using hidden variables. *bioRxiv*, pages 2024–04, 2024.

16. James D Hamilton. State-space models. *Handbook of econometrics*, 4:3039–3080, 1994.

17. Rudolph Emil Kalman. A new approach to linear filtering and prediction problems. 1960.

18. Kaare Brandt Petersen, Michael Syskind Pedersen, et al. The matrix cookbook. *Technical University of Denmark*, 7(15):510, 2008.

19. Paddy J Slator and Nigel J Burroughs. A hidden markov model for detecting confinement in single-particle tracking trajectories. *Biophysical journal*, 115(9): 1741–1754, 2018.

20. Jason Bernstein and John Fricks. Analysis of single particle diffusion with transient binding using particle filtering. *Journal of theoretical biology*, 401:109–121, 2016.

21. Francois Simon, Jean-Yves Tinevez, and Sven van Teeffelen. Extrack characterizes transition kinetics and diffusion in noisy single-particle tracks. *Journal of Cell Biology*, 222(5):e202208059, 2023.

22. Paddy J Slator, Christopher W Cairo, and Nigel J Burroughs. Detection of diffusion heterogeneity in single particle tracking trajectories using a hidden markov model with measurement noise propagation. *PloS one*, 10(10):e0140759, 2015.

23. Andrew J Berglund. Statistics of camera-based single-particle tracking. *Physical Review E*, 82(1):011917, 2010.

24. Morris L Eaton and ML Eaton. *Multivariate statistics: a vector space approach*, volume 512. Wiley New York, 1983.

25. Francois Simon, Lucien E Weiss, and Sven van Teeffelen. A guide to single-particle tracking. *Nature Reviews Methods Primers*, 4(1):66, 2024.

26. Theodor Förster. Energy migration and fluorescence. *Journal of biomedical optics*, 17(1):011002–011002, 2012.

27. B Wieb van der Meer. Förster theory. *FRET-Förster Resonance Energy Transfer*, pages 23–62, 2013.

28. Björn Hellenkamp, Sonja Schmid, Olga Doroshenko, Oleg Opanasyuk, Ralf Kühnemuth, Soheila Rezaei Adariani, Benjamin Ambrose, Mikayel Aznauryan, Anders Barth, Victoria Birkedal, et al. Precision and accuracy of single-molecule fret measurements—a multi-laboratory benchmark study. *Nature methods*, 15(9):669–676, 2018.

29. Taekjip Ha, Alice Y Ting, Joy Liang, W Brett Caldwell, Ashok A Deniz, Daniel S Chemla, Peter G Schultz, and Shimon Weiss. Single-molecule fluorescence spectroscopy of enzyme conformational dynamics and cleavage mechanism. *Proceedings of the National Academy of Sciences*, 96(3):893–898, 1999.

30. Wesley B Asher, Peter Geggier, Michael D Holsey, Grant T Gilmore, Avik K Pati, Jozsef Meszaros, Daniel S Terry, Signe Mathiasen, Megan J Kaliszewski, Mitchell D McCauley, et al. Single-molecule fret imaging of gpcr dimers in living cells. *Nature methods*, 18(4):397–405, 2021.

31. Lauren Ann Metskas and Elizabeth Rhoades. Single-molecule fret of intrinsically disordered proteins. *Annual Review of Physical Chemistry*, 71(1):391–414, 2020.

32. Aritra Chowdhury, Daniel Nettels, and Benjamin Schuler. Interaction dynamics of intrinsically disordered proteins from single-molecule spectroscopy. *Annual review of biophysics*, 52(1):433–462, 2023.

33. Simo Särkkä and Angel F García-Fernández. Temporal parallelization of bayesian smoothers. *IEEE Transactions on Automatic Control*, 66(1):299–306, 2020.

34. Zitong Zhou and Daniel M Tartakovsky. Markov chain monte carlo with neural network surrogates: Application to contaminant source identification. *Stochastic Environmental Research and Risk Assessment*, 35:639–651, 2021.
35. Martin Lindén, Vladimir Curić, Elias Amselem, and Johan Elf. Pointwise error estimates in localization microscopy. *Nature communications*, 8(1):1–9, 2017.

36. Carmine Di Rienzo, Enrico Gratton, Fabio Beltram, and Francesco Cardarelli. Spatiotemporal fluctuation analysis: a powerful tool for the future nanoscopy of molecular processes. *Biophysical journal*, 111(4):679–685, 2016.

37. Raibatak Das, Christopher W Cairo, and Daniel Coombs. A hidden markov model for single particle tracks quantifies dynamic interactions between lfa-1 and the actin cytoskeleton. *PLoS computational biology*, 5(11):e1000556, 2009.

38. Koen JA Martens, Bartosz Turkowyd, Johannes Hohlbein, and Ulrike Endesfelder. Temporal analysis of relative distances (tardis) is a robust, parameter-free alternative to single-particle tracking. *Nature Methods*, pages 1–8, 2024.

39. Joshua V Dillon, Ian Langmore, Dustin Tran, Eugene Brevdo, Srinivas Vasudevan, Dave Moore, Brian Patton, Alex Alemi, Matt Hoffman, and Rif A Saurous. Tensorflow distributions. *arXiv preprint arXiv:1711.10604*, 2017.

40. Gorka Muñoz-Gil, Giovanni Volpe, Miguel Angel Garcia-March, Erez Aghion, Aykut Argun, Chang Beom Hong, Tom Bland, Stefano Bo, J Alberto Conejero, Nicolás Firbas, et al. Objective comparison of methods to decode anomalous diffusion. *Nature communications*, 12(1):6253, 2021.

41. Makram Chahine, Ramin Hasani, Patrick Kao, Aaron Ray, Ryan Shubert, Mathias Lechner, Alexander Amini, and Daniela Rus. Robust flight navigation out of distribution with liquid neural networks. *Science Robotics*, 8(77):eadc8892, 2023.

# Appendix

## 1 Inferring the hidden variables

The high degree of structure present in the conditional probabilities leads to a specific form of covariance matrix that can occur. To effectively discuss this, instead of $x$ and $y$ to represent the observed and hidden variables consider $v = (v^1 v^2 \dots)^T = (x^1 x^2 \dots y^1 y^2 \dots)^T$. To index by time a subscript will be used instead; for example $v_1^2 = x_1^2$ and $v_1^{n+2} = y_1^2$. This allows the equation 4 to be written as:

$$L(\theta|X) = \int_Y I(y_1) \prod_{i=1}^{N} \prod_{j=1}^{n+m} f_i \left( \sum_{j=1}^{n+m} a_{ij} v_i^j \right) dY \quad (16)$$

Carrying out the product allows the integrand to be expressed as a single multivariate Gaussian. Writing $u = (v_1 v_2 \dots)$ and $d_i$ as the variance associated with Gaussian $f_i$, there exists some $(m+n)N \times (m+n)N$ matrix $\Sigma$ such that:

$$L(\theta|X) = \int_Y I(y_1) \left( \prod_{i=1}^{N} \prod_{j=1}^{n+m} \frac{1}{\sqrt{2\pi d_i}} \right) \exp \left( u^T \Sigma u \right) dY \quad (17)$$

The goal of expressing the integrand in this form is that it allows for the expression of the conditional likelihood of the hidden variables, given the observed variables. In particular, by separating the covariance matrix into blocks that pertain to $x$ and $y$, one may express the integrated Gaussian as:

$$\frac{1}{\sqrt{(2\pi)^{N(n+m)}|\Sigma|}} \exp \left( -\frac{1}{2} \begin{pmatrix} x \\ y \end{pmatrix}^T \begin{pmatrix} \Sigma_{xx} & \Sigma_{xy} \\ \Sigma_{yx} & \Sigma_{yy} \end{pmatrix}^{-1} \begin{pmatrix} x \\ y \end{pmatrix} \right)$$

Using this notation, $Y|x \sim N(\Sigma_{YX}\Sigma_{XX}X, \Sigma_{YY} - \Sigma_{YX}\Sigma_{XX}^{-1}\Sigma_{YY})$.

Proof: This can be seen by starting from the definition of conditional probability and substituting the above likelihood.

$$p(y|x) = \frac{p(x, y)}{p(x)} = \frac{1}{\sqrt{(2\pi)^{N_n} |\Sigma| |\Sigma_{XX}|}} \exp \left( -\frac{1}{2} \begin{pmatrix} x \\ y \end{pmatrix}^T \begin{pmatrix} \Sigma_{XX} & \Sigma_{XY} \\ \Sigma_{YX} & \Sigma_{YY} \end{pmatrix}^{-1} \begin{pmatrix} x \\ y \end{pmatrix} - x^T \Sigma_{XX}^{-1} x \right)$$

The elements of $\Sigma^{-1}$ can be found by Schur complements. Since these expressions are quite lengthy, the block form of $\Sigma^{-1}$ will be abbreviated as $\begin{pmatrix} C_{XX} & C_{XY} \\ C_{YX} & C_{YY} \end{pmatrix}$ and the specific terms will be written as needed. In this notation the exponent can be written as:

$$-\frac{1}{2} (y^T C_{YY} y + 2y^T C_{XY} X + x^T (C_{YY} - \Sigma_{XX}^{-1}) x)$$
The next step is to eliminate the $\Sigma_{XX}$ term in favour of C terms, first we have:

$$ 
\begin{aligned}
C_{YY} - \Sigma_{XX}^{-1} &= \Sigma_{XX}^{-1} + \Sigma_{XX}^{-1} \Sigma_{XY} (\Sigma_{YY} - \Sigma_{YX} \Sigma_{XX}^{-1} \Sigma_{XY})^{-1} \Sigma_{YX} \Sigma_{XX}^{-1} - \Sigma_{XX}^{-1} \\
&= -\Sigma_{XX}^{-1} (-\Sigma_{XY} (\Sigma_{YY} - \Sigma_{YX} \Sigma_{XX}^{-1} \Sigma_{XY})^{-1} \Sigma_{YX} \Sigma_{XX}^{-1}) \\
&= -\Sigma_{XX}^{-1} (-\Sigma_{XY} (\Sigma_{YY} - \Sigma_{YX} \Sigma_{XX}^{-1} \Sigma_{XY}) (\Sigma_{YY} - \Sigma_{YX} \Sigma_{XX}^{-1} \Sigma_{XY})^{-1} (\Sigma_{YY} - \Sigma_{YX} \Sigma_{XX}^{-1} \Sigma_{XY})^{-1} \Sigma_{YX} \Sigma_{XX}^{-1}) \\
&= C_{XY} C_{YY}^{-1} C_{YX}
\end{aligned}
$$

Using this, the exponent can be written out as:

$$ 
\begin{aligned}
& -\frac{1}{2} (y^T C_{YY} y + 2y^T C_{XY} x + x^T C_{XY} C_{YY}^{-1} C_{XY} x) \\
&= -\frac{1}{2} (y^T C_{YY} y + 2y^T C_{YY} C_{YY}^{-1} C_{YX} x + x^T C_{YX} C_{YY}^{-T} C_{YY} C_{YY}^{-1} C_{YX} x) \\
&= -\frac{1}{2} ((y + C_{YY}^{-1} C_{YX} x)^T C_{YY} (y + C_{YY}^{-1} C_{YX} x))
\end{aligned}
$$

By referring to the Schur complement decomposition of $\Sigma^{-1}$ the claim is verified. From collected data, the $\Sigma_{XX}$ portion can be obtained by analyzing statistics, the other components need a model that relates them to the quantities in the XX section to compute the $\Sigma_{YY}$ and $\Sigma_{YY}$. Once the new mean and covariance matrix are obtained, a Givens rotation QR decomposition provides an $O((Nm)^3)$ method for finding the determinant. If a banded cutoff is used, with b bands, this can be reduced to $O(Nmb^2)$.