2.1. Vector Autoregressive models

A VAR model can be used as a parametrized representation of a multichannel discrete-time series $ \boldsymbol{Y}\in {\unicode{x211D}}^{m\times N} $ for which stationarity of $ \boldsymbol{Y} $ is a fundamental requirement, implying that the model parameters are time-invariant over the analysis period (Farrar and Worden, Reference Farrar and Worden2012, p. 233). The VAR model is an extension of the conventional AR model, enabling the joint modeling of multiple, potentially coupled signals. Modeling each signal independently may neglect important cross-channel dependencies, whereas the VAR formulation captures both dynamics within each channel and the coupling between channels. The VAR model is formulated as:

(1) $$ \boldsymbol{y}\left[t\right]=-\sum \limits_{i=1}^{n_a}{\boldsymbol{A}}_i\cdot \boldsymbol{y}\left[t-i\right]+\boldsymbol{w}\left[t\right],\boldsymbol{w}\left[t\right]\sim \mathcal{N}\left(\mathbf{0},{\boldsymbol{\Sigma}}_w\right) $$

where $ \boldsymbol{y}\left[t\right]\in {\unicode{x211D}}^m $ is a column vector of multichannel observations in $ \boldsymbol{Y} $ at time $ t $ . The AR parameters at the given time lag $ i $ are represented in the matrix $ {\boldsymbol{A}}_i\in {\unicode{x211D}}^{m\times m} $ . The vector $ \boldsymbol{w}\left[t\right]\in {\unicode{x211D}}^m $ denotes the zero-mean normally and independently distributed innovations, which represent the part of the stochastic process that cannot be predicted from the previous $ {n}_a $ samples. These innovations are distributed with the innovation covariance matrix $ {\boldsymbol{\Sigma}}_w\in {\unicode{x211D}}^{m\times m} $ (Farrar and Worden, Reference Farrar and Worden2012, Ch. 7.11). This structure can be compactly written for the complete time-series from the $ {n}_a+1 $ observation to the $ {N}^{\mathrm{th}} $ observation, denoted as $ {\boldsymbol{Y}}_{\tau } $ , as:

(2) $$ {\boldsymbol{Y}}_{\tau }=\boldsymbol{\Theta} \boldsymbol{\Phi} +\boldsymbol{W} $$

where the innovations in $ \boldsymbol{W}=\left[\boldsymbol{w}\left[{n}_a+1\right]\hskip0.5em \boldsymbol{w}\left[{n}_a+2\right]\hskip0.5em \cdots \hskip0.5em \boldsymbol{w}\left[N\right]\right] $ are independent across $ t $ . The complete set of AR parameters is stacked in the matrix $ \boldsymbol{\Theta} \in {\unicode{x211D}}^{m\times d} $ with $ d={mn}_a $ , and the lagged observations are stacked in the regression matrix $ \boldsymbol{\Phi} \in {\unicode{x211D}}^{d\times \tau } $ where $ \tau =N-{n}_a $ such that:

(3) $$ \boldsymbol{\Theta} =\left[\begin{array}{cccc}{\boldsymbol{A}}_1& {\boldsymbol{A}}_2& \cdots & {\boldsymbol{A}}_{n_a}\end{array}\right]\mathrm{and}\boldsymbol{\Phi } =-\left[\begin{array}{cccc}\boldsymbol{y}\left[{n}_a\right]& \boldsymbol{y}\left[{n}_a+1\right]& \cdots & \boldsymbol{y}\left[N-1\right]\\ {}\boldsymbol{y}\left[{n}_a-1\right]& \boldsymbol{y}\left[{n}_a\right]& \cdots & \boldsymbol{y}\left[N-2\right]\\ {}\vdots & \vdots & \vdots & \vdots \\ {}\boldsymbol{y}\left[1\right]& \boldsymbol{y}\left[2\right]& \cdots & \boldsymbol{y}\left[N-{n}_a\right]\end{array}\right] $$

Classically, estimation of the AR parameters $ \boldsymbol{\Theta} $ and innovations covariance $ {\boldsymbol{\Sigma}}_w $ is performed by means of a Least Squares approach, where the objective is to minimize the sum of squared differences between the observed values and the model predictions (Hastie et al., Reference Hastie, Tibshirani and Friedman2009, pp. 11–12):

(4) $$ \underset{\boldsymbol{\Theta}}{\min }{\left\Vert {\boldsymbol{Y}}_{\tau }-\boldsymbol{\Theta} \boldsymbol{\Phi} \right\Vert}^2 $$

which leads to the solution:

(5) $$ \hat{\boldsymbol{\Theta}}={\boldsymbol{Y}}_{\tau }{\boldsymbol{\Phi}}^T{\left({\boldsymbol{\Phi} \boldsymbol{\Phi}}^T\right)}^{-1} $$

where the innovations can be estimated as:

(6) $$ \hat{\boldsymbol{W}}={\boldsymbol{Y}}_{\tau }-\hat{\boldsymbol{\Theta}}\boldsymbol{\Phi} $$

and the innovation covariance as:

(7) $$ {\hat{\boldsymbol{\Sigma}}}_w=\operatorname{cov}\left(\hat{\boldsymbol{w}}\left[t\right]\right) $$

This estimation assumes that $ \boldsymbol{\Theta} $ and $ {\boldsymbol{\Sigma}}_w $ are fixed, deterministic values, implicitly requiring that the available data are sufficient to estimate them accurately and that no parameter uncertainty remains—an assumption that is rarely satisfied in real-life applications. To account for this uncertainty, $ \boldsymbol{\Theta} $ and $ {\boldsymbol{\Sigma}}_w $ may instead be treated as random variables. However, this requires alternative estimation methods.

2.1.1. Bayesian-formulated VAR models

While the least squares solution in Equations (5)–(7) treats $ \boldsymbol{\Theta} $ and $ {\boldsymbol{\Sigma}}_w $ as fixed quantities, a Bayesian formulation recognizes that, in practice, these parameters are uncertain—particularly when they are estimated from short data segments, when operating conditions fluctuate, or when the signal-to-noise ratio is low. Under such conditions, the ordinary Least Squares estimator is prone to overfitting, producing parameter estimates that adapt too closely to noise rather than to the underlying system dynamics. By modeling $ \boldsymbol{\Theta} $ and $ {\boldsymbol{\Sigma}}_w $ as random variables, the Bayesian approach explicitly accounts for this estimation uncertainty and introduces prior structure that stabilizes the inference, especially in regimes where data are limited or highly variable. In the Bayesian VAR formulation, each observation $ \boldsymbol{y}\left[t\right] $ is conditioned on the model parameters $ \boldsymbol{\Theta} $ and $ {\boldsymbol{\Sigma}}_w $ , which are themselves random with hyperparameters $ \boldsymbol{\beta} $ , as well as on the lagged observations $ \boldsymbol{\phi} \left[t\right] $ corresponding to the $ {t}^{\mathrm{th}} $ column in $ \boldsymbol{\Phi} $ . This conditioned probability is therefore expressed as:

(8) $$ p\left(\boldsymbol{y}\left[t\right]\;|\;\boldsymbol{\Theta}, {\boldsymbol{\Sigma}}_w,\boldsymbol{\phi} \left[t\right]\right) $$

with parameter priors conditioned as:

(9) $$ p\left(\boldsymbol{\Theta}, {\boldsymbol{\Sigma}}_w\;|\;\boldsymbol{\beta} \right) $$

This captures that $ \boldsymbol{y}\left[t\right] $ depends both on the lagged observations $ \boldsymbol{\phi} \left[t\right] $ and on the random parameters $ \boldsymbol{\Theta} $ and, $ {\boldsymbol{\Sigma}}_w $ as illustrated in the Bayesian network in Figure 2.

Figure 2.

Bayesian network of the VAR model.

The likelihood of the complete time-series from the $ {n}_a+1 $ observation to the $ {N}^{\mathrm{th}} $ observation, $ {\boldsymbol{Y}}_{\tau } $ , is obtained by assuming conditional independence of the observations given the parameters:

(10) $$ p\left({\boldsymbol{Y}}_{\tau}\;|\;\boldsymbol{\Theta}, {\boldsymbol{\Sigma}}_w,\boldsymbol{\Phi} \right)=\prod \limits_{t={n}_a+1}^Np\left(\boldsymbol{y}\left[t\right]\;|\;\boldsymbol{\Theta}, {\boldsymbol{\Sigma}}_w,\boldsymbol{\phi} \left[t\right]\right) $$

By Bayes’ theorem (Murphy, Reference Murphy2012, p. 29), the joint posterior distribution of the parameters is then:

(11) $$ p\left(\boldsymbol{\Theta}, {\boldsymbol{\Sigma}}_w\;|\;{\boldsymbol{Y}}_{\tau },\boldsymbol{\Phi} \right)=\frac{p\left({\boldsymbol{Y}}_{\tau}\;|\;\boldsymbol{\Theta}, {\boldsymbol{\Sigma}}_w,\boldsymbol{\Phi} \right)p\left(\boldsymbol{\Theta}\;|\;{\boldsymbol{\Sigma}}_w\right)p\left({\boldsymbol{\Sigma}}_w\right)}{p\left({\boldsymbol{Y}}_{\tau}\;|\;\boldsymbol{\Phi} \right)} $$

where the dependence on the hyperparameter $ \boldsymbol{\beta} $ is implicitly assumed, and the denominator is the marginal likelihood acting as a normalizing constant.

Calculating the joint posterior in Equation (11) is often not feasible, and simulation-based methods are typically employed. Under the assumption of conjugate priors, where the product of the prior and the likelihood belongs to the same family as the posterior, analytical expressions for the posterior can be obtained, which simplifies the inference process (Gelman, Reference Gelman2013, p. 36). For this, it is assumed that the likelihood $ p\left({\boldsymbol{Y}}_{\tau}\;|\;\boldsymbol{\Theta}, {\boldsymbol{\Sigma}}_w,\boldsymbol{\Phi} \right) $ and the prior $ p\left(\boldsymbol{\Theta}\;|\;{\boldsymbol{\Sigma}}_w\right) $ are Gaussian, each modeled using the Matrix Normal ( $ \mathcal{M}\mathcal{N} $ ) distribution, and that the prior $ p\left({\boldsymbol{\Sigma}}_w\right) $ is Inverse-Wishart ( $ \mathcal{I}\mathcal{W} $ ) as follows:

(12) $$ {\boldsymbol{Y}}_{\tau} \mid \boldsymbol{\Theta}, {\boldsymbol{\Sigma}}_w,\boldsymbol{\Phi} \sim \mathcal{M}\mathcal{N}\left(\boldsymbol{\Theta} \boldsymbol{\Phi}, {\boldsymbol{\Sigma}}_w,\boldsymbol{I}\right) $$

(13) $$ \boldsymbol{\Theta} \mid {\boldsymbol{\Sigma}}_w\sim \mathcal{M}\mathcal{N}\left(\boldsymbol{\mu}, {\boldsymbol{\Sigma}}_w,\boldsymbol{V}\right) $$

(14) $$ {\boldsymbol{\Sigma}}_w\sim \mathcal{I}\mathcal{W}\left(\boldsymbol{\Psi}, \nu \right) $$

where $ \boldsymbol{I}\in {\unicode{x211D}}^{\tau \times \tau } $ is the identity matrix, the mean $ \boldsymbol{\mu} \in {\unicode{x211D}}^{m\times d} $ and the covariance $ \boldsymbol{V}\in {\unicode{x211D}}^{d\times d} $ are the hyperparameters of the prior on $ \boldsymbol{\Theta} $ , and the scale matrix $ \boldsymbol{\Psi} \in {\unicode{x211D}}^{m\times m} $ and degrees of freedom $ \nu \in \unicode{x211D} $ are the hyperparameters of the prior on $ {\boldsymbol{\Sigma}}_w $ . With this, the prior hyperparameter $ \boldsymbol{\beta} $ in Figure 2 consists of $ \left\{\boldsymbol{\mu}, \boldsymbol{V},\boldsymbol{\Psi}, \nu \right\} $ . Under these assumptions, conjugacy holds, and the posterior can be obtained analytically (Gelman, Reference Gelman2013, p. 584; Avendaño-Valencia, Reference Avendaño-Valencia2025). With the definitions in Equations. (12), (13), and (14), the posterior distribution in Equation (11) can be written as Gupta and Nagar (Reference Gupta and Nagar1999, Ch. 2–3):

(15) $$ {\displaystyle \begin{array}{l}p\left(\boldsymbol{\Theta}, {\boldsymbol{\Sigma}}_w|{\boldsymbol{Y}}_{\tau },\boldsymbol{\Phi} \right)\propto {\left|{\boldsymbol{\Sigma}}_w\right|}^{-\left(\nu +m+1\right)/2}\exp \left(-\frac{1}{2}\hskip0.1em \mathrm{tr}\left\{\boldsymbol{\Psi} \hskip0.1em {\boldsymbol{\Sigma}}_w^{-1}\right\}\right)\\ {}\times {\left|{\boldsymbol{\Sigma}}_w\right|}^{-d/2}\exp \left(-\frac{1}{2}\hskip0.1em \mathrm{tr}\left\{{\boldsymbol{\Sigma}}_w^{-1}\left(\boldsymbol{\Theta} -\boldsymbol{\mu} \right)\hskip0.1em {\boldsymbol{V}}^{-1}{\left(\boldsymbol{\Theta} -\boldsymbol{\mu} \right)}^{\top}\right\}\right)\\ {}\times {\left|{\boldsymbol{\Sigma}}_w\right|}^{-\tau /2}\exp \left(-\frac{1}{2}\hskip0.1em \mathrm{tr}\left\{{\boldsymbol{\Sigma}}_w^{-1}\left({\boldsymbol{Y}}_{\tau }-\boldsymbol{\Theta} \boldsymbol{\Phi} \right){\left({\boldsymbol{Y}}_{\tau }-\boldsymbol{\Theta} \boldsymbol{\Phi} \right)}^{\top}\right\}\right)\end{array}} $$

Expanding the third term, the likelihood trace, gives:

(16) $$ \mathrm{tr}\left\{{\boldsymbol{\Sigma}}_w^{-1}\left({\boldsymbol{Y}}_{\tau }-\boldsymbol{\Theta} \boldsymbol{\Phi} \right){\left({\boldsymbol{Y}}_{\tau }-\boldsymbol{\Theta} \boldsymbol{\Phi} \right)}^{\top}\right\}=\mathrm{tr}\left\{{\boldsymbol{\Sigma}}_w^{-1}{\boldsymbol{Y}}_{\tau }{\boldsymbol{Y}}_{\tau}^{\top}\right\}-2\hskip0.1em \mathrm{tr}\left\{{\boldsymbol{\Sigma}}_w^{-1}{\boldsymbol{Y}}_{\tau }{\boldsymbol{\Phi}}^{\top }{\boldsymbol{\Theta}}^{\top}\right\}+\mathrm{tr}\left\{{\boldsymbol{\Sigma}}_w^{-1}\boldsymbol{\Theta} \hskip0.1em {\boldsymbol{\Phi} \boldsymbol{\Phi}}^{\top }{\boldsymbol{\Theta}}^{\top}\right\}. $$

Hence, the $ \boldsymbol{\Theta} $ -dependent part of Equation (15) can be grouped as:

(17) $$ \mathrm{tr}\left\{{\boldsymbol{\Sigma}}_w^{-1}\left(\boldsymbol{\Theta} \left({\boldsymbol{\Phi} \boldsymbol{\Phi}}^{\top }+{\boldsymbol{V}}^{-1}\right){\boldsymbol{\Theta}}^{\top }-2\left({\boldsymbol{Y}}_{\tau }{\boldsymbol{\Phi}}^{\top }+\boldsymbol{\mu} {\boldsymbol{V}}^{-1}\right){\boldsymbol{\Theta}}^{\top}\right)\right\} $$

Let

(18) $$ {\boldsymbol{V}}_{\mathrm{post}}={\left({\boldsymbol{V}}^{-1}+{\boldsymbol{\Phi} \boldsymbol{\Phi}}^{\top}\right)}^{-1}\mathrm{and}{\boldsymbol{\mu}}_{\mathrm{post}}=\left(\boldsymbol{\mu} {\boldsymbol{V}}^{-1}+{\boldsymbol{Y}}_{\tau }{\boldsymbol{\Phi}}^{\top}\right){\boldsymbol{V}}_{\mathrm{post}} $$

Then the $ \boldsymbol{\Theta} $ -dependent part can be expressed as:

(19) $$ {\displaystyle \begin{array}{l}\mathrm{tr}\left\{{\boldsymbol{\Sigma}}_w^{-1}\left(\boldsymbol{\Theta} \left({\boldsymbol{\Phi} \boldsymbol{\Phi}}^{\top }+{\boldsymbol{V}}^{-1}\right){\boldsymbol{\Theta}}^{\top }-2\left(\boldsymbol{\mu} {\boldsymbol{V}}^{-1}+{\boldsymbol{Y}}_{\tau }{\boldsymbol{\Phi}}^{\top}\right){\boldsymbol{\Theta}}^{\top}\right)\right\}=\mathrm{tr}\left\{{\boldsymbol{\Sigma}}_w^{-1}\left(\boldsymbol{\Theta} -{\boldsymbol{\mu}}_{\mathrm{post}}\right){\boldsymbol{V}}_{\mathrm{post}}^{-1}{\left(\boldsymbol{\Theta} -{\boldsymbol{\mu}}_{\mathrm{post}}\right)}^{\top}\right\}\\ {}-\mathrm{tr}\left\{{\boldsymbol{\Sigma}}_w^{-1}{\boldsymbol{\mu}}_{\mathrm{post}}{\boldsymbol{V}}_{\mathrm{post}}^{-1}{\boldsymbol{\mu}}_{\mathrm{post}}^{\top}\right\}\end{array}} $$

Substituting this into Equation (15) and grouping the $ \boldsymbol{\Theta} $ -dependent and $ \boldsymbol{\Theta} $ -independent parts results in:

(20) $$ {\displaystyle \begin{array}{l}p\left(\boldsymbol{\Theta}, {\boldsymbol{\Sigma}}_w|{\boldsymbol{Y}}_{\tau },\boldsymbol{\Phi} \right)\propto {\left|{\boldsymbol{\Sigma}}_w\right|}^{-\left(\nu +\tau +m+1\right)/2}\hskip0.1em \exp \left(-\frac{1}{2}\hskip0.1em \mathrm{tr}\left\{{\boldsymbol{\Psi}}_{\mathrm{post}}{\boldsymbol{\Sigma}}_w^{-1}\right\}\right)\\ {}\times {\left|{\boldsymbol{\Sigma}}_w\right|}^{-d/2}\hskip0.1em \exp \left(-\frac{1}{2}\hskip0.1em \mathrm{tr}\left\{{\boldsymbol{\Sigma}}_w^{-1}\left(\boldsymbol{\Theta} -{\boldsymbol{\mu}}_{\mathrm{post}}\right)\hskip0.1em {\boldsymbol{V}}_{\mathrm{post}}^{-1}{\left(\boldsymbol{\Theta} -{\boldsymbol{\mu}}_{\mathrm{post}}\right)}^{\top}\right\}\right)\end{array}} $$

with updates

(21) $$ {\boldsymbol{V}}_{\mathrm{post}}={\left({\boldsymbol{V}}^{-1}+{\boldsymbol{\Phi} \boldsymbol{\Phi}}^{\top}\right)}^{-1} $$

(22) $$ {\boldsymbol{\mu}}_{\mathrm{post}}=\left(\boldsymbol{\mu} \hskip0.1em {\boldsymbol{V}}^{-1}+{\boldsymbol{Y}}_{\tau }{\boldsymbol{\Phi}}^{\top}\right)\hskip0.1em {\boldsymbol{V}}_{\mathrm{post}} $$

(23) $$ {\nu}_{\mathrm{post}}=\nu +\tau $$

(24) $$ {\boldsymbol{\Psi}}_{\mathrm{post}}=\boldsymbol{\Psi} \left({\boldsymbol{Y}}_{\tau }-{\boldsymbol{\mu}}_{\mathrm{post}}\boldsymbol{\Phi} \right){\left({\boldsymbol{Y}}_{\tau }-{\boldsymbol{\mu}}_{\mathrm{post}}\boldsymbol{\Phi} \right)}^{\top}\left({\boldsymbol{\mu}}_{\mathrm{post}}-\boldsymbol{\mu} \right)\hskip0.1em {\boldsymbol{V}}^{-1}\hskip0.1em {\left({\boldsymbol{\mu}}_{\mathrm{post}}-\boldsymbol{\mu} \right)}^{\top } $$

which constitute the posterior hyperparameters after the inference.

The posterior hyperparameters can be inferred in two ways. The first is the batch update, in which the prior hyperparameters are defined once, and all available observations are incorporated simultaneously, yielding the posterior hyperparameters in Equations (21)–(24). The second is the sequential update, where the posterior after each batch of data is reused as the prior when new observations arrive (Avendaño-Valencia, Reference Avendaño-Valencia2025). In practice, this means that each new segment of data contributes an incremental update to the hyperparameters, allowing the model to gradually refine its description of the underlying dynamics without revisiting past data. This recursive procedure is beneficial when data arrive incrementally. However, it implicitly assumes that all new observations are consistent with the underlying dynamics, which, in practice, may be violated due to transient deviations and non-stationary behavior. As a result, directly applying the recursive update in Equations (21)–(24) may lead to abrupt changes in the inferred hyperparameters. To mitigate this, a regularized update can be adopted, in which the posterior hyperparameters are combined with the previous prior through a weighted update governed by forgetting factors $ {\eta}_{\Sigma},{\eta}_{\Theta}\in \left[0,1\right] $ :

(25) $$ {\tilde{\boldsymbol{V}}}_{\mathrm{post}}=\left(1-{\eta}_{\Theta}\right)\boldsymbol{V}+{\eta}_{\Theta}{\boldsymbol{V}}_{\mathrm{post}} $$

(26) $$ {\tilde{\boldsymbol{\mu}}}_{\mathrm{post}}=\left(1-{\eta}_{\Theta}\right)\boldsymbol{\mu} +{\eta}_{\Theta}{\boldsymbol{\mu}}_{\mathrm{post}} $$

(27) $$ {\tilde{\nu}}_{\mathrm{post}}=\left(1-{\eta}_{\Sigma}\right)\nu +{\eta}_{\Sigma}{\nu}_{\mathrm{post}} $$

(28) $$ {\tilde{\boldsymbol{\Psi}}}_{\mathrm{post}}=\left(1-{\eta}_{\Sigma}\right)\boldsymbol{\Psi} +{\eta}_{\Sigma}{\boldsymbol{\Psi}}_{\mathrm{post}} $$

controlling the influence of newly observed data relative to previously inferred information (Prando et al., Reference Prando, Romeres and Chiuso2016; Søndergaard et al., Reference Søndergaard, Shaker and Nørregaard Jørgensen2024; Prasanthi et al., Reference Prasanthi, Shareef, Khalid and Selvaraj2025), which can be determined either heuristically as fixed weights or adaptively. Both the batch update and sequential update strategies rely on the same conjugate structure, ensuring analytical tractability.

Marginalizing $ \boldsymbol{\Theta} $ and $ {\boldsymbol{\Sigma}}_w $ .

The model parameters $ \boldsymbol{\Theta} $ and $ {\boldsymbol{\Sigma}}_w $ are represented as posterior distributions because they are treated as random variables. Further analysis often involves expectations with respect to these posteriors, which in general are approximated by sampling-based methods, but this can be computationally inefficient. In this case of conjugacy, however, both $ \boldsymbol{\Theta} $ and $ {\boldsymbol{\Sigma}}_w $ can be marginalized analytically, yielding closed-form expressions for the marginal likelihood and avoiding the need for sampling.

Starting from the joint model in Equations (12)–(14), the first step is to marginalize $ \boldsymbol{\Theta} $ , resulting in the partly marginalized likelihood of $ {\boldsymbol{Y}}_{\tau } $ :

(29) $$ p\left({\boldsymbol{Y}}_{\tau }|{\boldsymbol{\Sigma}}_w,\boldsymbol{\Phi} \right)=\int p\left({\boldsymbol{Y}}_{\tau }|\boldsymbol{\Theta}, {\boldsymbol{\Sigma}}_w,\boldsymbol{\Phi} \right)\hskip0.1em p\left(\boldsymbol{\Theta} |{\boldsymbol{\Sigma}}_w\right)\hskip0.1em d\boldsymbol{\Theta} $$

In its vectorized form, this is similar to a standard Gaussian model $ \boldsymbol{y}=\boldsymbol{Ax}+\boldsymbol{b} $ in Murphy (Reference Murphy2012, p. 119), in which:

(30) $$ p\left(\boldsymbol{x}\right)=\mathcal{N}\left({\boldsymbol{\mu}}_x,{\boldsymbol{\Sigma}}_x\right) $$

(31) $$ p\left(\boldsymbol{y}|\boldsymbol{x}\right)=\mathcal{N}\left(\boldsymbol{Ax}+\boldsymbol{b},{\boldsymbol{\Sigma}}_y\right) $$

which is equivalent to Equation (29) with the following:

(32) $$ {\displaystyle \begin{array}{l}\boldsymbol{y}={\overline{\boldsymbol{Y}}}_{\tau },\boldsymbol{x}=\overline{\boldsymbol{\Theta}},\boldsymbol{A}={\boldsymbol{\Phi}}^{\top}\otimes {\boldsymbol{I}}_{m\times m},\\ {}\boldsymbol{b}=0,{\boldsymbol{\Sigma}}_x=\boldsymbol{V}\otimes {\boldsymbol{\Sigma}}_w,{\boldsymbol{\Sigma}}_y={\boldsymbol{I}}_{\tau \times \tau}\otimes {\boldsymbol{\Sigma}}_w\end{array}} $$

Following the definition in Murphy (Reference Murphy2012, p. 120), the partly marginalized likelihood of $ {\boldsymbol{Y}}_{\tau } $ can be expressed as:

(33) $$ {\overline{\boldsymbol{Y}}}_{\tau}\mid {\boldsymbol{\Sigma}}_w,\boldsymbol{\Phi} \sim \mathcal{N}\left(\left({\boldsymbol{\Phi}}^{\top}\otimes {\boldsymbol{I}}_{m\times m}\right)\mathrm{vec}\left(\boldsymbol{\mu} \right),\left({\boldsymbol{I}}_{\tau \times \tau}\otimes {\boldsymbol{\Sigma}}_w\right)+\left({\boldsymbol{\Phi}}^{\top}\otimes {\boldsymbol{I}}_{m\times m}\right)\left(\boldsymbol{V}\otimes {\boldsymbol{\Sigma}}_w\right){\left({\boldsymbol{\Phi}}^{\top}\otimes {\boldsymbol{I}}_{m\times m}\right)}^{\top}\right) $$

and by following the Kronecker product of ordinary products (Harville, Reference Harville2006, p. 341), part of the covariance can be simplified:

(34) $$ \left({\boldsymbol{\Phi}}^{\top}\otimes {\boldsymbol{I}}_{m\times m}\right)\left(\boldsymbol{V}\otimes {\boldsymbol{\Sigma}}_w\right){\left({\boldsymbol{\Phi}}^{\top}\otimes {\boldsymbol{I}}_{m\times m}\right)}^{\top }=\left({\boldsymbol{\Phi}}^{\top}\boldsymbol{V}\boldsymbol{\Phi } \right)\otimes {\boldsymbol{\Sigma}}_w $$

so the complete covariance in Equation (33) simplifies into:

(35) $$ \left({\boldsymbol{I}}_{\tau \times \tau}\otimes {\boldsymbol{\Sigma}}_w\right)+\left({\boldsymbol{\Phi}}^{\top}\otimes {\boldsymbol{I}}_{m\times m}\right)\left(\boldsymbol{V}\otimes {\boldsymbol{\Sigma}}_w\right){\left(\boldsymbol{\Phi} \otimes {\boldsymbol{I}}_{m\times m}\right)}^{\top }=\left({\boldsymbol{I}}_{\tau \times \tau }+{\boldsymbol{\Phi}}^{\top}\boldsymbol{V}\boldsymbol{\Phi } \right)\otimes {\boldsymbol{\Sigma}}_w $$

Reshaping back to Matrix–Normal form (Harville, Reference Harville2006, p. 345) yields:

(36) $$ {\boldsymbol{Y}}_{\tau}\mid {\boldsymbol{\Sigma}}_w,\boldsymbol{\Phi} \sim \mathcal{M}\mathcal{N}\left(\boldsymbol{\mu} \boldsymbol{\Phi}, {\boldsymbol{\Sigma}}_w,{\boldsymbol{I}}_{\tau \times \tau }+{\boldsymbol{\Phi}}^{\top}\boldsymbol{V}\boldsymbol{\Phi } \right) $$

which is Equation (29) and expresses the partly marginalized likelihood of $ {\boldsymbol{Y}}_{\tau } $ given $ {\boldsymbol{\Sigma}}_w $ and $ \boldsymbol{\Phi} $ . To obtain the complete marginalized likelihood of $ {\boldsymbol{Y}}_{\tau } $ , the innovation covariance $ {\boldsymbol{\Sigma}}_w $ can also be marginalized as:

(37) $$ p\left({\boldsymbol{Y}}_{\tau }|\boldsymbol{\Phi} \right)=\int p\left({\boldsymbol{Y}}_{\tau }|{\boldsymbol{\Sigma}}_w,\boldsymbol{\Phi} \right)\hskip0.1em p\left({\boldsymbol{\Sigma}}_w\right)\hskip0.1em d{\boldsymbol{\Sigma}}_w $$

From Equation (36), the partly marginalized likelihood is Matrix–Normal and the prior on $ {\boldsymbol{\Sigma}}_w $ is Inverse–Wishart as in Equation (14). In vectorized form, the partly marginalized likelihood corresponds to:

(38) $$ {\overline{\boldsymbol{Y}}}_{\tau}\mid {\boldsymbol{\Sigma}}_w,\boldsymbol{\Phi} \sim \mathcal{N}\left(\left({\boldsymbol{\Phi}}^{\top}\otimes {\boldsymbol{I}}_{m\times m}\right)\mathrm{vec}\left(\boldsymbol{\mu} \right),\left({\boldsymbol{I}}_{\tau \times \tau }+{\boldsymbol{\Phi}}^{\top}\boldsymbol{V}\boldsymbol{\Phi } \right)\otimes {\boldsymbol{\Sigma}}_w\right) $$

and, combined with the Inverse-Wishart prior, the joint distribution $ p\left({\overline{\boldsymbol{Y}}}_{\tau },{\boldsymbol{\Sigma}}_w|\boldsymbol{\Phi} \right) $ belongs to the Normal-Inverse-Wishart ( $ \mathcal{N}\mathcal{I}\mathcal{W} $ ) family (Murphy, Reference Murphy2012, pp. 132–133) expressed as:

(39) $$ {\overline{\boldsymbol{Y}}}_{\tau },{\boldsymbol{\Sigma}}_w\mid \boldsymbol{\Phi} \sim \mathcal{N}\mathcal{I}\mathcal{W}\left(\left({\boldsymbol{\Phi}}^{\top}\otimes {\boldsymbol{I}}_m\right)\mathrm{vec}\left(\boldsymbol{\mu} \right),{\boldsymbol{I}}_{\tau }+{\boldsymbol{\Phi}}^{\top}\boldsymbol{V}\boldsymbol{\Phi }, \boldsymbol{\Psi}, \nu \right) $$

Consequently, the complete marginalized likelihood in Equation (37) is a multivariate Student- $ t $ distribution for $ {\overline{\boldsymbol{Y}}}_{\tau } $ (Murphy, Reference Murphy2012, p. 134). Reshaping back into matrix form gives a Matrix–Student– $ t $ distribution:

(40) $$ {\boldsymbol{Y}}_{\tau}\mid \boldsymbol{\Phi} \sim \mathcal{M}\mathcal{T}\left(\boldsymbol{\mu} \boldsymbol{\Phi}, \boldsymbol{\Psi}, {\boldsymbol{I}}_{\tau \times \tau }+{\boldsymbol{\Phi}}^{\top}\boldsymbol{V}\boldsymbol{\Phi } \right) $$

The corresponding log marginal likelihood is written as Gupta and Nagar (Reference Gupta and Nagar1999, pp. 133–134):

(41) $$ {\displaystyle \begin{array}{l}\log p\left({\boldsymbol{Y}}_{\tau }|\boldsymbol{\Phi} \right)=\log {\Gamma}_m\left(\frac{\nu +\tau }{2}\right)-\frac{m\tau}{2}\log \pi -\log {\Gamma}_m\left(\frac{\nu }{2}\right)\\ {}-\frac{m}{2}\log \mid \boldsymbol{S}\mid +\frac{\nu }{2}\log \mid \boldsymbol{\Psi} \mid -\frac{\nu +\tau }{2}\log \left|\boldsymbol{\Psi} +\left({\boldsymbol{Y}}_{\tau }-\boldsymbol{\mu} \boldsymbol{\Phi} \right){\boldsymbol{S}}^{-1}{\left({\boldsymbol{Y}}_{\tau }-\boldsymbol{\mu} \boldsymbol{\Phi} \right)}^{\top}\right|\end{array}} $$

where $ {\Gamma}_m\left(\cdot \right) $ denotes the multivariate gamma function and $ \boldsymbol{S}={\boldsymbol{I}}_{\tau \times \tau }+{\boldsymbol{\Phi}}^{\top}\boldsymbol{V}\boldsymbol{\Phi } $ .

2.2. Multi-model approach

For systems operating under time-varying conditions, the assumption of stationarity for the discrete time-series $ \boldsymbol{Y} $ is generally violated, as the statistical properties of the signal may change over time. A practical solution is to analyze the system on a local level by partitioning the full operation into operating regimes (Avendaño-Valencia, Reference Avendaño-Valencia2025), each approximately satisfying the stationarity assumption. Extending the Bayesian formulation introduced earlier (Figure 2), the VAR model can now be conditioned on the operating regime $ k\in \left[1,2,\cdots, K\right] $ , where $ K $ is the total number of operating regimes, as illustrated in Figure 3.

Figure 3.

Bayesian network of the VAR model conditioned on the operating regime $ k $ .

For each operating regime $ k $ , a separate set of prior hyperparameters $ \left\{{\boldsymbol{\mu}}_k,{\boldsymbol{V}}_k,{\boldsymbol{\Psi}}_k,{\nu}_k\right\} $ is introduced, such that the priors are conditioned on the operating regime $ k $ :

(42) $$ \boldsymbol{\Theta} \mid {\boldsymbol{\Sigma}}_w,k\sim \mathcal{M}\mathcal{N}\left({\boldsymbol{\mu}}_k,{\boldsymbol{\Sigma}}_w,{\boldsymbol{V}}_k\right),{\boldsymbol{\Sigma}}_w\mid k\sim \mathcal{I}\mathcal{W}\left({\boldsymbol{\Psi}}_k,{\nu}_k\right) $$

Given a segment belonging to a regime $ k $ , the posterior hyperparameters $ \left\{{\boldsymbol{\mu}}_k^{\left(\mathrm{post}\right)},{\boldsymbol{V}}_k^{\left(\mathrm{post}\right)},{\boldsymbol{\Psi}}_k^{\left(\mathrm{post}\right)},{\nu}_k^{\left(\mathrm{post}\right)}\right\} $ follow the updates in Equations (21)–(24) using the regime-specific priors. However, when the regime $ k $ is unknown, it must be inferred, which is detailed in the following.

2.3. Definition of local references

The regime-specific hyperparameters $ \left\{{\boldsymbol{\mu}}_k^{\left(\mathrm{post}\right)},{\boldsymbol{V}}_k^{\left(\mathrm{post}\right)},{\boldsymbol{\Psi}}_k^{\left(\mathrm{post}\right)},{\nu}_k^{\left(\mathrm{post}\right)}\right\} $ , representing the local references, can be defined either using prior knowledge of the operating regimes (supervised setting) or inferred directly from the data (unsupervised setting).

In the supervised setting, prior knowledge of the operating regimes is used to partition the data accordingly. The regime-specific hyperparameters representing the local references are then inferred from the grouped data, resulting in a set of reference models associated with each predefined regime. In the unsupervised setting, the operating regimes are identified directly from the data based on the log marginal likelihood in Equation (41). Initially, prior hyperparameters are specified, and an initial batch of data is used to infer posterior hyperparameters, defining the first local reference. As new batches of data are observed sequentially, the log marginal likelihood of each batch is evaluated with respect to the existing reference models $ k $ as:

(43) $$ {\displaystyle \begin{array}{l}\log p\left({\boldsymbol{Y}}_{\tau}^{(u)}|\boldsymbol{\Phi}, k\right)=\log {\Gamma}_m\left(\frac{\nu_k+\tau }{2}\right)-\frac{m\tau}{2}\log \pi -\log {\Gamma}_m\left(\frac{\nu_k}{2}\right)\\ {}-\frac{m}{2}\log \mid \boldsymbol{S}\mid +\frac{\nu_k}{2}\log \mid {\boldsymbol{\Psi}}_k\mid -\frac{\nu_k+\tau }{2}\log \left|{\boldsymbol{\Psi}}_k+\left({\boldsymbol{Y}}_{\tau}^{(u)}-{\boldsymbol{\mu}}_k\boldsymbol{\Phi} \right){\boldsymbol{S}}^{-1}{\left({\boldsymbol{Y}}_{\tau}^{(u)}-{\boldsymbol{\mu}}_k\boldsymbol{\Phi} \right)}^{\top}\right|\end{array}} $$

where $ \boldsymbol{S}={\boldsymbol{I}}_{\tau \times \tau }+{\boldsymbol{\Phi}}^{\top }{\boldsymbol{V}}_k\boldsymbol{\Phi} $ . Among the existing operating regimes, the one that maximizes the log marginal likelihood is selected as the most probable explanation of the new data batch:

(44) $$ \hat{c}=\arg \underset{k\in \left[1,\dots, K\right]}{\max}\log p\left({\boldsymbol{Y}}_{\tau}^{(u)}|\boldsymbol{\Phi}, k\right) $$

If the likelihood is sufficiently high, the corresponding reference model of the operating regime $ \hat{c} $ is updated using the new data. Conversely, if the likelihood is low, the data are not well explained by the existing references, and a new local reference is initialized based on the inferred posterior hyperparameters of that batch. This sequential procedure is repeated over the available data, resulting in a set of local references learned directly from the data without requiring prior knowledge of regime boundaries or the number of regimes.

2.3.1. Detection of anomalies

A large value of the log marginal likelihood given a segment $ {\boldsymbol{Y}}_t $ indicates that the new data segment is highly consistent with the statistical structure represented by regime $ \hat{c} $ . Conversely, low values indicate that the observed data are poorly explained by the corresponding regime-specific model. This may indicate that the segment $ {\boldsymbol{Y}}_t $ corresponds to an operating condition that deviates significantly from the defined operating regimes and thereby the regime-specific dynamics. However, individual segments that are poorly explained by the local references may arise due to transient deviations or isolated outliers. To improve robustness against such effects, the anomaly assessment is not performed on individual segments, but rather based on an accumulated statistic over a sliding window of length $ J $ segments. Specifically, the anomaly detection statistic is defined as:

(45) $$ {D}_t^{\left(\hat{c}\right)}=\sum \limits_{j=0}^J\log p\left({\boldsymbol{Y}}_{t-j}|\boldsymbol{\Phi}, \hat{c}\right) $$

which represents the cumulative log marginal likelihood of the most recent observations under regime $ \hat{c} $ . Low values of $ {D}_t^{\left(\hat{c}\right)} $ indicate that the recent sequence of observations is unlikely under the reference model, and therefore suggest a persistent deviation from the learned regime behavior. Based on this, the structural condition is assessed as:

(46) $$ \mathrm{Consistent}\ \mathrm{with}\ \mathrm{reference}\ \mathrm{if}:{D}_t^{\left(\hat{c}\right)}\ge {\vartheta}_{\hat{c}} $$

(47) $$ \mathrm{Anomaly}\ \mathrm{detected}\ \mathrm{if}:{D}_t^{\left(\hat{c}\right)}<{\vartheta}_{\hat{c}} $$

where $ {\vartheta}_{\hat{c}} $ denotes a regime-dependent threshold. The threshold $ {\vartheta}_c $ is determined empirically from the training data by estimating the distribution of $ {D}_t^{(c)} $ under each regime and selecting a lower-tail quantile. This provides a statistically interpretable criterion for anomaly detection, ensuring that the decision is consistent with the variability observed within each regime. This accumulated statistic can be interpreted as a measure of joint evidence over multiple observations, reducing sensitivity to isolated outliers while emphasizing persistent deviations.

As an extension to the above formulation, the evolution of the accumulated statistic $ {D}_t^{\left(\hat{c}\right)} $ over time can be monitored to identify gradual deviations from the reference condition. A consistent downward trend in $ {D}_t^{\left(\hat{c}\right)} $ indicates a progressive reduction in the log marginal likelihood, which may be interpreted as a gradual degradation of the system. This enables the framework to provide early-warning indicators of structural changes beyond binary anomaly detection.

3.1. Data acquisition

The CTV is equipped with four propulsion units, each consisting of a diesel engine, a gearbox, and a propeller set, as shown in Figure 5, where the propellers are facing forward. As the CTV is a catamaran, each hull contains two propulsion units. Accordingly, the data acquisition system measures the vibrational response of a single gearbox located in one of the hulls. The vibrational response of the gearbox is measured by a three-axis accelerometer mounted on top of the gearbox Figure 5. The accelerometer is oriented such that the x-axis corresponds to forward–backwards motion, the y-axis to lateral motion, and the z-axis to vertical motion. In this context, the gearbox is considered a coupled dynamical system, where the measured vibrational response arises from the interaction of multiple internal components, including gears, shafts, and bearings. The three-axis accelerometer therefore captures different projections of this coupled behavior, resulting in a multivariate signal.

Figure 5.

Location of the accelerometer on one of the propulsion units.

For this preliminary study, where the proposed multi-model framework is tested, data collected over 1 month during summer, in which the vibrational response was measured, is used to establish the pool of references. During this month, 14 days of operation were cancelled, resulting in 17 days of operation where the vibrational response was measured. Subsequently, a small dataset including 4 days of operation is used to test the multi-model on a dataset without any structural changes.

Approximately 6 months later, repairs were carried out on both drive shafts between the engines and gearboxes, and new sets of propellers were installed. A small dataset collected after these repairs is therefore used to test the detectability of anomalies, although data capturing gradual structural degradation would have been preferable. However, the limited size of the post-repair dataset prevents it from being used as a reference condition. Furthermore, the available post-repair dataset represents a step change in the system condition rather than gradual degradation. While this is not fully representative of typical damage evolution, it provides a controlled change in system behavior that enables evaluation of the multi-model’s ability to detect deviations from the learned reference condition. The complete dataset is available through Zenodo (Drangsfeldt, Reference Drangsfeldt2026). The details of the data acquisition system and the dimensions of the datasets are detailed in Table 1.

Table 1.

Overview of the data acquisition

Table 1. Long description

The table consists of two columns: Description and Value.

* Row 1: Accelerometer is M M F K S 9 4 3 B 1 0 0.

* Row 2: D A Q is D E W E dash 4 3 A.

* Row 3: Software is D E W E Soft X.

* Row 4: Data dimension to establish references is 3 times 2,697,650,000.

* Row 5: Data dimension for test without structural changes is 3 times 486,733,301.

* Row 6: Data dimension to test detectability is 3 times 499,266,354.

* Row 7: Sampling frequency is 5000 H z.

3.2. Comparison to a traditional SHM approach

As a benchmark representing a traditional SHM approach, one of the methodologies presented in Drangsfeldt and Avendaño-Valencia (Reference Drangsfeldt and Avendaño-Valencia2024) is adopted. This approach is based on using the parameters of VAR models as DSFs, implicit mitigation of EOV through PCA, and the squared Mahalanobis distance as a measure of correspondence to a global reference model. This benchmark serves as a reference for evaluating the limitations of global reference-based approaches under discrete and varying operational conditions, where regime-dependent variability is not explicitly accounted for.

First, the dataset used to establish the reference model is segmented into intervals of 5-second duration, resulting in $ 107906 $ segments, each represented as $ \boldsymbol{Y}\in {\unicode{x211D}}^{3\times 25000} $ . A VAR model is fitted to each segment using Least Squares estimation as defined in Equations (5) and (7). The estimated AR parameters, $ \boldsymbol{\Theta} $ , constitute the DSFs. The complete set of DSFs from all $ 107906 $ segments defines the global reference:

(48) $$ \mathcal{M}=\left[\begin{array}{cccc}{\boldsymbol{\Theta}}_1& {\boldsymbol{\Theta}}_2& \cdots & {\boldsymbol{\Theta}}_{107906}\end{array}\right] $$

where $ \mathcal{M} $ represents the global reference state of the gearbox. Following the PCA procedure for implicit mitigation of EOV, each $ \boldsymbol{\Theta} $ is first vectorized, denoted as $ \overline{\boldsymbol{\theta}}=\mathrm{vec}\left(\boldsymbol{\Theta} \right) $ , and subsequently normalized to obtain $ \tilde{\boldsymbol{\theta}} $ , using the mean and standard deviation estimated from the vectorized reference $ \overline{\mathcal{M}} $ . A Singular Value Decomposition is then applied to the vectorized, normalized reference:

(49) $$ \tilde{\mathcal{M}}=\left[\begin{array}{cccc}{\tilde{\boldsymbol{\theta}}}_1& {\tilde{\boldsymbol{\theta}}}_2& \cdots & {\tilde{\boldsymbol{\theta}}}_{107906}\end{array}\right]=\boldsymbol{U}\hskip0.1em \boldsymbol{S}\hskip0.1em {\boldsymbol{P}}^{\top } $$

where $ \boldsymbol{P} $ denotes the right singular vector matrix. This is used to transform each incoming, vectorized, and normalized AR parameter vector $ \tilde{\boldsymbol{\theta}} $ into Principal Component (PC) space:

(50) $$ \boldsymbol{Z}=\tilde{\boldsymbol{\theta}}\hskip0.1em \boldsymbol{P} $$

As proposed in Drangsfeldt and Avendaño-Valencia (Reference Drangsfeldt and Avendaño-Valencia2024), the PCs accounting for up to 90% of the accumulated variance are discarded, as they are assumed to primarily represent EOV. Similarly, PCs beyond 99% of the accumulated variance are removed to reduce the influence of noise. Finally, the squared Mahalanobis distance is used to assess whether an incoming observation corresponds to an outlier relative to the transformed global reference:

(51) $$ {d}^2\left(\boldsymbol{Z},{\mathcal{M}}^{\star}\right) $$

where $ {\mathcal{M}}^{\star } $ denotes the global reference projected into the reduced PC space. Then, by treating every 5-second segment as an incoming observation, the squared Mahalanobis distance is shown in Figure 6. Under the assumption that the retained PCs are approximately Gaussian, the squared Mahalanobis distance follows a chi-square distribution, $ {d}^2\sim {\chi}_r^2 $ , with degrees of freedom equal to the number of retained PCs. Therefore, outlier (anomaly) detection is performed as:

(52) $$ {d}^2\left(\boldsymbol{Z},{\mathcal{M}}^{\star}\right)\ge {\chi}_{r,1-\alpha}^2\to {H}_1:\mathrm{Anomaly}\ \mathrm{is}\ \mathrm{detected} $$

(53) $$ \mathrm{Otherwise}\to {H}_0:\mathrm{No}\;\mathrm{anomaly}\ \mathrm{is}\ \mathrm{detected} $$

where $ r $ denotes the number of retained PCs, and $ {\chi}_{r,1-\alpha}^2 $ is the $ \left(1-\alpha \right) $ -quantile of the chi-square distribution. The horizontal threshold shown in Figure 6 represents a confidence level of 99%, corresponding to $ \alpha =0.01 $ .

Figure 6.

The squared Mahalanobis distance for each 5-second segment follows the benchmark approach for anomaly detection.

Figure 6. Long description

A scatter plot with a logarithmic Y axis labeled Squared Mahalanobis distance d squared (Z, M super asterisk) ranging from 10 super 1 to 10 super 4. The X axis is labeled 5-second segments ranging from 0 to 150,000. A horizontal dashed black line represents a threshold at approximately 10 super 2.

* The first phase, References, is shaded dark gray and spans from segment 0 to approximately 107,000. It contains blue data points mostly clustered below the threshold, with several vertical spikes reaching 10 super 3 and one extreme outlier near 10 super 4 at segment 62,000.

* The second phase, Validation, is shaded light gray from segment 107,000 to 127,000. It continues with blue data points showing similar density and variance to the reference phase.

* The third phase, After repairs, is shaded light orange from segment 127,000 to the end. The data points in this section are colored orange. The variance appears slightly reduced compared to previous phases, with most points tightly clustered around the threshold line and fewer high-magnitude outliers.

As a result, 46.8% of the 5-second segments from the dataset acquired after the repairs are classified as anomalies. However, 12.6% of the segments from the dataset prior to the repairs are also classified as anomalies, indicating a relatively high false positive rate for the global reference-based approach adopted from Drangsfeldt and Avendaño-Valencia (Reference Drangsfeldt and Avendaño-Valencia2024).

3.3. Establishing the references of the multi-model

The first step of the proposed multi-model framework in the supervised setting is to establish an accurate representation of the posterior hyperparameters for each regime. As detailed in Section 2.1.1, this can be achieved either through a batch update or a sequential update process. A batch update, however, presents two main challenges. First, the matrix inversion required in Equation (21) becomes computationally expensive when the batch is large, since the regression matrix $ \boldsymbol{\Phi} $ grows proportionally with the data length. Second, maintaining the stationarity assumption for $ \boldsymbol{Y} $ becomes increasingly difficult when a long data segment is used. In contrast, a sequential update is less computationally demanding because it processes smaller segments iteratively. However, the iterative procedure can be extensive, and some segments may continue to violate the stationarity assumption. To address the challenge of possible non-stationarities in the data obtained, a regularized updating process is adopted, as detailed in Section 2.1.1. This approach introduces forgetting factors to control the influence of newly observed data relative to previously inferred information, thereby mitigating the effect of short data segments that may violate the stationarity assumption.

First, the dataset used to establish the local references is segmented into segments of 5 seconds. This results in $ \boldsymbol{Y}\in {\unicode{x211D}}^{3\times 25000} $ , where $ \boldsymbol{Y} $ is associated with its corresponding operational mode according to the DPR. Then, flat priors for each operational mode are adopted initially, and then updated sequentially using the segmented $ \boldsymbol{Y} $ . These prior hyperparameters are therefore defined such that the corresponding prior distributions are flat:

(54) $$ {\boldsymbol{\Psi}}_k={\boldsymbol{I}}_{m\times m}\&{\nu}_k=d $$

(55) $$ {\boldsymbol{\mu}}_k={\mathbf{0}}_{m\times d}\&{\boldsymbol{V}}_k=1000\cdot {\boldsymbol{I}}_{d\times d} $$

The forgetting factors for the regularized updating process are selected as 0.01 and 0.05, respectively, reflecting a balance between responsiveness to evolving system dynamics and robustness to transient deviations and noise. These values were chosen to ensure stable convergence of the inferred hyperparameters while maintaining sufficient adaptability to changes in the observed behavior. The first part of the updating strategy is shown in Figure 7, displaying the first four coefficients in $ \boldsymbol{\Theta} $ for operational mode 3. The underlying grey Credible Interval (CI) and mean labelled as No updating correspond to the case where every 5-second segment in $ \boldsymbol{Y} $ follows the batch updating strategy where no updating is performed. In this case, the prior hyperparameter $ {\boldsymbol{\mu}}_k $ is estimated using Least Squares estimation and the prior hyperparameter $ {\boldsymbol{V}}_k $ is set to $ 0.001\cdot {\boldsymbol{I}}_{d\times d} $ . Additionally, prior hyperparameters $ {\boldsymbol{\Psi}}_k $ and $ {\nu}_k $ are set as in Equation (54). With the defined forgetting factors, the posterior gradually adapts to new observations while avoiding overfitting to minor variations observed when every 5-second segment follows the batch updating strategy.

Figure 7.

Example of the updating strategy for Mode 3 after having observed 144 sets of 5 seconds duration, all belonging to Mode 3.

Figure 7. Long description

A legend at the top identifies four data series: No updating C I (gray shaded area), No updating mean (black dotted line), Updated 95% C I (blue shaded area), and Updated mean (solid blue line). All four panels share an x-axis labeled No. of 5 s sets ranging from 0 to 144 and a y-axis labeled Value.

* Panel 1 (theta sub 1): The updated mean shows a steep exponential decay from 0 to a steady state of negative 0.5. The non-updated mean fluctuates around negative 0.55 with a wider gray C I.

* Panel 2 (theta sub 2): The updated mean shows a gradual linear decrease from 0 to negative 0.07. The non-updated mean shows high-frequency oscillations between 0 and negative 0.1.

* Panel 3 (theta sub 3): The updated mean shows an asymptotic increase from 0, leveling off at 0.32. The non-updated mean is relatively flat, fluctuating around 0.32 with significant noise.

* Panel 4 (theta sub 4): The updated mean shows an asymptotic increase from 0, leveling off at 0.42. The non-updated mean fluctuates around 0.45 with a gray C I that narrows slightly over time.

Lastly, the model order $ {n}_a $ is selected based on a systematic model selection procedure. Specifically, Least Squares estimations are performed on multiple data segments representing different operational modes, and the resulting models are evaluated using both the Residual Sum of Squares over the Series Sum of Squares (RSS/SSS) and the Bayesian Information Criterion (BIC). These criteria provide a balance between model fit and model complexity. Based on this analysis, a model order of $ {n}_a=20 $ was found to provide a suitable trade-off and is therefore adopted for all operational modes.

3.4. Procedure of the multi-model

After establishing the posterior hyperparameters for each regime, which serve as local references, the procedure follows Equation (44) for regime predictions among the existing reference models, where every 5-second segment is treated as a new observation. Following Equation (44), the log marginal likelihood of observing every 5-second segment given each reference operational mode is calculated, and determines the classified operational modes. The maximum log marginal likelihood for every 5-second segment is shown in Figure 8 with its classified operational mode. Here, every 5-second segment in both the dataset used for establishing the references and for testing without structural changes is treated as a new observation. It is observed that, despite significant variability for some of the operational modes, the 5-second segments in the testing dataset align well with the 5-second segments in the dataset used to define the references.

Figure 8.

Maximum log marginal likelihood of observing every 5-second segment with the associated predicted operational mode, $ p\left({\boldsymbol{Y}}_{\tau}^{(u)}|\boldsymbol{\Phi}, \hat{c}\right) $ , where every segment in both the reference and test dataset is treated as a new observation.

Figure 8. Long description

A multi-panel figure with seven horizontal scatter plots. A legend at the top indicates dark gray shading for References and light gray shading for Validation. Each plot shares a common y-axis label: log p ( Y sub tau super (u) | Phi, c-hat) times 10 super 5, and a common x-axis label: 5 s segments.

* Mode 1: Data points fluctuate between negative 3.4 and negative 3.0. The reference region ends near segment 4700.

* Mode 2: Data shows a dense band between negative 3.5 and negative 3.0, extending to approximately 6500 segments.

* Mode 3: Likelihood values are concentrated between negative 3.0 and negative 2.5, with significant downward spikes. The reference region ends at segment 45000.

* Mode 4: A very stable horizontal trend near negative 4.0 with occasional sharp downward drops. The reference region ends at segment 27000.

* Mode 5: Data fluctuates around negative 3.7, with the reference region concluding at segment 4700.

* Mode 6: Values are mostly stable near negative 3.5 with downward noise. The reference region ends at segment 15000.

* Mode 7: Shows a distinct shift in likelihood levels around segment 1400, where the reference region ends and validation begins, with values rising from negative 3.0 toward negative 2.8.

Finally, the maximum log marginal likelihood for every segment, as shown in Figure 8, also serves as an anomaly indicator as detailed in Section 2.3.1. This is elaborated and analyzed in the following section.