# Online Student-$t$ Processes with an Overall-local Scale Structure for
Modelling Non-stationary Data

Taole Sha  Michael Minyi Zhang

###### Abstract

Time-dependent data often exhibit characteristics, such as non-stationarity
and heavy-tailed errors, that would be inappropriate to model with the typical
assumptions used in popular models. Thus, more flexible approaches are
required to be able to accommodate such issues. To this end, we propose a
Bayesian mixture of student-$t$ processes with an overall-local scale
structure for the covariance. Moreover, we use a sequential Monte Carlo (SMC)
sampler in order to perform online inference as data arrive in real-time. We
demonstrate the superiority of our proposed approach compared to typical
Gaussian process-based models on real-world data sets in order to prove the
necessity of using mixtures of student-$t$ processes.

##  1 Introduction

In modelling dynamical systems, it is common that the data will exhibit non-
stationarity, where the trend changes across the input space. Kernel methods
like the Gaussian process (GP) are a popular choice of prior distribution over
real-valued functions in Bayesian models of time series data (Rasmussen and
Williams,, 2005). However, in the non-stationary time series setting that this
paper focuses on, they face several challenges: 1.) The calculation of the
likelihood in GP inference requires inverting an $N\times N$ matrix, which
generally incurs a computational complexity of $\mathcal{O}(N^{3})$ where $N$
is the number of observations; 2.) Updating the model in real-time is not
trivial; 3.) Stationarity is often assumed by covariance kernels, while non-
stationary kernels typically lead to computationally intractable GPs,
especially when the sample size is large. As a related stochastic process, the
student-$t$ process (TP) has been treated as an attractive alternative prior
distribution over function space compared to the GP (Shah et al.,, 2014). The
TP is a more general elliptical process, where the likelihood of observations
decreases in their distance from the mode, which is a reasonable assumption
for the prior. It also has heavy tails controlled by the degree of freedom
parameter, allowing more modelling flexibility.

Similar to the GP, the TP has consistent marginals and closed-form
conditionals which make it as convenient as the GP to use in statistical
modelling without any additional computational cost. However, TPs are still
liable to suffer from the aforementioned three issues that GPs face when
modelling real-world data. Hence, we introduce in this paper a mixture of TPs
with an SMC sampler, so that we may take advantage of the additional
flexibility of a mixture-of-experts model with a convenient online inference
algorithm. To derive the TP, we assume a latent GP and integrate out an
inverse gamma prior on the kernel amplitude and the noise parameters.
Moreover, we add an additional parameter of the noise term to control the
heteroscedasticity. Lastly, we model the level of heavy-tailedness by
automatically controlling the TPs’ degree of freedom using an efficient slice
sampling scheme.

Our paper proceeds as follows: In Section 2, we discuss some previous work
about online GP models. We introduce our way of handling noisy data using TP
in Section 3. The online TP inference algorithm is detailed in Section 4. We
use the experiment results to compare it with GP-based models in Section 5.
Finally, we conclude the paper in Section 6 with a discussion of future work.

##  2 Related Work

The Gaussian process is a typical choice of prior used in Bayesian methods for
modelling time series and non-linear regression problems (Rasmussen and
Williams,, 2005). A GP distributed function,
$f\sim\mathcal{GP}(\mu(\cdot),\Sigma(\cdot,\cdot)$, is defined by a mean
function, $\mu(\cdot)$, and a covariance function, $\Sigma(\cdot,\cdot)$, with
a property that GPs are multivariate normally distributed conditioned on a
finite set of points: $f(x)\sim\mathcal{N}(\mu(x),\Sigma(x,x^{\prime}))$.
While the GP is a convenient choice of prior due to this multivariate normal
property, as it leads to tractable posterior inference in many classes of
models, GPs suffer from the typical cubic computational that other kernel
methods face.

Numerous scalable methods have been developed to tackle the computational
issue of GPs: Sparse inducing point methods are a popular technique for
reducing the computational complexity of GP methods (Snelson and Ghahramani,,
2006; Titsias,, 2009; Bauer et al.,, 2016). In the sparse GP methods, they
form a low-rank approximation of the kernel function using a collection of $M$
“pseudo-inputs” which reduce the computational complexity of the GP to
$O(NM^{2})$ from $O(N^{3})$. Product-of-expert models employ a block diagonal
approximation of the full covariance matrix in order to reduce the complexity
of the full covariance matrix inversion to individually inverting each smaller
block (Deisenroth and Ng,, 2015; Cohen et al.,, 2020). While not necessarily
faster, mixture-of-expert models use a mixture of GPs to model functions with
greater flexibility compared to a single GP (Rasmussen and Ghahramani,, 2001;
Meeds and Osindero,, 2005).

For fast online GP methods, Csató and Opper, (2002) used variational inference
to approximate the posterior in a sparse online GP model, however, the
hyperparameters are assumed to be fixed in their method. Nguyen-tuong et al.,
(2008) proposed a product-of-experts local GP method for online fitting, where
the weights are based on the distance of the new observation to the local
models. Though, in such methods, ignoring the correlation between experts when
adopting the local assumption can lead to poor uncertainty quantification. Bui
et al., (2017) developed a sparse variational GP regression approach that
allows for online updating of the hyperparameters, called OSVGP. However,
OSVGP has a tendency to be numerically unstable and, empirically, is liable to
underfit the data. Stanton et al., (2021) developed an exact sparse online
model called WISKI, where a structured and sparse covariance matrix
approximation developed by Wilson and Nickisch, (2015) is used, leading to
constant computational complexity with respect to the number of observations.

Regarding SMC methods in GPs, Svensson et al., (2015) proposed an SMC sampler
with the purpose of marginalizing the kernel hyperparameters and Gramacy and
Polson, (2011) proposed an SMC sampler for sequential design in GPs. While
these SMC methods allow for updating the GP model sequentially, they cannot
account for non-stationarity in the data, nor are they able to limit the
computational cost of the model as the complexity still scales $O(N^{3})$.
However, Zhang and Williamson, (2019) proposed an importance sampling method
for scaling up a mixture-of-experts GP model to an average complexity of
$O(N^{3}/K^{2})$ for non-stationary data. Later, Zhang et al., (2023) and
Härkönen et al., (2022) developed an online SMC and SMC2 sampler for mixture
of GPs. But despite the advances in online mixtures of GPs, little attention
has been paid to online mixtures of the student-$t$ process.

##  3 Student-$t$ Process for Noisy Data

Student-$t$ priors have long been used in Bayesian linear regression problems,
where we may desire modelling sparse regression coefficients or heavy-tailed
errors Fernández and Steel, (1999); Tipping, (2001); West, (1984); Geweke,
(1993). Vanhatalo et al., (2009) introduced a more robust method of GP
regression where the latent function was GP distributed but the observation
likelihood was assumed to be a student-$t$ distribution. However, they
estimated all the parameters in this model using a Laplace approximation to
the posterior distribution instead of performing exact Bayesian inference and
therefore can properly capture the posterior uncertainty. Later, Jylänki et
al., (2011) used an expectation propagation algorithm for posterior inference
in the same model. Again, expectation propagation is only an approximate
method for posterior inference that cannot exactly capture the underlying
uncertainty.

When modelling the noisy time series data with $D$-dimensional input
$\mathbf{x}_{i}\in\mathbb{R}^{D}$ and output $y_{i}\in\mathbb{R}$, where times
$i=1,2,\ldots,N$, student-$t$ processes (TPs) are an attractive alternative to
the Gaussian processes (Shah et al.,, 2014). Since student-$t$ distributions
are not closed under addition, we cannot analytically obtain a latent TP and
independent student-t noise for modelling regression problems. As a result,
this model formulation is not convenient for the statistical practitioner
compared to the GP. Instead, we may obtain a TP by incorporating the noise
into the kernel function. According to Shah et al., (2014), Zhang and Yeung,
(2010) wrongly assumed the noise to be independent when raising this model.
Shah et al., (2014) stated that the noise term is uncorrelated but dependent,
and argued it to behave similarly to a sum of a latent TP with independent
noise. Tang et al., (2017) combined both a student-$t$ process model with
student-$t$ noise, but, again, used only a Laplace approximation for the
posterior instead of performing exact inference.

However, directly incorporating the noise term into the kernel may not be
sufficiently flexible for modelling real data. In our model, we handle noisy
data using an additional heteroscedastic parameter for each mixture. It is
assumed that $y_{i}=f(\mathbf{x}_{i})+\sigma_{0}\epsilon_{i}$, the output is
generated by a latent zero-mean Gaussian Process $f(\mathbf{x}_{i})$ and a
Gaussian noise term where $\epsilon_{i}\sim N(0,\sigma_{1}^{2})$.
$\sigma_{0}^{2}$ is an overall scale parameter for both the covariance
function $\sigma_{0}^{2}\Sigma(\cdot,\cdot)$ and the noise term, while
$\sigma_{1}^{2}$ is a scale parameter for the noise term to control the
heteroscedasticity. The GP and the noise are not independent here since they
share the same overall scale. Any finite samples drawn from the GP at
locations $\mathbf{X}=[\mathbf{x}_{1},\cdots,\mathbf{x}_{N}]^{T}$ jointly
follow a multivariate normal distribution:

$\displaystyle
f(\mathbf{X})|\mathbf{X},\sigma_{0}^{2}\sim\mathcal{N}_{N}(0,\sigma_{0}^{2}\Sigma_{\mathbf{X}\mathbf{X}^{\prime}}).$
(1)

For noisy observations $\mathbf{y}=[y_{1},\cdots,y_{N}]^{T}$, the data is
generated by:

$\displaystyle\mathbf{y}|\mathbf{X},f,\sigma_{0}^{2}\sim\mathcal{N}_{N}(f(\mathbf{X}),\sigma_{0}^{2}\sigma_{1}^{2}\mathbf{I}).$
(2)

Due to the conjugacy between Gaussians, marginally:

$\displaystyle\mathbf{y}|\mathbf{X},\sigma_{0}^{2}\sim\mathcal{N}_{N}(0,\sigma_{0}^{2}(\Sigma_{\mathbf{X}\mathbf{X}^{\prime}}+\sigma_{1}^{2}\mathbf{I})).$
(3)

By marginalizing an inverse Gamma prior on $\sigma_{0}^{2}$ out, we can also
arrive at the target multivariate student-t distribution with degree of
freedom $\nu$:

$\displaystyle\sigma_{0}^{2}\sim\mbox{Inv\text{-}Gamma}\left(\frac{\nu}{2},\frac{\nu}{2}\right),$
$\displaystyle\mathbf{y}|\mathbf{X}\sim\mathcal{T}_{N}(\nu,0,\Sigma_{\mathbf{X}\mathbf{X}^{\prime}}+\sigma_{1}^{2}\mathbf{I}).$
(4)

The TP’s log marginal likelihood is:

$\displaystyle\log
P(\mathbf{y}|\nu,\Sigma_{\mathbf{X}\mathbf{X}^{\prime}},\sigma_{1}^{2})=-\frac{N}{2}\log(\nu\pi)$
$\displaystyle~{}~{}~{}~{}~{}-\frac{1}{2}\log(|\Sigma_{\mathbf{X}\mathbf{X}^{\prime}}+\sigma_{1}^{2}\mathbf{I}|)+\log\left(\frac{\Gamma(\frac{\nu+N}{2})}{\Gamma(\frac{\nu}{2})}\right)$
$\displaystyle~{}~{}~{}~{}~{}-\frac{\nu+N}{2}\log\left(1+\frac{\mathbf{y}^{T}(\Sigma_{\mathbf{X}\mathbf{X}^{\prime}}+\sigma_{1}^{2}\mathbf{I})^{-1}\mathbf{y}}{\nu}\right).$
(5)

When making predictions, the posterior predictive distribution of $N^{*}$
target outputs $\mathbf{y}^{*}$ given new inputs $\mathbf{X}^{*}$ is:

$\displaystyle\tilde{\phi}_{2}=\Sigma_{\mathbf{X}^{*}\mathbf{X}}(\Sigma_{\mathbf{X}\mathbf{X}^{\prime}}+\sigma_{1}^{2}\mathbf{I})^{-1}\mathbf{y},$
$\displaystyle\beta_{1}=\mathbf{y}^{T}(\Sigma_{\mathbf{X}\mathbf{X}^{\prime}}+\sigma_{1}^{2}\mathbf{I})^{-1}\mathbf{y},$
$\displaystyle\tilde{\mathbf{K}}_{22}=(\Sigma_{\mathbf{X}^{*}\mathbf{X}^{{*}^{\prime}}}+\sigma_{1}^{2}\mathbf{I})$
$\displaystyle~{}~{}~{}~{}~{}-\Sigma_{\mathbf{X}^{*}\mathbf{X}^{\prime}}(\Sigma_{\mathbf{X}\mathbf{X}^{\prime}}+\sigma_{1}^{2}\mathbf{I})^{-1}\Sigma_{\mathbf{X}\mathbf{X}^{*}},$
$\displaystyle\mathbf{y}^{*}|\mathbf{y},\mathbf{X},\mathbf{X}^{*}\sim\mathcal{T}(\nu+N,\tilde{\phi}_{2},\frac{\nu+\beta_{1}}{\nu+N}\tilde{\mathbf{K}}_{22}).$
(6)

##  4 Online Student-$t$ Processes for Non-stationary Data

The data generating process for our proposed model is:

$\displaystyle\textbf{x}_{i}\sim\mathcal{T}(\boldsymbol{\mu}_{z_{i}},\boldsymbol{\Psi}_{z_{i}},\nu_{z_{i}}),$
$\displaystyle\alpha\sim\mbox{Gamma}(a_{0},b_{0}),~{}~{}z_{i}|\alpha\sim\mbox{CRP}(\alpha),$
$\displaystyle\boldsymbol{\theta}_{k}\sim
log\mathcal{N}(m_{0},s_{0}^{2}\mathbf{I}),~{}~{}\nu_{k}\sim\mbox{Gamma}(2,0.1),$
$\displaystyle h_{k}\sim\mathcal{N}(0,k_{0}^{2}),~{}~{}k_{0}^{2}\sim\mbox{Inv-
Gamma}\left(\frac{1}{2},\frac{1}{2}\right),$
$\displaystyle\mathbf{y}_{k}|\mathbf{X}_{k},\boldsymbol{\theta}_{k},\sim\mathcal{T}(\nu_{k},0,\mathbf{K}_{\boldsymbol{\theta}_{k}}+|h_{k}|\mathbf{I}).$
(7)

where the $i$-th input $\mathbf{x}_{i}$ comes from an infinite Dirichlet
process Gaussian-inverse Wishart mixture model (Antoniak,, 1974):
$\mathcal{N}(\mathbf{M}_{z_{i}},\mathbf{C}_{z_{i}})$. The latent parameters
$(\mathbf{M}_{z_{i}},\mathbf{C}_{z_{i}})$, for $z_{i}\in\\{1,2,\ldots\\}$, are
integrated out over a normal-inverse Wishart prior,
$\mathcal{NIW}(\boldsymbol{\mu}_{z_{i}},\lambda_{z_{i}},\boldsymbol{\Psi}_{z_{i}},\nu_{z_{i}})$.
The $x_{i}$ marginally follows a student-$t$ distribution
$\mathcal{T}(\boldsymbol{\mu}_{z_{i}},\boldsymbol{\Psi}_{z_{i}},\nu_{z_{i}})$.

The outputs $y_{i}$ from cluster $k$ are denoted as $\mathbf{y}_{k}$, which we
obtain:

$\displaystyle\sigma_{k}^{2}|\nu_{k}\sim\mbox{Inv-
Gamma}\left(\frac{\nu_{k}}{2},\frac{\nu_{k}}{2}\right)$
$\displaystyle\mathbf{y}_{k}|\sigma_{k}^{2},\mathbf{X}_{k}\sim\mathcal{N}(0,\sigma_{k}^{2}(\mathbf{K}_{\boldsymbol{\theta}_{k}}+|h_{k}|\mathbf{I})),$
$\displaystyle\int
P(\mathbf{y}_{k}|\sigma_{k}^{2},-)P(\sigma_{k}^{2})d\sigma_{k}^{2}\sim
T(\nu_{k},0,\mathbf{K}_{\boldsymbol{\theta}_{k}}+|h_{k}|\mathbf{I}).$ (8)

It is assumed to be a sum of a GP and dependent Gaussian noise as mentioned in
Section 3. The covariance is determined by an overall scale parameter
$\sigma_{k}^{2}$ for both the kernel and the noise, a local scale parameter
$h_{k}$ for the noise only to control this heteroscedasticity, and kernel
parameters $\boldsymbol{\theta}_{k}$. The overall scale parameter
$\sigma_{k}^{2}$ is integrated out over an inverse gamma prior, and a TP can
be derived. When the $i^{th}$ streaming data $(\mathbf{x}_{i},y_{i})$ comes,
we assign it to cluster $k$ according to the predictive distribution of the
DP, the Chinese restaurant process (Aldous,, 1985):

$\displaystyle
P(z_{i}=k|\alpha,\mathbf{X}_{k})\propto\left\\{\begin{array}[]{ll}N_{k}^{\prime}\cdot\mathcal{T}(\boldsymbol{\mu}_{k}^{\prime},\boldsymbol{\Psi}_{k}^{\prime},\nu_{k}^{\prime})&k\in
K^{+}.\\\
\alpha\cdot\mathcal{T}(\boldsymbol{\mu}_{0},\boldsymbol{\Psi}_{0},\nu_{0})&\text{o.w.}\end{array}\right.$
(11)

$K^{+}$ refers to the existing clusters, and all $(\cdot)^{\prime}$ represent
summary statistics calculated with first $i-1$ observations. The student-t
likelihood’s parameters
$(\mathbf{\mu}_{k}^{{}^{\prime}},\mathbf{\Psi}_{k}^{{}^{\prime}},\nu_{k}^{{}^{\prime}})$
of inputs can be updated by:

$\displaystyle\boldsymbol{\mu}_{k}^{\prime}=$
$\displaystyle\frac{\lambda_{0}\boldsymbol{\mu}_{0}+N_{k}^{\prime}\bar{\mathbf{x}}_{k}}{\lambda_{k}^{\prime}},\bar{\mathbf{x}}_{k}^{\prime}=\frac{\sum_{i^{\prime}:(z_{i^{\prime}}=k,i^{\prime}<i)}\mathbf{x}_{i^{\prime}}}{N_{k}^{\prime}},$
$\displaystyle N_{k}^{\prime}=$
$\displaystyle\sum_{i^{\prime}=1}^{i-1}I(z_{i^{\prime}}=k),\lambda_{k}^{\prime}=\lambda_{0}+N_{k}^{\prime},$
$\displaystyle\nu_{k}^{\prime}=$ $\displaystyle\nu_{0}+N_{k}^{\prime}-D+1,$
$\displaystyle\boldsymbol{\Psi}_{k}^{\prime}=$
$\displaystyle\frac{\lambda_{k}^{\prime}+1}{\lambda_{k}^{\prime}\nu_{k}^{\prime}}\left(\boldsymbol{\Psi}_{0}+\mathbf{S}_{k}^{\prime}+\mathbf{S}_{\bar{\mathbf{x}}_{k}}^{\prime}\right)$
$\displaystyle\mathbf{S}_{k}^{\prime}=$
$\displaystyle\sum_{i^{\prime}:(z_{i^{\prime}}=k,i^{\prime}<i)}\left(\mathbf{x}_{i^{\prime}}-\bar{\mathbf{x}}^{\prime}_{k}\right)\left(\mathbf{x}_{i^{\prime}}-\bar{\mathbf{x}}^{\prime}_{k}\right)^{T}$
$\displaystyle\mathbf{S}_{\bar{\mathbf{x}}_{k}}^{\prime}=$
$\displaystyle\frac{\lambda_{0}N_{k}^{\prime}}{\lambda_{k}^{\prime}}\left(\bar{\mathbf{x}}^{\prime}_{k}-\boldsymbol{\mu}_{0}\right)\left(\bar{\mathbf{x}}^{\prime}_{k}-\boldsymbol{\mu}_{0}\right)^{T}.$
(12)

Also, a Gamma prior is placed on the Dirichlet process concentration parameter
$\alpha$. We can use a variable augmentation scheme to sample its full
conditional posterior up to observation $i$ (Escobar and West,, 1995).

$\displaystyle\rho|\alpha\sim$
$\displaystyle\text{Beta}(\alpha+1,i),K=|\\{k:N_{k}>0\\}|$
$\displaystyle\frac{\pi_{\alpha}}{1-\pi_{\alpha}}=$
$\displaystyle\frac{a_{0}+K-1}{N(b_{0}-\log\rho)}$
$\displaystyle\alpha|\mathbf{z}_{1:i},\pi_{\alpha},\rho=$
$\displaystyle(1-\pi_{\alpha})$ $\displaystyle\cdot$
$\displaystyle\text{Gamma}(\alpha_{0}+K-1,b_{0}-\log\rho)$ $\displaystyle+$
$\displaystyle\pi_{\alpha}\cdot\text{Gamma}(\alpha_{0}+K,b_{0}-\log\rho).$
(13)

The $\mbox{Gamma}(2,0.1)$ prior is commonly used when inferring the degree of
freedom, which puts mass on a large range of reasonable values for the degrees
of freedom (Juárez and Steel,, 2010). We sample the degrees of freedom
parameter through an efficient variable augmentation scheme. Given the latent
overall scale $\sigma_{k}^{2}$, the degree of freedom $\nu_{k}$ will be
independent of all other parameters and data. Due to the conjugacy between the
Gaussian likelihood and the inverse Gamma prior, we can directly Gibbs sample
the $\sigma_{k}^{2}$ from its full conditional. Then, conditioned on
$\sigma^{2}_{k}$, we sample $\nu_{k}$ using the slice sampler from
$P(\nu_{k}|\sigma_{k}^{2})$ (Neal,, 2003; Damien et al.,, 1999).

$\displaystyle\sigma_{k}^{2}|\mathbf{X}_{k},\mathbf{y}_{k},\nu_{k}\sim\mbox{Inv-
Gamma}(\alpha^{\prime}_{\sigma^{2}},\beta^{\prime}_{\sigma^{2}}),$
$\displaystyle\alpha^{\prime}_{\sigma^{2}}=\frac{\nu_{k}+N_{k}^{\prime}}{2}$
$\displaystyle\beta^{\prime}_{\sigma^{2}}=\frac{\nu_{k}+\mathbf{y}_{k}^{T}(\mathbf{K}_{\boldsymbol{\theta}_{k}}+|h_{k}|\mathbf{I})^{-1}\mathbf{y}_{k}}{2}$
(14)

We assume a hierarchical structure on the local heteroscedasticity parameter,
$|h_{k}|$, where global scale $k_{0}^{2}$ is shared over all mixtures. Here,
we will share scale data from other clusters to inform the posterior sampling
of $h_{k}$. Because $h_{k}$ has a normal prior, we can again sample the full
conditional of $k_{0}^{2}$ in closed form:

$\displaystyle k_{0}^{2}|h_{1},\ldots,h_{K}\sim\mbox{Inv-
Gamma}\left(\frac{K+1}{2},\frac{1+\sum_{i=1}^{K}h_{i}^{2}}{2K}\right).$ (15)

Then we sample $h_{k}$ and the TP parameters $\boldsymbol{\theta}_{k}$ using
the elliptical slice sampler (ESS), which is an efficient sampling algorithm
for non-conjugate models with Gaussian priors (Murray et al.,, 2010).

###  4.1 SMC for Online TP-MOE

In our proposed method, we use a sequential Monte Carlo sampler in order to
update the model as new data arrive (Del Moral et al.,, 2006). SMC follows
from importance sampling (IS) and sequential importance sampling (SIS)
algorithms in Monte Carlo methods, where IS and SIS sample the parameter of
interest from a proposal distribution in order to approximate an intractable
distribution:

$\displaystyle\int
P(X|\theta)P(\theta)\mbox{d}\theta\approx\frac{1}{J}\sum_{j=1}^{J}w^{(j)}\delta_{\theta^{(j)}}.$
(16)

However, IS and SIS suffer from the particle degeneracy problem where one
proposal weight, $w^{(j)}$, dominates the rest of the proposals. In SMC
methods, we resample the particles with probability equal to the proposal
weight. In this way, we replenish the sampler with particles that have high
weight and remove particles that have low weight.

For $j=1,\ldots,J$ particles, the particles
$(\mathbf{z}^{(j)},\boldsymbol{\theta}^{(j)},\mathbf{h}^{(j)},\alpha^{(j)})$
are updated as described before when a new observation arrives. Then, we
calculate the particle weights, which results in a posterior weighted sample
TP product-of-experts models. Initially when $i=1$, the particle $j$’s weight
is:

$\displaystyle w_{1}^{(j)}\propto
P(y_{1}|z_{1}^{(j)},\mathbf{x}_{1},\boldsymbol{\theta}^{(j)},h^{(j)},\nu^{(j)})P(\mathbf{x}_{1}|z_{1}^{(j)},\alpha^{(j)}).$
(17)

Then the updating procedure for $i>1$ is shown in Algorithm 1.

Input: New observation $(\mathbf{x}_{i},y_{i})$

for  _$j=1,\cdots,J$ in parallel_ do

Sample $z_{i}^{(j)}=k$ from $P(z_{i}^{(j)}|\alpha^{(j)},\mathbf{X}_{1:i-1})$

Sample $\alpha^{(j)}$ from the full conditional
$P(\alpha^{(j)}|\mathbf{z}_{1:i})$

Sample $\theta_{k}^{(j)}$ and $h_{k}^{(j)}$ jointly by using the elliptical
slice sampler

Sample $(k_{0}^{2})^{(j)}$ from
$P((k_{0}^{2})^{(j)}|h_{1}^{(j)},\ldots,h_{K}^{(j)})$

Sample $(\sigma_{k}^{2})^{(j)}$ from
$P(\sigma_{k}^{2}|\mathbf{X}_{k},\mathbf{y}_{k},\nu_{k}^{(j)})$

Sample $\nu_{k}^{(j)}$ by using the slice sampler

Update particle weight:  $\displaystyle w_{i}^{(j)}$
$\displaystyle=w_{i-1}^{(j)}P(\mathbf{x}_{i}|\alpha^{(j)},z_{i}^{(j)})$
$\displaystyle\times\frac{P\left(\mathbf{y}_{1:i}|\mathbf{X}_{1:i},\mathbf{\theta}_{k,i}^{(j)},h_{k}^{(j)},\nu_{k}^{(j)}\right)}{P\left(\mathbf{y}_{1:i-1}|\mathbf{X}_{1:i-1},\mathbf{\theta}_{k}^{\prime(j)},h_{k}^{\prime(j)},\nu_{k}^{\prime(j)}\right)}$
(18)

end for

Normalize weights:  $\displaystyle
w_{i}^{(j)}:=\frac{w_{i}^{(j)}}{\sum_{j=1}^{J}{w_{i}^{(j)}}}$ if  _$N_{eff}
<\frac{J}{2}$_ then

Resample particles
$(\mathbf{z}_{1:i}^{(\mathbf{j}^{*})},\mathbf{\theta}_{k}^{(\mathbf{j}^{*})},h_{k}^{(\mathbf{j}^{*})},\nu_{k}^{(\mathbf{j}^{*})},\alpha^{(\mathbf{j}^{*})})$,
where $\mathbf{j}^{*}\sim\mbox{Multinomial}(J,w_{i}^{(1)},\cdots,w_{i}^{(j)})$

Set $w_{i}^{(j)}:=\frac{1}{J}\mbox{ for }j=1,\cdots,J$

end if

Output: Particle weights $(w_{i}^{(1)},\cdots,w_{i}^{(j)})$ and particles
$(\mathbf{z}_{1:i}^{(1:J)},\mathbf{\theta}^{(1:J)},\mathbf{h}^{(1:J)},\mathbf{\nu}^{(1:J)},\alpha^{(1:J)})$

Algorithm 1 SMC Sampler for TP-MOE

for  _$j=1,\cdots,J$_ do

Predict new observations on particle j with  $\displaystyle p_{k}\propto
N_{k}\cdot P(\mathbf{x}_{*}|z_{*}=k,\mathbf{X}_{k},-)$ $\displaystyle
P(y_{*}^{(j)}|\mathbf{y},\mathbf{X},\mathbf{x}_{*},-)=$
$\displaystyle~{}~{}~{}~{}~{}\sum_{k\in\mathbf{K}^{+}}{p_{k}\cdot
P(y_{k*}^{(j)}|\mathbf{y}_{k},\mathbf{X}_{k},\mathbf{x}_{*},-)}$ (19)

end for

Average predictions:
$P(\bar{y}_{*}|\mathbf{y},\mathbf{X},\mathbf{x}_{*})=\sum_{j=1}^{J}{w_{i}^{(j)}P(y_{*}^{(j)}|\mathbf{y},\mathbf{X},-)}$

Algorithm 2 TP-MOE Prediction

The computational complexity is dominated by the inversion of a $N_{k}\times
N_{k}$ matrix. If we assume that the average size of $N_{k}$ is $N/K$, the
number of data divided by the number of clusters, the n the computational
complexity will be $\mathcal{O}(JN^{3}/K^{2})$. Under the basic setting of our
sampler, the complexity of the sampler still grows as new data arrive so the
method cannot truly be considered “online”. To this end, we adopt the
“minibatched” stochastic approximation that is widely used as a method for
substantially reducing the computational complexity of posterior inference
(Zhang et al.,, 2023; Zhang and Williamson,, 2019; Minsker et al.,, 2014;
Srivastava et al.,, 2015). A subsample of size $B$ from the mixture with
$N_{k}$ observations is drawn uniformly without replacement, then their
likelihood is calculated and upweighted by $N_{k}/B$ power to approximate the
full likelihood. The stochastic approximation method leads us to:

$\displaystyle\mathbf{u}_{k}^{(j)}\\!=\\!(u_{1},\ldots,u_{B})$
$\displaystyle~{}~{}~{}~{}~{}\sim\text{HyperGeometric}\big{(}B,\big{\\{}i:z_{i}^{(j)}=k\big{\\}}\big{)},$
$\displaystyle(\mathbf{y}_{\mathbf{u}_{k}},\mathbf{X}_{\mathbf{u}_{k}})=\big{(}y_{u},\mathbf{x}_{u}:u\in\mathbf{u}_{k}^{(j)}\big{)}.$
$\displaystyle
P\big{(}\mathbf{y}_{\mathbf{u}_{k}^{(j)}}|\mathbf{X}_{\mathbf{u}_{k}^{(j)}},\boldsymbol{\theta}_{k}^{(j)},h_{k},\sigma_{k}^{2}\big{)}$
$\displaystyle~{}~{}~{}~{}~{}\sim\mathcal{N}\left(0,\sigma_{k}^{2}(\mathbf{K}_{\boldsymbol{\theta}_{k}^{(j)}}+\frac{N_{k}|h_{k}|}{B}\mathbf{I})\right),$
$\displaystyle
P\big{(}\mathbf{y}_{\mathbf{u}_{k}^{(j)}}|\mathbf{X}_{\mathbf{u}_{k}^{(j)}},\boldsymbol{\theta}_{k}^{(j)},h_{k},\nu_{k}\big{)}$
$\displaystyle~{}~{}~{}~{}~{}\sim\mathcal{T}\left(\nu,0,\mathbf{K}_{\boldsymbol{\theta}_{k}^{(j)}}+\frac{N_{k}|h_{k}|}{B}\mathbf{I}\right).$
(20)

With minibatching, the complexity is reduced to
$\mathcal{O}(J\text{min}\\{N_{k},B\\}^{3}/K^{2})$. As each particle can be
updated independently, the parallel computation can be adopted to further
reduce the complexity to $\mathcal{O}(\text{min}\\{N_{k},B\\}^{3}/K^{2})$.
Then, we calculate the effective sample size,
$N_{eff}=1/\sum_{j=1}^{J}(w_{i}^{(j)})^{2}$, based on the particle weights. If
it is lower than a threshold set by the user, typically $J/2$, the particles
are resampled to only preserve the high-weighted ones. We make predictions
using a weighted average detailed in Algorithm 2.

##  5 Experiments

The choice of hyperparameters could significantly influence the TP-MOE’s
performance. According to Zhang and Williamson, (2019), larger $J$ and $B$
will lead to better performance but more computation time, while in contrast,
increasing $K$ will decrease both model performance and computation time. In
this section, we proceed to study the advantages of the heavy tails by
implementing the TP-MOE and other Gaussian-based models on different non-
stationary datasets and analysing their performances in terms of one-step-
ahead predictions. The GP models include a Gaussian mixture-of-experts model
(GP-MOE) (Zhang et al.,, 2023), a sparse online GP method using the Woodbury
identity and structured kernel interpolation (WISKI) (Stanton et al.,, 2021),
and an online sparse variational GP method (OSVGP) (Bui et al.,, 2017) 111 The
implementation for GP-MOE is available at
https://github.com/michaelzhang01/GPMOE. The code of OSVGP and WISKI are
available at: https://github.com/wjmaddox/online_gp. Our code will be
submitted in the supplementary material..

For the experiments, we sequentially predict the next future observation and
update the model with the real data point. The one-step predictive mean
squared error (MSE) is adopted to evaluate the results. The 5 datasets used
include: 1.) An accelerometer measurement of a motorcycle crash (N=94). 2.)
The price of Brent crude oil (N=100). 3.) The annual carbon dioxide output in
Canada (N=215). 4.) The annual water level of the Nile River data (N=100). 5.)
The exchange rate between the Euro and the US Dollar (N=200) 222The motorcycle
dataset can be found in the R package VarReg. The Brent, Canada CO2, and Nile
River datasets are available at: https://github.com/alanturing-institute/TCPD.
The EUR-USD dataset is available in the R package priceR.. The first three
exhibit non-stationarity in both length-scale and noise, while the Nile River
dataset shows only time-varying mean values and the exchange rate dataset is a
series of non-stationary noise. They have been pre-processed to have zero mean
and unit variance.

To make the results comparable, The TP-MOE and the GP-MOE share the same
particle number $J=100$ and the same $16$ cores used on a shared memory
process based on OpenMP, and the number of inducing points for all models is
set to be $50$. The OSVGP’s number of optimization iterations is set to the
default value of $1$. The radial basis function kernel for all models is:

$\displaystyle\Sigma(\mathbf{x},\mathbf{x}^{\prime})=\exp\left\\{-\frac{\theta}{2}\sum_{d=1}^{D}(x_{d}-x_{d}^{\prime})^{2}\right\\}$
(21)

The plots of four algorithms’ sample runs are shown in Figure 1-5, which
contain data points, one-step predictive mean (plotted with solid red lines)
and 95% predictive interval (plotted with dashed black lines). The data points
in TP-MOE’s and GP-MOE’s plots are coloured according to the cluster
assignment given by the particle with the highest weight. The results in terms
of the predictive MSE are listed in Table 1.

Table 1: One-step Predictive MSE. One Standard Error Reported in Parentheses | Motorcycle | Brent | Canada | Nile | EUR-USD
---|---|---|---|---|---
TP-MOE | 0.363 (0.028) | 0.146 (0.014) | 0.015 (0.003) | 0.738 (0.017) | 1.028 (0.013)
GP-MOE | 0.381 (0.038) | 0.160 (0.019) | 0.016 (0.004) | 0.752 (0.025) | 1.004 (0.009)
WISKI | 0.631 (0.000) | 0.220 (0.000) | 0.048 (0.000) | 0.767 (0.000) | 1.061 (0.000)
OSVGP | 0.998 (0.002) | 0.782 (0.021) | 0.711 (0.030) | 0.908 (0.008) | 1.019 (0.003)

From the comparisons in Table 1, we observe that our TP-MOE performs better
than the GP-based models and achieves lower predictive MSE. According to the
plots of sample runs (Figure 1-3), we can see that the MOE models can better
capture the heterogeneity of the underlying function better than the
stationary models. Moreover, we can see that the TP-MOE produces tighter
predictive credible intervals compared to the GP-MOE, which sometimes produces
overly conservative predictive intervals which suggests that the TP-MOE has
better uncertainty quantification capabilities. The OSVGP tends to underfit as
expected, while the WISKI cannot quantify the uncertainty as well as the
mixture-of-experts models due to its assumption of stationarity.

(a) TP-MOE

(b) GP-MOE

(c) WISKI

(d) OSVGP

Figure 1: Sample Runs on the Motorcycle Dataset. N = 94.

(a) TP-MOE

(b) GP-MOE

(c) WISKI

(d) OSVGP

Figure 2: Sample Runs on the Brent Dataset. N = 100.

(a) TP-MOE

(b) GP-MOE

(c) WISKI

(d) OSVGP

Figure 3: Sample Runs on the Canada Dataset. N = 215.

For the Nile River data which exhibits only non-stationary mean values, the
zero-mean assumption makes it hard for the TP-MOE and the GP-MOE to model the
trend. Despite the model misspecification, the TP-MOE still achieves the best
predictive MSE among four models according to Table 1. Also, based on the
comparison in Figure 4, the TP-MOE’s 95% predictive intervals are the most
consistent with the trend. The heavy-tailed property helps it be more robust
to the model misspecification.

(a) TP-MOE

(b) GP-MOE

(c) WISKI

(d) OSVGP

Figure 4: Sample Runs on the Nile River Dataset. N = 100.

However, when modelling the time-varying noise in the EUR-USD dataset, the GP-
MOE handles this task the best (Table 1). Figure 5 reveals the potential
reason for the failure of TP-MOE here. the predictive means are not very
stable, and tend to capture some trends that may not really exist. The OSVGP
which performs poorly in previous datasets even achieves better results this
time, because it usually underfits the streaming data and maintains stable
predictions.

(a) TP-MOE

(b) GP-MOE

(c) WISKI

(d) OSVGP

Figure 5: Sample Runs on the EUR-USD Dataset. N = 200.

Hence we conclude that, compared with the GP-based models, our TP-MOE can
better fit the data that exhibit non-stationarity in length scale and noise,
providing more accurate predictive means and most of the time better
uncertainty quantification. Also, the heavy tails help it be more robust to
the model misspecification. However, if the noise dominates the series, the
kernelized covariance is less useful to capture such a trend.

##  6 Conclision

Heavy-tailed data sets appear in a wide variety of applied settings. However,
devising models that can adequately handle their noise structure is not
trivial. In this paper, we build a Bayesian mixture of student-$t$ processes
model with an overall-local scale structure for noisy data, which can be
inferred by an SMC online algorithm. We have shown that TP-MOE has advantages
over the Gassian-based models when facing commonly encountered non-stationary
data.

In future work, we are interested in applying the TP-MOE in optimization and
reinforcement learning tasks. For such tasks, the learning, prediction, and
decision making aspects of the model occur in sparse, noisy environments that
require heavy-tailed models in order for a learning agent to properly handle
the problem at hand. Modelling data with a mixture of Student-$t$ processes is
a natural method for dealing with non-stationarity and heavy-tailed errors yet
their popularity has still eluded the machine learning community. We seek to
fill that gap with the method proposed in this paper.

## References

  * Aldous, (1985) Aldous, D. (1985).  Exchangeability and related topics.  École d’Été de Probabilités de Saint-Flour XIII—1983, pages 1–198.
  * Antoniak, (1974) Antoniak, C. E. (1974).  Mixtures of Dirichlet processes with applications to Bayesian nonparametric problems.  Annals of Statistics, pages 1152–1174.
  * Bauer et al., (2016) Bauer, M., van der Wilk, M., and Rasmussen, C. E. (2016).  Understanding probabilistic sparse Gaussian process approximations.  Advances in Neural Information Processing Systems, 29.
  * Bui et al., (2017) Bui, T. D., Nguyen, C., and Turner, R. E. (2017).  Streaming sparse gaussian process approximations.  In Guyon, I., Luxburg, U. V., Bengio, S., Wallach, H., Fergus, R., Vishwanathan, S., and Garnett, R., editors, Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc.
  * Cohen et al., (2020) Cohen, S., Mbuvha, R., Marwala, T., and Deisenroth, M. (2020).  Healing products of Gaussian process experts.  In International Conference on Machine Learning, pages 2068–2077. PMLR.
  * Csató and Opper, (2002) Csató, L. and Opper, M. (2002).  Sparse On-Line Gaussian Processes.  Neural Computation, 14(3):641–668.
  * Damien et al., (1999) Damien, P., Wakefield, J., and Walker, S. (1999).  Gibbs sampling for Bayesian non-conjugate and hierarchical models by using auxiliary variables.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 61(2):331–344.
  * Deisenroth and Ng, (2015) Deisenroth, M. and Ng, J. W. (2015).  Distributed Gaussian processes.  In International Conference on Machine Learning, pages 1481–1490. PMLR.
  * Del Moral et al., (2006) Del Moral, P., Doucet, A., and Jasra, A. (2006).  Sequential Monte Carlo samplers.  Journal of the Royal Statistical Society Series B: Statistical Methodology, 68(3):411–436.
  * Escobar and West, (1995) Escobar, M. D. and West, M. (1995).  Bayesian density estimation and inference using mixtures.  Journal of the American Statistical Association, 90(430):577–588.
  * Fernández and Steel, (1999) Fernández, C. and Steel, M. F. (1999).  Multivariate student-$t$ regression models: Pitfalls and inference.  Biometrika, 86(1):153–167.
  * Geweke, (1993) Geweke, J. (1993).  Bayesian treatment of the independent student-$t$ linear model.  Journal of Applied Econometrics, 8(S1):S19–S40.
  * Gramacy and Polson, (2011) Gramacy, R. B. and Polson, N. G. (2011).  Particle learning of Gaussian process models for sequential design and optimization.  Journal of Computational and Graphical Statistics, 20(1):102–118.
  * Härkönen et al., (2022) Härkönen, T., Wade, S., Law, K., and Roininen, L. (2022).  Mixtures of Gaussian process experts with SMC2.  arXiv preprint arXiv:2208.12830.
  * Juárez and Steel, (2010) Juárez, M. A. and Steel, M. F. J. (2010).  Model-based clustering of non-gaussian panel data based on skew-t distributions.  Journal of Business & Economic Statistics, 28(1):52–66.
  * Jylänki et al., (2011) Jylänki, P., Vanhatalo, J., and Vehtari, A. (2011).  Robust $g$aussian process regression with a student-$t$ likelihood.  Journal of Machine Learning Research, 12(11).
  * Meeds and Osindero, (2005) Meeds, E. and Osindero, S. (2005).  An alternative infinite mixture of Gaussian process experts.  Advances in Neural Information Processing Systems, 18.
  * Minsker et al., (2014) Minsker, S., Srivastava, S., Lin, L., and Dunson, D. (2014).  Scalable and robust bayesian inference via the median posterior.  In Xing, E. P. and Jebara, T., editors, Proceedings of the 31st International Conference on Machine Learning, volume 32 of Proceedings of Machine Learning Research, pages 1656–1664, Bejing, China. PMLR.
  * Murray et al., (2010) Murray, I., Adams, R., and MacKay, D. (2010).  Elliptical slice sampling.  In Teh, Y. W. and Titterington, M., editors, Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics, volume 9 of Proceedings of Machine Learning Research, pages 541–548, Chia Laguna Resort, Sardinia, Italy. PMLR.
  * Neal, (2003) Neal, R. M. (2003).  Slice sampling.  Annals of Statistics, 31(3):705–767.
  * Nguyen-tuong et al., (2008) Nguyen-tuong, D., Peters, J., and Seeger, M. (2008).  Local gaussian process regression for real time online model learning.  In Koller, D., Schuurmans, D., Bengio, Y., and Bottou, L., editors, Advances in Neural Information Processing Systems, volume 21. Curran Associates, Inc.
  * Rasmussen and Ghahramani, (2001) Rasmussen, C. and Ghahramani, Z. (2001).  Infinite mixtures of Gaussian process experts.  Advances in Neural Information Processing Systems, 14.
  * Rasmussen and Williams, (2005) Rasmussen, C. E. and Williams, C. K. I. (2005).  Gaussian Processes for Machine Learning.  The MIT Press.
  * Shah et al., (2014) Shah, A., Wilson, A., and Ghahramani, Z. (2014).  Student-t Processes as Alternatives to Gaussian Processes.  In Kaski, S. and Corander, J., editors, Proceedings of the Seventeenth International Conference on Artificial Intelligence and Statistics, volume 33 of Proceedings of Machine Learning Research, pages 877–885, Reykjavik, Iceland. PMLR.
  * Snelson and Ghahramani, (2006) Snelson, E. and Ghahramani, Z. (2006).  Sparse Gaussian processes using pseudo-inputs.  Advances in Neural Information Processing Systems, 18:1259–1266.
  * Srivastava et al., (2015) Srivastava, S., Cevher, V., Dinh, Q., and Dunson, D. (2015).  WASP: Scalable Bayes via barycenters of subset posteriors.  In Lebanon, G. and Vishwanathan, S. V. N., editors, Proceedings of the Eighteenth International Conference on Artificial Intelligence and Statistics, volume 38 of Proceedings of Machine Learning Research, pages 912–920, San Diego, California, USA. PMLR.
  * Stanton et al., (2021) Stanton, S., Maddox, W., Delbridge, I., and Gordon Wilson, A. (2021).  Kernel interpolation for scalable online gaussian processes.  In Banerjee, A. and Fukumizu, K., editors, Proceedings of The 24th International Conference on Artificial Intelligence and Statistics, volume 130 of Proceedings of Machine Learning Research, pages 3133–3141. PMLR.
  * Svensson et al., (2015) Svensson, A., Dahlin, J., and Schön, T. B. (2015).  Marginalizing Gaussian process hyperparameters using sequential Monte Carlo.  In 2015 IEEE 6th International Workshop on Computational Advances in Multi-Sensor Adaptive Processing (CAMSAP), pages 477–480. IEEE.
  * Tang et al., (2017) Tang, Q., Niu, L., Wang, Y., Dai, T., An, W., Cai, J., and Xia, S.-T. (2017).  Student-$t$ process regression with student-$t$ likelihood.  In International Joint Conferences on Artificial Intelligence, pages 2822–2828.
  * Tipping, (2001) Tipping, M. E. (2001).  Sparse Bayesian learning and the relevance vector machine.  Journal of Machine Learning Research, 1(Jun):211–244.
  * Titsias, (2009) Titsias, M. (2009).  Variational learning of inducing variables in sparse Gaussian processes.  In Artificial Intelligence and Statistics, pages 567–574. PMLR.
  * Vanhatalo et al., (2009) Vanhatalo, J., Jylänki, P., and Vehtari, A. (2009).  Gaussian process regression with student-$t$ likelihood.  Advances in Neural Information Processing Systems, 22.
  * West, (1984) West, M. (1984).  Outlier models and prior distributions in Bayesian linear regression.  Journal of the Royal Statistical Society Series B: Statistical Methodology, 46(3):431–439.
  * Wilson and Nickisch, (2015) Wilson, A. and Nickisch, H. (2015).  Kernel interpolation for scalable structured gaussian processes (kiss-gp).  In Bach, F. and Blei, D., editors, Proceedings of the 32nd International Conference on Machine Learning, volume 37 of Proceedings of Machine Learning Research, pages 1775–1784, Lille, France. PMLR.
  * Zhang et al., (2023) Zhang, M. M., Dumitrascu, B., Williamson, S. A., and Engelhardt, B. E. (2023).  Sequential gaussian processes for online learning of nonstationary functions.  IEEE Transactions on Signal Processing, 71:1539–1550.
  * Zhang and Williamson, (2019) Zhang, M. M. and Williamson, S. A. (2019).  Embarrassingly parallel inference for gaussian processes.  Journal of Machine Learning Research, 20(169):1–26.
  * Zhang and Yeung, (2010) Zhang, Y. and Yeung, D. (2010).  Multi-task learning using generalized t process.  In Teh, Y. W. and Titterington, M., editors, Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics, volume 9 of Proceedings of Machine Learning Research, pages 964–971, Chia Laguna Resort, Sardinia, Italy. PMLR.