[1,2]Katie Buchhorn

[1]School of Mathematical Sciences, Queensland University of Technology,
George Street, Brisbane, 4000, Queensland, Australia

2]Centre for Data Science (CDS), Queensland University of Technology, George
Street, Brisbane, 4000, Queensland, Australia

# Bayesian Design for Sampling Anomalous Spatio-Temporal Data

<EMAIL_ADDRESS>Kerrie Mengersen<EMAIL_ADDRESS>Edgar Santos-
Fernandez<EMAIL_ADDRESS>James McGree
<EMAIL_ADDRESS>*  [

###### Abstract

Data collected from arrays of sensors are essential for informed decision-
making in various systems. However, the presence of anomalies can compromise
the accuracy and reliability of insights drawn from the collected data or
information obtained via statistical analysis. This study aims to develop a
robust Bayesian optimal experimental design (BOED) framework with anomaly
detection methods for high-quality data collection. We introduce a general
framework that involves anomaly generation, detection and error scoring when
searching for an optimal design. This method is demonstrated using two
comprehensive simulated case studies: the first study uses a spatial dataset,
and the second uses a spatio-temporal river network dataset. As a baseline
approach, we employed a commonly used prediction-based utility function based
on minimising errors. Results illustrate the trade-off between predictive
accuracy and anomaly detection performance for our method under various design
scenarios. An optimal design robust to anomalies ensures the collection and
analysis of more trustworthy data, playing a crucial role in understanding the
dynamics of complex systems such as the environment, therefore enabling
informed decisions in monitoring, management, and response.

###### keywords:

Anomaly Detection, Bayesian Design, Optimal Experimental Design, Robust
Design, Sensor Data, Spatio-temporal Model, Spatial Model

##  1 Introduction

Optimal experimental design can be considered intelligent data collection [1],
and while high-quality data underpins many functions of modern society, there
remains a lack of methodologies that account for data quality within a design
framework. In-situ sensor technology has revolutionised data collection across
various domains, including environmental monitoring of air, water and soil
quality. Bayesian optimal experimental design (BOED) provides a natural
mechanism for optimising and automating the data collection process, as a
model-based framework aiming to maximise the amount of information gathered
from such experiments. In this paper, we develop a robust BOED framework that
combines complex spatial and spatio-temporal modelling with anomaly detection
methods, ensuring the data collected and analysed at design sites provide
useful information (e.g. accurate predictions) as well as automated and
reliable anomaly detection.

The principles of optimal design are evident across various domains, including
economics [2], ecology [3], social sciences [4, 5], physics [6, 7], and
healthcare [8], among other quantitative fields. Optimal experimental design
serves as a means to attain maximum information with respect to an
experimental objective, and is particularly important when such experiments
are costly in terms of equipment and time. Improving data quality through
optimal design methodologies remains a largely unexplored area of research.
When designing microarray experiments, [9] considered the quality of gene
expression data alongside other image processing analyses. The research by
[10] focused on optimal design in production systems. It balanced sunk costs
and revenue increments due to process quality improvements, based on Taguchi
cost functions. Even in clinical trial settings, where optimal design
methodologies are widely applied, there has been limited research on
experimental objectives aimed at enhancing the quality and quantity of data
collected from (for example) questionnaires [11]. While there is evidence for
design features that improve data completeness [12], there is a recognised
need for further research to evaluate these strategies. In this work, we
account for an often overlooked consideration: potential anomalies in the data
collected from the optimal design.

The expanding body of research on optimal design for environmental monitoring
underscores the significance of innovative statistical methodologies and
optimal designs for complex ecosystems with multifaceted objectives. In
particular, [13] developed a spatial statistical methodology for designing
cost-effective national air pollution monitoring networks, accounting for the
nonstationary nature of atmospheric processes. Adaptive design approaches have
also been proposed for monitoring coral reefs [14, 15]. In survey design, [16]
employed a Bayesian approach for a multi-attribute valuation of environmental
actions for landscape conservation and improvement. Authors [17] presented a
methodology for designing optimal solvent blends for minimising environmental
impacts, balancing operational and environmental constraints. Later, [18]
reviewed optimal design in supply chain management for chemical processes,
focusing on energy efficiency, waste management, and sustainable water
management. These studies collectively emphasise the importance of Bayesian
approaches and the necessity for balancing accuracy, cost, and environmental
impact in environmental monitoring and management. Despite the evident demand,
wide-scale adoption of BOED for complex, real-world ecosystems has remained
relatively limited due to significant computational challenges, particularly
in adaptive scenarios or when dealing with complex models [19, 20].

River and stream ecosystems, which are crucial for both ecological habitats
and economic activities, are increasingly threatened by climate change,
pollution, and human activities. High-frequency data from sensors offer
opportunities to understand and manage spatio-temporal dynamics of stream
attributes. Unlike many spatial applications, data collected on stream
networks exhibit highly complex and multi-layered spatial dependencies,
influenced by factors like climate gradients, organism mobility, and within-
network physiochemical and biological processes [21]. River network system
modelling therefore must account for complex covariance relationships, also
due to the branching network topology and unidirectional flow of water,
resulting in the passive movement of materials, nutrients, and organisms
downstream. Only recently has BOED been applied to river network systems, but
for a spatial model only [22]. In this paper we explore an extension of BOED
for spatio-temporal river network models.

Anomalies in sensor data can indicate critical events such as extreme weather
conditions, pollution, or equipment failure. Rapid automated methods for
identifying and determining the probable source of such anomalies are
therefore vital in safeguarding lives and valuable assets. In particular,
distinguishing between anomalies that arise from sensor malfunctions, herein
defined as technical anomalies, from extreme river events is essential.
Existing approaches for anomaly detection in river network environments
include autoregressive models, e.g. autoregressive integrated moving average
(ARIMA) [23], machine learning methods involving artificial neural networks
(ANN), random forests [24] and long short-term memory (LSTM) [25]. ARIMA
captures temporal dependencies only, shown to be suitable in understanding
trends and cyclical patterns in water quality parameters. The machine learning
techniques considered both univariate and multivariate dependencies. However,
such methods lack the ability to provide a holistic understanding of water
quality dynamics. Simultaneously accounting for spatial and temporal variation
was considered in subsequent work by [26], showcasing the effectiveness of
spatially aware models (posterior predictive distributions, finite mixtures,
and Hidden Markov Models) in capturing extreme river events occurring across
multiple locations. Also, [27] use a graph neural network model to capture
inter-sensor relationships as a learned graph, and employs graph attention-
based forecasting to predict future sensor behaviour. However, attributing the
origin of an anomaly to the correct sensor is challenging using this method. A
recent feature-based method called _oddstream_ , suitable for nonstationary
streaming time-series data, is proposed by [28]. Oddstream uses kernel density
estimates on a 2-dimensional projection with extreme value theory to calculate
a threshold boundary to detect outliers in time series data.

Bayesian design principles with optimised anomaly detection represent a novel
and promising avenue for advancing the real-world application in experimental
design. In this paper, we consider the context of in-situ, online quality
control systems commonly used in system monitoring. In complex systems like
river networks, strategically placed sensors can better identify anomalous
data, discrepancies can be promptly detected and addressed, ensuring the
reliability of the collected data.

###  1.1 Bayesian design

BOED is a statistical approach that integrates Bayesian inference into the
process of optimal experimental design to quantify (and make decisions under)
uncertainty in a coherent and systematic manner. In a broad sense, a design
$\boldsymbol{d}$ refers to a systematic arrangement and selection of variables
or factors to explore (and optimise) across all possible configurations or
values, denoted by the design space, $\mathcal{D}$. The Bayesian design
process begins with the formulation of a statistical model that represents the
underlying system or phenomenon being studied. This model typically includes a
prior distribution, $p(\boldsymbol{\theta})$, which can be informed by
previous studies, expert opinion, or theoretical considerations, and a
likelihood function, $p(\boldsymbol{y}|\boldsymbol{\theta},\boldsymbol{d})$,
which describes how the observed data are generated given the parameters. A
utility function, $u(\cdot)$, is used to quantify the value of different
experimental outcomes from a given design. The utility function can take
various forms, depending on the specific goals of the experiment, such as
maximizing information gain or minimizing cost or error. In the first case,
the utility function may be defined in terms of the reduction in uncertainty
about the parameters [e.g., decrease in the posterior variance 29, 30]. In the
second case, the utility may be a function of the financial cost, time, or
other resources associated with the experiment, or it may be related to the
accuracy of experimental results. The overall objective in BOED is to choose
the design, $\boldsymbol{d}$, from the design space, $\mathcal{D}$, that
maximises the expected utility, which is computed as the weighted average of
the utility over all possible outcomes $(\Theta,\mathcal{Y})$. The Bayesian
optimal design $\boldsymbol{d}_{*}$ can therefore be expressed as,

$\displaystyle\boldsymbol{d}_{*}$ $\displaystyle=\mbox{arg
max}_{\boldsymbol{d}\in\mathcal{D}}\mathbb{E}[u(\boldsymbol{d},\boldsymbol{\theta},\boldsymbol{y})]$
$\displaystyle=\mbox{arg
max}_{\boldsymbol{d}\in\mathcal{D}}\int_{\boldsymbol{y}\in\mathcal{Y}}\int_{\boldsymbol{\theta}\in\Theta}u(\boldsymbol{d},\boldsymbol{\theta},\boldsymbol{y})\
p(\boldsymbol{\theta},\boldsymbol{y}|\boldsymbol{d})\ d\boldsymbol{\theta}\
d\boldsymbol{y}.$ (1)

It is often challenging to find Bayesian optimal designs for realistic
problems, as the expected utility is typically intractable and the design
space may be high-dimensional. In the next section, we introduce techniques to
handle such challenges as well as the statistical models considered throughout
this paper.

##  2 Background

This section begins with a description of the motivating case study, namely
optimising the design specifying the location of in-situ sensors across a
river network. Details of spatio-temporal models used for river networks are
given, followed by an overview of the coordinate exchange algorithm used for
optimisation throughout this paper.

###  2.1 Bayesian river network model

Generally, spatio-temporal linear regression models are formulated as:

$\boldsymbol{y}=\mathbf{X}\boldsymbol{\beta}+\boldsymbol{v}+\boldsymbol{\epsilon},$
(2)

where $\boldsymbol{y}$ is a stacked response vector of size $n\times 1$, with
$n=S\times T$ (for $S$ spatial locations and $T$ time points), $\mathbf{X}$ is
a $n\times p$ design matrix of $p$ covariates, $\boldsymbol{\beta}$ is a
$p\times 1$ vector or regression coefficients, and $\boldsymbol{v}$ is a
stacked vector of spatio-temporal autocorrelated random effects, and
$\boldsymbol{\epsilon}$ is the unstructured error term with,
$\text{var}(\boldsymbol{\epsilon})=\sigma_{0}^{2}\bf{I}$, where $\bf{I}$ is
the identity matrix and $\sigma_{0}^{2}$ is called the nugget effect.

####  2.1.1 Spatial covariance for river networks

In geostatistics, it is standard to adopt covariance models based on the
Euclidean distance between two spatial locations $x_{i}$ and $x_{j}$, such as
the exponential (Equation 3), Gaussian (Equation 4), and spherical functions
(Equation 5):

$\displaystyle\text{Exponential model:}\quad
C_{\text{Euc}}(h_{e};\sigma^{2}_{e},\rho_{e})=\sigma^{2}_{e}\exp\left(-\frac{3h_{e}}{\rho_{e}}\right),$
(3) $\displaystyle\text{Gaussian model:}\quad
C_{\text{Euc}}(h_{e};\sigma^{2}_{e},\rho_{e})=\sigma^{2}_{e}\exp\left(-3\left(\frac{h_{e}}{\rho_{e}}\right)^{2}\right),$
(4) $\displaystyle\text{Spherical model:}\quad
C_{\text{Euc}}(h_{e};\sigma^{2}_{e},\rho_{e})=\sigma^{2}_{e}\left(1-\left(\frac{3d}{2\rho_{e}}\right)+\left(\frac{d^{3}}{2\rho_{e}^{3}}\right)\right)\mathbbm{1}\left(\frac{h_{e}}{\rho_{e}}\leq
1\right).$ (5)

Here, $h_{e}$ represents the Euclidean distance between locations,
$\sigma^{2}_{e}$ is the partial sill and $\rho_{e}$ is the spatial range
parameter. Consider two locations, $s_{i}$ and $s_{j}$, on a river network,
$\mathcal{S}$. These models may not effectively describe spatial correlation
between $s_{i}$ and $s_{j}$, considering the unique branching topology,
connectivity, and water flow volume of the river network within which they
exist. To address this, tail-up (Equation 6) and tail-down (Equation 7) models
have been proposed [31]. Exponential tail-up models restrict autocorrelation
to flow-connected sites only, using spatial weights based on the flow volume
and branching structure. Exponential tail-down models consider both flow-
connected and flow-unconnected sites, allowing spatial dependence even between
unconnected locations:

$\displaystyle C_{\text{TU}}(h;\sigma^{2}_{u},\alpha_{u})$
$\displaystyle=\begin{cases}\sigma^{2}_{u}\exp\left(-\frac{3h}{\alpha_{u}}\right)&\text{if
sites are flow-connected,}\\\ 0&\text{if sites are flow-
unconnected,}\end{cases}$ (6) $\displaystyle
C_{\text{TD}}(h;\sigma^{2}_{d},\alpha_{d})$
$\displaystyle=\sigma^{2}_{d}\exp\left(-\frac{3h}{\alpha_{d}}\right),$ (7)

where $h$ is the distance between $s_{i}$ and $s_{j}$ measured _along_ the
stream network. For the purely spatial case, when $T=1$ in Equation (2), the
vector of random effects $\boldsymbol{v}$ is only spatially related with
$\boldsymbol{\Sigma}=\text{cov}(\boldsymbol{v})$. A covariance mixture
approach for capturing the unique spatial patterns across stream networks
combines these models, accounting for Euclidean (Euc), tail-up (TU) and tail-
down (TD) components as follows,

$\boldsymbol{\Sigma}=\mathbf{C}_{\text{Euc}}+\mathbf{W}\bf{C}_{\text{TU}}+\mathbf{C}_{\text{TD}},$
(8)

where $\mathbf{C}_{Euc}$, $\mathbf{C}_{\text{TU}}$ and
$\mathbf{C}_{\text{TD}}$ are matrices derived from their corresponding
covariance functions, and $\mathbf{W}$ is a symmetric squared matrix
indicating the weights between sites used to proportionally allocate (i.e.
split) the tail-up moving average function at river junctions. Denote the
vector of spatial parameters by
$\boldsymbol{\gamma}=(\sigma^{2}_{e},\alpha_{e},\sigma^{2}_{u},\alpha_{u},\sigma^{2}_{d},\alpha_{d},\sigma^{2}_{0})$.

####  2.1.2 Spatio-temporal river network model

In this spatio-temporal analysis of stream networks, we consider repeated
measures at times $t=1,...,T$ at fixed spatial locations $s=1,...,S$.
Following [32], for continuous response variables we assume the following
Bayesian hierarchical model,

$p(\boldsymbol{y}_{1},\boldsymbol{y}_{2},...,\boldsymbol{y}_{T}|\boldsymbol{\theta},\boldsymbol{d})=\prod_{t=2}^{T}p(\boldsymbol{y}_{t}|\boldsymbol{y}_{t-1},\boldsymbol{\gamma},\mathbf{X}_{t},\mathbf{X}_{t-1},\boldsymbol{\beta},\boldsymbol{\Phi},\boldsymbol{\Sigma},\boldsymbol{d})p(\boldsymbol{y}_{1}|\boldsymbol{\theta},\boldsymbol{d}),$
(9)

where $\boldsymbol{y}_{1}$ is the process at $t=1$, $\boldsymbol{d}$ is the
design, and

$\displaystyle
p(\boldsymbol{y}_{t}|\boldsymbol{y}_{t-1},\boldsymbol{\gamma},\mathbf{X}_{t},\mathbf{X}_{t-1},\boldsymbol{\beta},\boldsymbol{\Phi},\boldsymbol{\Sigma},\boldsymbol{d})$
$\displaystyle=\mathcal{N}(\boldsymbol{\mu}_{t},\boldsymbol{\Sigma}+\sigma_{0}^{2}\mathbf{I})$
(10) $\displaystyle\boldsymbol{\mu}_{t}$
$\displaystyle=\mathbf{X}_{t}\boldsymbol{\beta}+\boldsymbol{\Phi}(\boldsymbol{y}_{t-1}-\mathbf{X}_{t-1}\boldsymbol{\beta}),$
(11)

where $\mathbf{X}_{t}$ is the design matrix of covariates and an intercept.
Here, we let $\boldsymbol{\theta}$ collectively represent the model parameters
including spatial parameters, $\boldsymbol{\gamma}$, regression parameters,
$\boldsymbol{\beta}$, and autoregressive parameters, $\boldsymbol{\phi}$. This
dynamical model assumes constant spatial correlation over time intervals and
relies on first-order Markovian dependence, where
$\boldsymbol{\Sigma}=\text{cov}(\boldsymbol{v}_{t})\ \forall t=2,...,T$ is the
$S\times S$ spatial covariance matrix, e.g., formulated by the exponential
function in Equation (8), and $\boldsymbol{\Phi}$ is the $S\times S$
transition matrix determining the amount of temporal autocorrelation, e.g.,
assuming the same temporal autocorrelation for all spatial locations. In this
case, the diagonal elements of $\boldsymbol{\Phi}$ are all equal to $\phi$ and
all the off-diagonal ones are set to zero, i.e.,

$\boldsymbol{\Phi}=\begin{bmatrix}\phi&0&0&\cdots&0\\\ 0&\phi&0&\cdots&0\\\
0&0&\phi&\cdots&0\\\ \vdots&\vdots&\vdots&\ddots&\vdots\\\
0&0&0&\cdots&\phi\\\ \end{bmatrix}.$

An alternative approach to model spatio-temporal autocorrelation involves
constructing the full space-time covariance function. However the dynamical
model described above that models the evolution of a spatial process is more
computationally efficient, by operating with spatial covariance matrices
instead of joint space-time matrices, since the computational bottleneck in
these methods is the inversion of large covariance matrices [32]. In fact, the
authors of [32] show these two models are equivalent mathematically.

###  2.2 Finding the optimal design

Finding the Bayesian optimal design in Equation (1) requires evaluation of the
integrals and an approach to search the design space for the maximisation
problem. The evaluation of the expected utility,
$\text{U}(\boldsymbol{d}):=\mathbb{E}[u(\boldsymbol{d},\boldsymbol{\theta},\boldsymbol{y})]$,
is analytically intractable for the linear mixed-model in Equation (2), as
with most cases. In practice, an approximation to the expected utility,
$\hat{\text{U}}(\boldsymbol{d})$, is typically evaluated using Monte Carlo
integration, as follows:

$\hat{\text{U}}(\boldsymbol{d})=\frac{1}{B}\sum^{B}_{b=1}u(\boldsymbol{d},\boldsymbol{\theta}^{(b)},\boldsymbol{y}^{(b)}),$
(12)

with draws from the prior $\boldsymbol{\theta}^{(b)}\sim
p(\boldsymbol{\theta})$ and then the likelihood $\boldsymbol{y}^{(b)}\sim
p(\boldsymbol{y}|\boldsymbol{\theta}^{(b)},\boldsymbol{d})$.

For complex experimental settings, the Approximate Coordinate Exchange (ACE)
[33] is an effective approach in searching the design space. The ACE method
utilises the approximate coordinate exchange algorithm, and simplifies the
optimisation of (an approximation to) the expected utility in high-dimensional
design spaces through Gaussian process (GP) regression models. Assuming a GP
prior, a one-dimensional emulator for the $ij$-th coordinate of
$\tilde{\text{U}}(\boldsymbol{d})$ is constructed as the mean of the posterior
predictive distribution conditioned on a few evaluations of
$\hat{\text{U}}(\boldsymbol{d})$. This emulation is calculated for $i=1,...,n$
and $j=1,...,k$, and then the entire process is iterated for a total of
$N_{1}$ times. See Appendix A.2 for further details. To mitigate the impact of
a poor emulator, the proposal design $\boldsymbol{d}^{*}_{ij}$ is only
accepted as the next design with a probability $p^{*}$, determined by a
Bayesian hypothesis test. Here, $p^{*}$ is calculated as the posterior
probability that the expected utility for the proposed design,
$\boldsymbol{d}^{*}_{ij}$ is greater than that of the current design.
Typically, larger Monte Carlo sample sizes are used for these tests than for
the construction of the emulator, enhancing the precision of the
approximation, $\hat{\text{U}}_{ij}(\boldsymbol{d})$. Unless otherwise stated,
we use Monte Carlo sample sizes $B_{1}=1500$ and $B_{2}=1000$, respectively,
in the following experiments. This approach is adept at handling generalised
linear and nonlinear models, which significantly broadens the capability and
practicality of BOED for a variety of objectives, providing a strong basis for
our methods application.

##  3 Method

In this section, we introduce our novel approach that incorporates anomaly
detection into the Bayesian optimal design framework to improve the ability to
detect anomalies across a network. The method is introduced generally,
followed by a detailed discussion of each component within the context of
environmental monitoring. Finally, an example implementation is presented
which is later utilised in the case studies.

###  3.1 Bayesian design with anomaly detection

The below algorithm serves as a general schema for Bayesian design,
incorporating anomaly generation and detection techniques, while introducing a
dual-purpose utility aiming to balance between the objectives of the design
and the efficacy of anomaly detection.

Input: Design, anomaly generator, anomaly detector, utility function.

Output: Vector of utilities.

1 Generate data $\boldsymbol{y}$ (and possibly
$\boldsymbol{y}_{\text{train}}$) based on design $\boldsymbol{d}$.

2 Generate anomalies and contaminate data $\boldsymbol{y}$ to produce
$\boldsymbol{y}_{\text{anom}}$.

3 Detect anomalies in $\boldsymbol{y}_{\text{anom}}$ (possibly using
$\boldsymbol{y}_{\text{train}}$), and remove (or impute) predicted anomalies
in $\boldsymbol{y}_{\text{anom}}$ to produce
$\boldsymbol{y}_{\text{cleaned}}$.

4 Evaluate performance of the anomaly detector.

5 Evaluate objective.

6 Evaluate dual-purpose utility.

Algorithm 1 General schema for Bayesian design with anomaly detection

####  3.1.1 General schema

Algorithm 1 begins by taking inputs: the design parameters, an anomaly
generator, an anomaly detector, and a utility function. Data $\boldsymbol{y}$,
and optionally, a training set $\boldsymbol{y}_{\text{train}}$, are generated
from the likelihood function based on the given design $\boldsymbol{d}$. Then,
anomalies are generated to contaminate the data $\boldsymbol{y}$ to create the
data variant $\boldsymbol{y}_{\text{anom}}$. Anomalies are detected given
$\boldsymbol{y}_{\text{anom}}$, potentially utilising
$\boldsymbol{y}_{\text{train}}$ for training. Detected anomalies are then
removed or imputed from $\boldsymbol{y}_{\text{anom}}$ to produce
$\boldsymbol{y}_{\text{cleaned}}$. Finally, the utility is evaluated based on
the accuracy of the anomaly detector and experimental design objective. The
following sections discusses each of these components in more detail.

####  3.1.2 Anomaly generation

While data collected in environmental systems exhibit complex patterns of
spatio-temporal dependence, here we assume that sensor-related anomalies occur
independently of space, but may persist over time (indicating the potential
need for sensor maintenance or battery failure). Consider an indicator matrix,
$\mathbf{A}$, with elements, $a_{ti}\in\\{0,1\\}$, where anomalies are denoted
as 1. Persistent anomalies (i.e., consecutive anomalous observations across
time) are simulated by setting $a_{\text{A};ti},\ldots,a_{\text{A};(t+l)i}=1$,
with,

$\displaystyle a_{\text{A};ti}\sim$
$\displaystyle\text{Bernoulli}(p_{\text{A}})$ anomaly frequency $\displaystyle
l\sim$ $\displaystyle\text{Poisson}(\lambda_{\text{A}})$ length of anomaly

In the spatial case, when $T=1$, we set $l=1$.

####  3.1.3 Anomaly detection

In this section, we introduce two different anomaly detection methods, noting
that other anomaly detection algorithms may also be used [34, 35, 36, 37].

_Spatial anomaly detection:_ The approach described in Algorithm 2 effectively
identifies anomalies in spatial data by evaluating the mean of each point’s
nearest neighbours and setting thresholds based on the standard deviation of
the sensor from the training data, flagging points as anomalous if they fall
outside these defined bounds.

Input: training data, $\boldsymbol{y}_{\text{train}}$, and test data,
$\boldsymbol{y}_{\text{anom}}$.

Output: Indicator vector $\boldsymbol{a}\in\mathbb{R}^{n}$ for $n$ spatial
locations, with $a_{i}=1$ indicating an anomaly for sensor $i$.

1 for  _$i\in 1:n$_ do

2 Find the k-nearest neighbours of ${y}_{\text{anom},i}$ based on Euclidean
distance, $\text{NN}_{k}({y}_{\text{anom},i})$.

3 Calculate $\bar{{y}}_{i}$ as the mean of the k-nearest neighbours,
$\bar{{y}}_{i}=\frac{1}{k}\sum_{y\in\text{NN}_{k}({y}_{\text{anom},i})}y$.

4 Define thresholds $\tau^{+}=\bar{y}_{i}+3\sigma_{i}$ and
$\tau^{-}=\bar{y}_{i}-3\sigma_{i}$, where $\sigma_{i}$ is the standard
deviation of the training set $\boldsymbol{y}_{\text{train},i}$.

5 if  _$\tau^{-} <y_{i}<\tau^{+}$ _ then

6 Set $a_{i}=0$

7

8 else

9 Set $a_{i}=1$

10

11

Algorithm 2 Anomaly Detection in Spatial Data

_Spatio-temporal anomaly detection_ : Oddstream [38] is an anomaly detection
method that uses extreme value theory and feature-based time series analysis
for early detection of anomalies in non-stationary streaming time series data.
The Oddstream method in the Appendix as Algorithm 4 starts with the training
phase where it extracts features, normalises them, applies PCA, estimates the
probability density of the subspace, and determines the anomalous threshold
therein. The online anomaly detection phase then uses the threshold to
identify anomalies in ‘windows’ of the test data.

_Performance metrics for anomaly detection_ : For assessing the performance of
the anomaly detection method, we establish the following metrics based on the
confusion matrix $\mathbf{M}$:

$\displaystyle\text{specificity}(\mathbf{M})$
$\displaystyle=\frac{\text{TN}}{\text{TN}+\text{FP}},$ (13)

where TN (True Negatives) represent the number of correctly identified normal
instances, while FP (False Positives) denote the misclassification of normal
instances as anomalies. For our experimental aim, high specificity is
desirable as it means the anomaly detection algorithm is effective at
correctly identifying and retaining “good” data, based on the learned
relationships between data, $\boldsymbol{y}$, collected at design,
$\boldsymbol{d}$. If specificity is low, the algorithm is mistakenly flagging
too many normal data points as anomalies, leading to the loss of valuable
information. This can be particularly problematic in scenarios where the cost
of False Positives is high, such as the travel and time costs associated with
checking remote in-situ sensors.

For completeness, we also consider the following performance metrics in the
assessment and comparison of the optimised designs,

$\displaystyle\text{sensitivity}(\mathbf{M})$
$\displaystyle=\frac{\text{TP}}{\text{TP}+\text{FN}},\text{ and}$
$\displaystyle\text{accuracy}(\mathbf{M})$
$\displaystyle=\frac{\text{TN}+\text{TP}}{\text{TN}+\text{TP}+\text{FP}+\text{FN}},$

where TP (True Positives) represent correctly classified anomalies and FN
(False Negatives) represent anomalies incorrectly classified as normal
instances. There is a trade-off between sensitivity and specificity.
Increasing specificity usually leads to a decrease in sensitivity.

Ultimately, we use the Matthews correlation coefficient [MCC 39] for an
overall evaluation as it provides an informative measure of classifier
performance, particularly useful in cases where the classes are imbalanced. A
larger value of the MCC indicates better performance of the classifier in
terms of both sensitivity and specificity, producing a score between -1 and
+1, defined by,

$\text{MCC}(\mathbf{M})=\frac{\text{TP}\times\text{TN}-\text{FP}\times\text{FN}}{\sqrt{(\text{TP}+\text{FP})(\text{TP}+\text{FN})(\text{TN}+\text{FP})(\text{TN}+\text{FN})}}.$
(14)

####  3.1.4 Dual-purpose utility

The utility function introduced herein is a compromise between the accuracy of
the posterior predictive mean at $\hat{n}$ specified locations (inverse root
mean squared error; irmse), and the specificity (sp) of anomaly detection at
$n$ sampled data locations. This compromise can be expressed as the product of
the two components,

$\displaystyle\text{U}({\boldsymbol{d}})$
$\displaystyle=\mathbb{E}[\text{u}_{\text{irmse}}({\boldsymbol{y}},{\boldsymbol{d}})]\times\mathbb{E}[\text{u}_{\text{sp}}(\mathbf{A}_{\text{true}},\hat{\mathbf{A}}_{\text{anom}})],$
(15)

with,

$\displaystyle\text{u}_{\text{irmse}}({\boldsymbol{y}},{\boldsymbol{d}})$
$\displaystyle=\bigg{(}\frac{1}{\hat{n}}\sum^{\hat{n}}_{i=1}(\mathbb{E}(\hat{y}_{i}|\boldsymbol{y},\hat{\boldsymbol{d}})-\hat{y}_{i})^{2}\bigg{)}^{-\frac{1}{2}},$
(16)
$\displaystyle\text{u}_{\text{sp}}(\mathbf{A}_{\text{true}},\hat{\mathbf{A}}_{\text{anom}})$
$\displaystyle=\text{specificity}(\mathbf{M}),\text{ for confusion matrix
$\mathbf{M}$,}$

where the $i$-th element of the posterior predictive mean is written as
$\mathbb{E}(\hat{y}_{i}|\boldsymbol{y},\hat{\boldsymbol{d}})$ and specificity
is defined in Equation (13). In this context, a design, $\boldsymbol{d}$,
represents the specification of physical sensor locations, while design,
$\hat{\boldsymbol{d}}$, denotes the design of prediction locations. Utilising
the inverse root mean squared error metric implies that design scenarios with
various $\hat{n}$ will result in utilities being measured on differing scales.
Note that multiple approaches exist to compromise between these objectives,
alternative methods may also be considered, including the incorporation of
other anomaly detection performance metrics.

###  3.2 Example implementation

An example implementation of our general approach is specified in Algorithm 3.
In line 1, we generate a true anomaly indicator matrix,
$\mathbf{A}_{\text{true}}$, specifying where actual anomalies occur in the
data collected at design sites, $\boldsymbol{y}$, based on the specified
anomaly generation method. Then, we simulate anomalies on the supposed data
collected, $\boldsymbol{y}_{\text{anom}}$, by adding normally distributed
noise based on the specified mean and standard deviation, in this case.
Training data is also generated from the likelihood function. The anomaly
detection algorithm is then applied to predict anomalies, given the design
locations, $\boldsymbol{d}$, anomalous data, $\boldsymbol{y}_{\text{anom}}$,
and training data, $\boldsymbol{y}_{\text{train}}$, resulting in a predicted
anomaly indicator matrix, $\hat{\mathbf{A}}_{\text{anom}}$. The success of the
anomaly detection is evaluated as $\text{u}_{\text{sp}}$, defined as
specificity computed from the confusion matrix. The predicted anomalies are
then imputed from the contaminated data $\boldsymbol{y}_{\text{anom}}$ to
create $\boldsymbol{y}_{\text{cleaned}}$. Crucially, it is
$\boldsymbol{y}_{\text{cleaned}}$ that is used to compute the posterior
predictive mean,
$\mathbb{E}(\hat{\boldsymbol{y}}|\hat{\boldsymbol{d}},\boldsymbol{y}_{\text{cleaned}})$
in the approximation to $\text{u}_{\text{irmse}}$ defined in Equation (16). In
the following section, we examine the errors associated with the posterior
predictive mean given $\boldsymbol{y}$, as well as
$\boldsymbol{y}_{\text{anom}}$, and $\boldsymbol{y}_{\text{cleaned}}$.
Finally, an approximation to the combined utility defined in Equation (15) is
computed as the product of the expected predictive errors (given potentially
anomalous data) and the expected specificity (performance of the anomaly
detection).

Input: $\boldsymbol{y}_{\text{train}}$ drawn from the likelihood function,
anomaly generation parameters corresponding to length, $\lambda_{A}$,
frequency, $p_{\text{A}}$, mean, $\mu_{\text{A}}$, and standard deviation,
$\sigma_{\text{A}}$.

Output: vector of utility samples.

1

2for  _$b\in 1:B$_ do

3 Draw $\boldsymbol{y}^{(b)}$ from the likelihood function.

4 Randomly generate indicator matrix,
$\mathbf{A}_{\text{true}}^{(b)}|\lambda_{A},p_{\text{A}}$, with entries:
$a_{ti}=\begin{cases}1,&\text{if element }(t,i)\text{ represents a true
anomaly},\\\ 0,&\text{otherwise},\end{cases}$ for $t=1,\ldots,T$ and
$i=1,\ldots,n$ per Section 3.1.2.

5 Randomly generate proposed contaminated data,
$\boldsymbol{y}^{(b)}_{\text{anom}}|\mu_{\text{A}},\sigma_{\text{A}}$, as:
$y_{\text{anom};ti}=\begin{cases}y_{ti},&\text{if }a_{ti}=0,\\\
y_{ti}+z,&\text{otherwise, with
}z\sim\mathcal{N}(\mu_{\text{A}},\sigma_{\text{A}}^{2})\end{cases}$ for
elements $a_{ti}\in\mathbf{A}_{\text{true}}^{(b)}$ and
$y_{ti}\in\boldsymbol{y}^{(b)}$.

6 Detect anomalies via a method outlined in Section 3.1.3; obtain predicted
anomaly indicator matrix $\hat{\mathbf{A}}_{\text{anom}}^{(b)}$.

7 Remove anomalies to obtain $\boldsymbol{y}^{(b)}_{\text{cleaned}}$ for which
the corresponding elements of $\hat{\mathbf{A}}_{\text{anom}}^{(b)}$ equal 0.

8

9

10Compute combined utility which can be approximated as,
${\text{U}}(\boldsymbol{d})\approx\hat{\text{U}}_{\text{irmse}}({\boldsymbol{d}})\times\hat{\text{U}}_{\text{sp}}(\mathbf{A}_{\text{true}},\hat{\mathbf{A}}_{\text{anom}}),$
via Monte Carlo approximation,
$\hat{\text{U}}_{\text{irmse}}({\boldsymbol{d}})=\frac{1}{B}\sum^{B}_{b=1}\text{u}_{\text{irmse}}({\boldsymbol{y}^{(b)}_{\text{cleaned}}},{\boldsymbol{d}}).$
and with
$\mathbf{A}_{\text{true}}=\odot^{B}_{b=1}\mathbf{A}^{(b)}_{\text{true}}$ and
$\hat{\mathbf{A}}_{\text{anom}}=\odot^{B}_{b=1}\hat{\mathbf{A}}_{\text{anom}}^{(b)}$
where $\odot$ represents vertical concatenation.

Algorithm 3 Implementation of Bayesian design with anomaly detection

##  4 Case studies

In this section, we investigate the proposed method on a spatial example and
then our motivating spatio-temporal river network example, each with a
corresponding anomaly detection method.

###  4.1 Spatial simulation

In this design problem, we consider $n=6$ spatial locations
$\boldsymbol{x}_{i}=(x_{i1},x_{i2})^{\top}$ within a specified 2-dimensional
area for observing responses $\boldsymbol{y}=(y_{1},...,y_{n})^{\top}$. A
Gaussian process model is fit to $\boldsymbol{y}$ for prediction at $\hat{n}$
predefined unobserved locations
$\hat{\boldsymbol{x}_{i}}=(\hat{x}_{i1},\hat{x}_{i2})^{\top}$, corresponding
to responses
$\hat{\boldsymbol{y}}=(\hat{y}_{1},...,\hat{y}_{\hat{n}})^{\top}$. The
combined vector of observed and unobserved responses is denoted as
$\tilde{\boldsymbol{y}}=(\boldsymbol{y}^{\top},\hat{\boldsymbol{y}}^{\top})^{\top}$,
where $\tilde{n}=n+\hat{n}$. The objective is to optimise the design
$\boldsymbol{d}:=\\{\boldsymbol{x}_{i}\\}_{i=1}^{n}$ for effective prediction
of $\hat{\boldsymbol{y}}$, despite anomalies in $\boldsymbol{y}$. The model is
a zero-mean Gaussian process, described by:

$p(\tilde{\boldsymbol{y}}|\sigma^{2},\rho,\tau^{2})\sim\mathcal{N}\big{(}\boldsymbol{0},\sigma^{2}\tilde{\boldsymbol{\Sigma}}\big{)}.$
(17)

Here,
$\tilde{\boldsymbol{\Sigma}}=\tilde{\mathbf{C}}+\tau^{2}\mathbf{I}_{\tilde{n}}$,
with $\tilde{n}\times\tilde{n}$ identity matrix $\mathbf{I}_{\tilde{n}}$, and
$\tilde{\mathbf{C}}$ being an $\tilde{n}\times\tilde{n}$ correlation matrix
partitioned as follows:

$\tilde{\mathbf{C}}=\begin{bmatrix}{\mathbf{C}}&{\mathbf{S}}\\\
{\mathbf{S}}^{\top}&{\bf\hat{C}}\end{bmatrix}.$

In this partition, ${\mathbf{C}}$ is the correlation matrix for design
locations $\boldsymbol{x}_{i}$, ${\bf\hat{C}}$ for prediction locations
$\hat{\boldsymbol{x}_{i}}$, and ${\mathbf{S}}$ between locations
$\boldsymbol{x}_{i}$ and $\hat{\boldsymbol{x}_{i}}$. The elements of
$\tilde{\mathbf{C}}$ are based on $C_{\text{Euc}}(h_{e};\sigma^{2},\rho)$ in
Equation (3), with priors:

$\displaystyle\sigma^{-2}$
$\displaystyle\sim\text{Gamma}\big{(}\frac{3}{2},\frac{1}{2}\big{)}$
$\displaystyle\rho$
$\displaystyle\sim\text{Uniform}\big{(}1,\frac{3}{2}\big{)}$

A small, fixed noise term $\tau=1\times 10^{-5}$ is assumed. The posterior
predictive mean is given by:

$\mathbb{E}(\hat{\boldsymbol{y}}|\hat{\boldsymbol{d}},\boldsymbol{y})={\mathbf{S}}^{\top}({\mathbf{C}}+\tau^{2}\mathbf{I}_{n})^{-1}\boldsymbol{y}.$
(18)

Consider a random starting configuration of the sensors, $\boldsymbol{d}_{1}$,
shown in Figure 1. Anomalies are generated according to Section 3.1.2 with
$p_{\text{A}}=0.10$, $\mu_{\text{A}}=5$, and $\sigma_{\text{A}}^{2}=10$. For
each design, training data with is drawn independently from the likelihood
function, in this example we set $B_{\text{train}}=100$. As per Algorithm 3,
the approximation to the expected utility is
$\hat{\text{U}}(\boldsymbol{d})=5.59$, corresponding to an anomaly detection
specificity of 96.4%. Data generated at $\boldsymbol{d}_{1}$ is shown in
Figure 2 across the first $B=100$ draws of the prior distribution, noting that
sensor 2 accounts for the largest false positive rate of 6% among the sensors.
Adjusting the position of sensor 3 to be nearer to sensor 2 effectively
reduces the occurrence of false positives in this scenario, shown in Figure 1
as $\boldsymbol{d}_{2}$. The design $\boldsymbol{d}_{2}$ reduces the rate of
false positives in Sensor 2, corresponding to an overall anomaly detection
specificity of 97.8% but decreases the utility approximation to
$\hat{\text{U}}(\boldsymbol{d})=4.98$ (due to an increase in false negatives).

Figure 1: Random design configuration, $\boldsymbol{d}_{1}$, and an adjusted
configuration, $\boldsymbol{d}_{2}$, where sensor 3 is relocated closer to
sensor 2. The background is coloured by a realisation of the Gaussian process
generated by Equation 17.  Figure 2: Variation in data sets $\boldsymbol{y}$,
$\boldsymbol{y}_{\text{anom}}$, and $\boldsymbol{y}_{\text{cleaned}}$
generated in Algorithm 3 for the spatial model, and the anomaly detection
method described in Algorithm 2. The sensor locations are specified by
$\boldsymbol{d}_{1}$ in Figure 1. Sensor 2 shows a false positive rate of 6%,
with an overall anomaly detection specificity of 96.4%.

We further examine the data variations $\boldsymbol{y}$,
$\boldsymbol{y}_{\text{anom}}$ and $\boldsymbol{y}_{\text{cleaned}}$ generated
in Algorithm 3, shown in Figure 2. Traditional optimal design procedures
operate under the assumption that the data collected is anomaly-free. This
perfect data is then used for making predictions, as indicated by the
posterior predictive mean,
$\mathbb{E}(\hat{\boldsymbol{y}}|\hat{\boldsymbol{d}},\boldsymbol{y})$.
However, practical data collection experiences suggest that this assumption is
often flawed. Figure 3 illustrates the posterior predictive mean for the
predefined $\hat{n}$ prediction locations, based on data of varying quality
collected at $\boldsymbol{d}_{1}$. The top part of the figure demonstrates the
scenario under the assumption of perfect data. The middle set of posterior
predictions,
$\mathbb{E}(\hat{\boldsymbol{y}}|\hat{\boldsymbol{d}},\boldsymbol{y}_{\text{anom}})$,
represents a real-life situation where anomalies in the measured data visibly
affect the posterior predictions. The bottom part of the figure displays the
posterior predictions,
$\mathbb{E}(\hat{\boldsymbol{y}}|\hat{\boldsymbol{d}},\boldsymbol{y}_{\text{cleaned}},\hat{\boldsymbol{d}})$,
following the automatic detection and removal of anomalies in Algorithm 3,
showcasing the potential improvements in data quality and prediction accuracy.

Figure 3: Illustration of posterior predictive means at specified $\hat{n}$
prediction locations using data collected at $\boldsymbol{d}_{1}$ under
different quality conditions. The top section reflects the assumption of
perfect data collection,
$\mathbb{E}(\hat{\boldsymbol{y}}|\hat{\boldsymbol{d}},\boldsymbol{y})$, the
middle shows the impact of anomalies in the data,
$\mathbb{E}(\hat{\boldsymbol{y}}|\hat{\boldsymbol{d}},\boldsymbol{y}_{\text{anom}})$,
and the bottom demonstrates the posterior predictions after automatic anomaly
removal,
$\mathbb{E}(\hat{\boldsymbol{y}}|\hat{\boldsymbol{d}},\boldsymbol{y}_{\text{cleaned}})$,
highlighting the influence of data quality on prediction accuracy.

|  |
---|---|---
|  |

(a) ACE trace plot for utility, u (proposed method). (b) ACE trace plot for
utility, $\text{u}_{\text{irmse}}$ (prediction errors only).

(c) ACE trace plot for utility, $\text{u}_{\text{sp}}$ (anomaly detection
performance only). (d) Prediction design, $\hat{\boldsymbol{d}}$, and the
optimal design, $\boldsymbol{d}_{*}$.

(e) Prediction design, $\hat{\boldsymbol{d}}$, the optimal design,
$\boldsymbol{d}_{*\text{irmse}}$. (f) Optimal design,
$\boldsymbol{d}_{*\text{sp}}$.

Figure 4: Top row: trace plots for the Approximate Coordinate Exchange (ACE)
algorithm, optimised according to their respective utilities, defined in
Equation 15. Coloured lines indicate different random starts. Bottom row:
spatial locations of corresponding optimal design points within the design
space, and the prediction design locations.

We ran the ACE algorithm with $N_{1}=30$, utilising the proposed dual-purpose
utility function from Algorithm 3, to find the optimal design
$\boldsymbol{d}_{*}$. The optimal design process outlined above was carried
out under various scenarios of anomaly contamination, as detailed in Table 1.
For brevity, we refer to non-anomalous data as “good” data. This table
displays the percentage of anomalies in each scenario, the corresponding
percentage of anomalies removed, and the percentage of good data retained. The
findings demonstrate successful removal of most anomalies (67.7%-80.4%) and
high retention of good data (94.0% - 99.3%), enhancing the overall quality of
the retained data set.

The trace plots for the ACE optimisation are displayed in Figure 4. As a
baseline method, we executed the optimisation algorithm outlined above for
minimising prediction error, $\text{u}_{\text{irmse}}$, resulting in the
optimal design, $\boldsymbol{d}_{*\text{irmse}}$. For comparative purposes, we
also optimised only for anomaly detection performance, $\text{u}_{\text{sp}}$,
resulting in optimal design, $\boldsymbol{d}_{*\text{sp}}$. The locations of
these optimal designs are depicted in Figure 4, with anomalies generated at a
proportion of 10% (Scenario 3). The design locations for
$\boldsymbol{d}_{*\text{irmse}}$ appear evenly spaced, consistent with other
patterns for a utility based on prediction accuracy. The observed clustering
of the design locations $\boldsymbol{d}_{*\text{sp}}$ is associated with the
anomaly detection process, which relies on neighbouring site values to compute
the expected mean. We also observe less uniform spacing in
$\boldsymbol{d}_{*}$ compared to the baseline optimal design,
$\boldsymbol{d}_{*\text{irmse}}$.

Table 1: Definition of anomaly detection scenarios defined by the proportion of generated anomalies, $p_{\text{A}}$. The expected utility, percentage of anomalies removed, and percentage of good data retained is reported for each scenario. Scenario | $p_{\text{A}}$ | % anomalies | $\boldsymbol{d}_{*}$
---|---|---|---
| $\hat{\text{U}}(\boldsymbol{d}_{*})$ |  % anomalies removed |  % good data retained
0 | 0.000 | 0.0 | 8.03 | - | 98.7
1 | 0.010 | 1.0 | 7.71 | 67.7 | 99.3
2 | 0.050 | 5.0 | 6.91 | 80.4 | 97.1
3 | 0.100 | 10.0 | 5.53 | 77.7 | 94.0

In Table 4.1, we assess the performance of our proposed optimal design,
$\boldsymbol{d}_{*}$, with the comparative optimal designs in terms of
prediction accuracy, $\text{u}_{\text{irmse}}$. In columns 1-3, we assume the
application of an automated anomaly detection procedure on the data at design
sites to produce $\boldsymbol{y}_{\text{cleaned}}$, as outlined in Algorithm
3. The last column in the table represents a practical scenario, where sensors
are placed according to the baseline optimal design for minimising prediction
errors, $\boldsymbol{d}_{*\text{irmse}}$, but the data collected at these
locations are anomalous, $\boldsymbol{y}_{\text{anom}}$.

As expected, the design optimised for prediction accuracy,
$\boldsymbol{d}_{*\text{irmse}}$, demonstrates the highest prediction
accuracy. The numbers in parenthesis indicate the standard deviation of the
utility values. The first row indicates the impact that removing False
Positives has on the error, since no anomalies are generated and
$\boldsymbol{y}_{\text{anom}}$ is considered perfect data (without anomalies).
In this case, the utility decreases marginally from 8.16 to 8.08. Scenario 1
indicates the shift away from the effectiveness of the baseline design, to the
effectiveness of an automatic anomaly detection procedure, at a contamination
rate of 1%. Interestingly, the proposed optimal design, $\boldsymbol{d}_{*}$,
demonstrates near-optimal prediction accuracy when anomalies are present. Note
that the anomaly detection design, $\boldsymbol{d}_{*\text{sp}}$, shows
considerably lower prediction accuracy, even compared to the last column with
no anomalies removed. This suggests that the predictive capability is largely
affected by more clustered sites. Particularly in the last row, a
contamination rate of 10% anomalies from data collected at the baseline
design, $\boldsymbol{d}_{*\text{irmse}}$, still results in lower error than
cleaned data collected at $\boldsymbol{d}_{*\text{sp}}$.

Table 2: Comparison of optimal designs $\boldsymbol{d}_{*}$,
$\boldsymbol{d}_{*\text{irmse}}$ and $\boldsymbol{d}_{*\text{sp}}$, under
differing levels of anomaly contamination. The left and middle columns
calculate the inverse root mean squared errors for each optimal design, using
an automatic anomaly detection process. The rightmost column illustrates a
real-world scenario with sensors placed according to
$\boldsymbol{d}_{*\text{irmse}}$ encountering anomalous data,
$\boldsymbol{y}_{\text{anom}}$.

Scenario |  $\text{u}_{\text{irmse}}(\boldsymbol{d}_{*})$ w/ |  |  |
---|---|---|---|---
$\text{u}_{\text{irmse}}(\boldsymbol{d}_{*\text{irmse}})$ w/ |  |  |
|  $\text{u}_{\text{irmse}}(\boldsymbol{d}_{*\text{sp}})$ w/ |  |  |
|  $\text{u}_{\text{irmse}}(\boldsymbol{d}_{*\text{irmse}})$ w/ |  |  |
| $\mathbb{E}(\hat{\boldsymbol{y}}|\hat{\boldsymbol{d}},\boldsymbol{y}_{\text{anom}})$ |  |  |
0 | 8.00 (4.73) |  8.08 (4.80) | 2.59 (1.66) |  8.16 (4.70)
1 |  7.81 (4.78) |  7.88 (4.80) | 2.57 (1.68) | 7.75 (4.92)
2 |  7.37 (4.52) |  7.43 (4.68) | 2.60 (1.77) | 6.44 (4.51)
3 |  5.84 (4.52) |  5.90 (4.62) | 2.73 (1.90) | 4.61 (5.14)

Table 3: Performance comparison of anomaly detection for the proposed optimal design, $\boldsymbol{d}_{*}$, the prediction-only optimal design, $\boldsymbol{d}_{*\text{irmse}}$, and the anomaly-only optimal design, $\boldsymbol{d}_{*\text{sp}}$, across various scenarios. Design | Scenario | Specificity | Accuracy | Sensitivity | MCC
---|---|---|---|---|---
$\boldsymbol{d}_{*}$ | 0 | 98.7 | 98.7 | - | -
$\boldsymbol{d}_{*\text{irmse}}$ | 0 | 98.9 | 98.9 | - | -
$\boldsymbol{d}_{*\text{sp}}$ | 0 | 99.7 | 99.7 | - | -
$\boldsymbol{d}_{*}$ | 1 | 99.3 | 99.0 | 67.7 | 57.5
$\boldsymbol{d}_{*\text{irmse}}$ | 1 | 97.9 | 97.7 | 81.8 | 47.5
$\boldsymbol{d}_{*\text{sp}}$ | 1 | 99.5 | 99.2 | 70.2 | 63.7
$\boldsymbol{d}_{*}$ | 2 | 97.1 | 96.2 | 80.4 | 66.7
$\boldsymbol{d}_{*\text{irmse}}$ | 2 | 95.3 | 94.6 | 83.0 | 60.7
$\boldsymbol{d}_{*\text{sp}}$ | 2 | 97.9 | 96.8 | 75.3 | 68.5
$\boldsymbol{d}_{*}$ | 3 | 94.0 | 92.3 | 77.7 | 63.6
$\boldsymbol{d}_{*\text{irmse}}$ | 3 | 92.6 | 91.3 | 80.5 | 61.9
$\boldsymbol{d}_{*\text{sp}}$ | 3 | 95.3 | 93.4 | 76.6 | 66.8

Table 3 illustrates the performance of anomaly detection among different
optimal designs: $\boldsymbol{d}_{*}$, $\boldsymbol{d}_{*\text{irmse}}$, and
$\boldsymbol{d}_{*\text{sp}}$ across various scenarios. Results in the table
were computed with $B=15,000$ Monte Carlo draws, resulting in $B\times
n=90,000$ binary classifications. In scenarios where anomalies are present, we
observe that the design optimised for specificity has the highest specificity,
as expected. The specificity is generally high across all designs in scenarios
0, 1 & 2 due to a proportionally high number of True Negatives, given the
class imbalance with an anomaly rate of 0%, 1% & 5%, respectively. For a more
balanced evaluation of the designs, we also consider the Matthews Correlation
Coefficient (MCC), defined in Equation (14). We observe that
$\boldsymbol{d}_{*}$ consistently outperforms
$\boldsymbol{d}_{*\text{irmse}}$, with a gain ranging between 1.4 - 10.5% in
the MCC. The comparison study indicates a trade-off between predictive
accuracy and anomaly detection performance in the dual-purpose utility,
$\boldsymbol{d}_{*}$, as expected.

###  4.2 Spatio-temporal river simulation

In this case study, we explore the spatio-temporal river network model
outlined in Section 2.1.2, simulating a river network with 300 segments with
the R software package SSN [40]. The aim of this design scenario is to
maximise prediction accuracy at $\hat{n}=56$ predetermined locations denoted
by design, $\hat{\boldsymbol{d}}$, while considering integrity of the data.
Recall that $\boldsymbol{y}$ is the stacked vector of observations from the
design locations over time, $T$, and denote $\hat{\boldsymbol{y}}$ as the
stacked vector of predictions predictions across the network over time. The
matrices $\mathbf{X}$ and $\hat{\mathbf{X}}$ represent the space-time design
matrices for covariates, and $\boldsymbol{\beta}$ is the vector of regression
coefficients. Using the simple kriging approach, the posterior predictive mean
is given by,

$\mathbb{E}(\hat{\boldsymbol{y}}|\hat{\boldsymbol{d}},\boldsymbol{y})=\hat{\mathbf{X}}\boldsymbol{\beta}+{\mathbf{S}}^{\top}\mathbf{C}^{-1}(\boldsymbol{y}-\mathbf{X}\boldsymbol{\beta}),$
(19)

where the covariance matrix $\mathbf{S}$, of dimension $n\times T$ by
$\hat{n}\times T$, denotes the covariance between the observation and
prediction points in space and time. We assumed the tail-up exponential model
as per Equation 6, with the following priors:

$\displaystyle\sigma_{0}^{2}$ $\displaystyle\sim\text{Uniform}(1\times
10^{-5},1\times 10^{-4})$ $\displaystyle\sigma_{u}^{2}$
$\displaystyle\sim\text{Uniform}\left(1,\frac{3}{2}\right)$
$\displaystyle\alpha_{u}$
$\displaystyle\sim\text{Uniform}\left(2,\frac{5}{2}\right)$
$\displaystyle\phi$ $\displaystyle\sim\text{Uniform}(0,1)$

This design problem positions $n=14$ sensors within the river network,
$\mathcal{S}$, collecting data at $T=50$ distinct time intervals. Figure 5
illustrates randomly generated design, denoted as $\boldsymbol{d}_{1}$, and
the prediction design, $\hat{\boldsymbol{d}}$. We refer to the $i$-th
coordinate of design $\boldsymbol{d}$ as $\boldsymbol{d}(i)$. The intricate
nature of the river network’s spatial domain is evident, with the branching
network topology embedded in a 3D terrestrial landscape. To accommodate for
this complexity, modifications were essential in the deployment of the ACE
methodology. For this, we randomly assigned the design space for design
coordinate, $\boldsymbol{d}(i)$, to a corresponding network ‘path’. A network
path is a sequence of segments spanning from the most upstream reach of a
river network to the outlet. In this way, the design is governed by a singular
parameter, distance upstream. Refer to Figure 5 for a visual representation of
a path within the river network.

|
---|---

(a) Prediction locations, $\hat{\boldsymbol{s}}$, of design,
$\hat{\boldsymbol{d}}$, across the river network (b) Random start design,
$\boldsymbol{d}_{1}$, and the associated design space for sensor $i$ within
the design.

Figure 5: Visualisation of design configurations across a dendritic river
system, comprising 300 randomly generated segments. The watercourse flows from
the upper reaches toward the lower basin, converging at the river outlet
positioned at coordinates (0,0).

Consider the data variants generated in Algorithm 3. The Oddstream anomaly
detection algorithm was run with a window length of $w=50$, chosen as a hyper
parameter. Training data were independently generated from the likelihood
function based on the given design. Accordingly, each training set was
comprised of $T\times n$ observations. For computational efficiency, anomalies
were detected for all $B$ Monte Carlo draws at once. To ensure this detection
remained independent for each Monte Carlo draw, the window length, $w$, was
set to $T$. Recall that Oddstream yields an output series $(a_{1},...,a_{n})$,
where $a_{i}=1$ indicates an anomaly detected at the $i$-th sensor (otherwise,
$a_{i}=0$), for a given input window. As such, this approach groups anomaly
detection by windows of size $w$, making it effective for very large datasets
and persistent anomaly patterns. The predicted anomalies are then imputed with
the posterior predictive mean,
$\mathbb{E}(\hat{y}_{i}|\boldsymbol{y},\hat{\boldsymbol{d}})$, for the window.
See Figure 6 for an illustration of the anomaly detection method and the
resulting data variants.

Figure 6: Data variants, $\boldsymbol{y}$, $\boldsymbol{y}_{\text{anom}}$, and
$\boldsymbol{y}_{\text{cleaned}}$ generated in Algorithm 3, for the spatio-
temporal river network model, and the anomaly detection method described in
Algorithm 4 (see Appendix). The parameter $w$ represents the window length
utilised in the anomaly detection method.

ACE was run to find optimal designs; see trace plots in Figure 7. At the
random start, $k$, each design coordinate is assigned to a fixed ‘path’. The
number of random starts was set to $K=10$. When employing an anomaly detection
method that utilises windowing, the assessment of performance becomes grouped
by the length of the window, reducing the number of individual data points
available for evaluating performance by a factor of $w$. Hence, we may observe
increased fluctuations in the trace plot, prompting us to raise the number of
iterations to $N_{1}=40$. The locations of the optimal designs are shown in
Figure 7, highlighting a similar pattern to the previous case study, even
within a more complex spatial domain. In $\boldsymbol{d}_{*\text{sp}}$,
optimising for anomaly detection specificity reveals clustered optimal design
locations, particularly around the river outlet. Conversely,
$\boldsymbol{d}_{*\text{irmse}}$ exhibits a more uniform distribution across
the river network. The proposed method’s $\boldsymbol{d}_{*}$ shows a mix of
clustered sites around the river outlet and evenly dispersed sites across the
river network.

|  |
---|---|---
|  |

(a) ACE trace plot for utility, u (proposed method). (b) ACE trace plot for
utility, $\text{u}_{\text{irmse}}$ (prediction errors only).

(c) ACE trace plot for utility, $\text{u}_{\text{sp}}$ (anomaly detection
performance only). (d) Locations given by optimal design,
$\boldsymbol{d}_{*}$, across the river network.

(e) Locations given by optimal design, $\boldsymbol{d}_{*\text{irmse}}$,
across the river network. (f) Locations given by optimal design,
$\boldsymbol{d}_{*\text{sp}}$, across the river network.

Figure 7: Top row: trace plots for the Approximate Coordinate Exchange (ACE) algorithm, optimised according to their respective utilities, defined in Equation 15. Bottom row: locations of corresponding design points across the river network, for the spatio-temporal model. Table 4: Scenarios of various levels of anomaly generation in a spatio-temporal context, defined by the anomaly proportion, $p_{\text{A}}$, and the anomaly length, $\lambda_{\text{A}}$, for in the anomaly generation methodology in Section 3.1.2. Corresponding average percentages of generated anomalies for each scenario are shown on the right. Scenario | $p_{\text{A}}$ | $\lambda_{\text{A}}$ | % anomalies
---|---|---|---
0 | 0.000 | 0.00 | 0.0
1 | 0.001 | 1.25 | 0.6
2 | 0.010 | 1.25 | 6.5
3 | 0.010 | 1.75 | 16.5
Table 5: Summary of performance metrics for the optimal designs, $\boldsymbol{d}_{*}$, across different scenarios. Displayed are the expected utility, $\hat{\text{U}}(\boldsymbol{d}_{*})$, the proportion of anomalies automatically removed, and the proportion of good data retained. These results illustrate the balance between anomaly mitigation and data preservation achieved by the design across various levels of data contamination. Scenario | $\boldsymbol{d}_{*}$
---|---
$\hat{\text{U}}(\boldsymbol{d}_{*})$ |  % anomalies removed |  % good data retained
0 | 9.21 | - | 99.1
1 | 8.44 | 91.2 | 99.1
2 | 8.94 | 93.4 | 98.2
3 | 9.65 | 96.5 | 97.8

The methodology described above was executed across the scenarios detailed in
Table 4, representing varying anomaly generation levels in a spatio-temporal
context. Table 5 presents the proposed expected utility,
$\hat{\text{U}}(\boldsymbol{d}_{*})$, the percentage of automatically removed
(true) anomalies, and the percentage of retained good data. Notably, across
all scenarios (0 to 3), both the percentage of anomalies automatically
detected and the percentage of retained good data are notably high.

For instance, as the percentage of automatically removed anomalies increases
from 91.2% to 96.%, the percentage of preserved good data remains high, over
98% in all cases. The findings demonstrate the ability of the optimal design
to effectively identify anomalies while retaining a significant volume of good
data, thereby enhancing data quality and showcasing resilience against varying
contamination levels.

Table 6: Efficacy of the proposed optimal design, $\boldsymbol{d}_{*}$, versus
the comparative designs, $\boldsymbol{d}_{*\text{irmse}}$ and
$\boldsymbol{d}_{*\text{sp}}$, in terms of prediction error across various
scenarios. Data collected at each design are assumed to be either cleaned from
detected anomalies, $\boldsymbol{y}_{\text{cleaned}}$, or anomalous,
$\boldsymbol{y}_{\text{anom}}$.

In river network in-situ data, various types of anomalies can occur,
indicating abnormal patterns or irregularities in the data. These anomalies
can include sudden spikes or drops in water temperature, unexpected changes in
flow rates, unusual levels of pollutants, or irregularities in sediment
deposition. Additionally, anomalies might arise from sensor malfunctions, data
transmission errors, or human interference. Detecting these anomalies is
crucial for maintaining the accuracy and reliability of the data collected
from river network sensors, enabling timely responses to environmental
changes, pollution incidents, or equipment failures. Automated anomaly
detection methods are essential for identifying these irregularities promptly
and ensuring the integrity of the collected data.

The case study analysis highlighted the distinct strengths among the optimal
design configurations: $\boldsymbol{d}_{*}$ (proposed optimal design),
$\boldsymbol{d}_{*\text{irmse}}$ (prediction-focused optimal design), and
$\boldsymbol{d}_{*\text{sp}}$ (optimal design accounting for anomalies). In
both case studies the optimal design, $\boldsymbol{d}_{*}$, demonstrated high
predictive accuracy similar to $\boldsymbol{d}_{*\text{irmse}}$. Conversely,
the comparable optimal design, $\boldsymbol{d}_{*\text{sp}}$, demonstrated
high anomaly detection performance but revealed limited predictive ability,
especially with increased anomalies. The proposed design $\boldsymbol{d}_{*}$
significantly improved anomaly detection performance compared to the baseline
design $\boldsymbol{d}_{*\text{irmse}}$. Overall, we found that
$\boldsymbol{d}_{*}$ offers a balanced performance in both prediction and
anomaly detection.

In the context of anomaly detection for data quality assurance, high
specificity means that the anomaly detection algorithm is highly effective at
correctly identifying non-anomalous data points as normal. This is important
because it minimises the risk of mistakenly discarding or altering normal
data, amidst the complexity of potential anomalies in data collected from in-
situ sensors, which may have diverse origins. However, overemphasising
specificity can lead to a model that misses anomalies in the data, potentially
missing critical anomalies in the process. In cases like extreme weather
events, where anomaly certainty is crucial, metrics like accuracy may be more
relevant. In our context, the focus on specificity helped to filter out some
outliers, but it’s essential to consider other anomaly detection metrics
depending on the scenario. Future work may also extend the approach to allow
for the potential of missing data within the design.

While our focus has been on BOED for prediction and anomaly detection for
river network sensor data, there are broader implications and applications of
this novel method. Spatial and spatio-temporal models are extensively applied
across a wide range of research fields, from economics to environmental
science, urban planning, epidemiology, and meteorology. Optimal design is
crucial in collecting and analysing data, providing insights into complex
dynamics and interactions in space and time. Intelligent design with anomaly
detection is an exciting new concept that can provide invaluable tools for
researchers and policymakers more broadly, from tracking the spread of
diseases, planning more effective urban development strategies, to better
predicting extreme weather patterns.

##  6 Acknowledgments

This project was supported by the Australian Research Council (ARC) Linkage
Project (LP180101151) “Revolutionising water-quality monitoring in the
information age”.

##  7 Data and Code Accessibility

The design optimisations were performed using the R package acebayes [33].
This code implementation can be found at
https://github.com/KatieBuc/design_anom.

##  Appendix A Extended Data

###  A.1 Oddstream Anomaly Detection

// Training Phase

Input: Collection of time series data,
$\boldsymbol{y}_{\text{train}}\in\mathbb{R}^{l_{\text{train}}\times n}$,
representative of the typical system behavior.

Output: Anomalous threshold, $\tau^{*}$.

1 Extract $m$ features from each time series in
$\boldsymbol{y}_{\text{train}}$ and create a feature matrix $\mathbf{M}$ of
size $n\times m$.

2 Normalise the columns of $\mathbf{M}$ to produce $\mathbf{M}^{*}$.

3 Apply principal component analysis to $\mathbf{M}^{*}$.

4 Create a 2D space using the first two principal components; each data point
in this space corresponds to a time series in $\boldsymbol{y}_{\text{train}}$.

5 Estimate the probability density, $\hat{f}_{2}$, of the 2D space using
kernel density estimation with a bivariate Gaussian kernel.

6 Draw a large number of extremes from the estimated probability density
function, $\hat{f}_{2}$, and form an empirical distribution in the
$\Psi$-transform space, with
$\Psi|\hat{f}_{2}(x)|=\begin{cases}(-2ln(\hat{f}_{2}(x))-2ln(2\pi))^{\frac{1}{2}},&\hat{f}_{2}(x)<(2\pi)^{-1}\\\
0,&\text{otherwise}.\end{cases}$

7 Fit a Gumbel distribution to the $\Psi$-transformed values.

8 Determine the anomalous threshold contour, $\tau^{*}$, (i.e. where the most
extreme of the $n$ typical samples generated from $f_{2}$ will lie to some
given level of probability, $\tau_{F}$) by equating the CDF of the fitted
Gumbel distribution to $\tau_{F}$.

9

// Online Anomaly Detection Phase

Input: Matrix of data as a sliding window
$\boldsymbol{y}_{\text{anom},[j-w,j]}\in\mathbb{R}^{w\times n}$, anomalous
threshold $\tau^{*}$.

Output: Indicator vector, $\boldsymbol{a}_{j}$, of anomalous series in
$\boldsymbol{y}_{\text{anom},[j-w,j]}$.

10 Extract $m$ features from the time series
$\boldsymbol{y}_{\text{anom},[j-w,j]}$ to create feature matrix
$\mathbf{M}_{\text{test}}$.

11 Project $\mathbf{M}_{\text{test}}$ onto the 2D PC space obtained from
typical data.

12 Calculate probability densities for the projected data points with respect
to $\hat{f}_{2}$ (from the training phase).

13 Identify time series where $\hat{f}_{2}(y_{i})<\tau^{*}$ for
$i=1,2,\ldots,n$ and return $\boldsymbol{a}_{j}=(a_{j1},...,a_{jn})$ with
$a_{ji}=1$ for an anomaly in the $i$-th sensor (else $a_{ji}=0$), given window
input $\boldsymbol{y}_{\text{anom},[j-w,j]}$.

14 Repeat steps 1–4 for each new time window.

Algorithm 4 Oddstream: Anomaly Detection in Spatio-temporal Data

###  A.2 Approximate coordinate exchange

We specify a design, $\boldsymbol{d}\in\mathcal{D}\subset\mathbb{R}^{n\times
k}$. Phase I of the algorithm utilises cyclic ascent to maximise the
approximation $\hat{\text{U}}(\boldsymbol{d})$ to the expected utility. For
each coordinate, $x_{ij}$, a one-dimensional emulator of
$\tilde{\text{U}}(\boldsymbol{d})$ is constructed as the mean of the posterior
predictive distribution conditioned on a few evaluations of
$\hat{\text{U}}(\boldsymbol{d})$, assuming a Gaussian process (GP) prior. The
emulator for the $ij$-th coordinate is constructed as follows:

  1. 1.

Select a one-dimensional space-filling design with points,
$x^{q}_{ij}\in\mathcal{D}_{ij}$, for $q=1,\ldots,Q$.

  2. 2.

Construct designs, $\boldsymbol{d}^{q}_{ij}$, identical to the current design
but with the $ij$-th coordinate replaced by a point in the one-dimensional
space-filling design, x^q_i = (x_i1, …, x_ij-1, x^q_ij, x_ij+1, …, x_ik)^⊤.

  3. 3.

Evaluate $\tilde{\text{U}}(\boldsymbol{d}^{q}_{ij})$ for $q=1,\ldots,Q$.

  4. 4.

Fit a GP regression model to the data
$\\{\boldsymbol{x}^{q}_{i},\tilde{\text{U}}(\boldsymbol{d}^{q}_{ij})\\}_{q=1}^{Q}$
and set emulator $\hat{\text{U}}_{ij}(x)$ as the resulting predictive mean.

  5. 5.

Find x^*_ij = arg max_x ∈D_ij^U_ij(x) and produce design
$\boldsymbol{d}^{*}_{ij}$ with $i$-th row
$(\boldsymbol{x}^{*}_{i})^{\top}=(x_{i1},\ldots,x_{ij-1},x^{*}_{ij},x_{ij+1},\ldots,x_{ik})^{\top}$.

The approach to maximising $\hat{\text{U}}_{ij}(x)$ is robust to multi-modal
emulators and is efficient due to the minimal computational overhead in
evaluating the predictive mean.

## References

  *   * Bon et al. [2023] Bon, J.J., Bretherton, A., Buchhorn, K., Cramb, S., Drovandi, C., Hassan, C., Jenner, A.L., Mayfield, H.J., McGree, J.M., Mengersen, K., et al.: Being Bayesian in the 2020s: opportunities and challenges in the practice of modern applied Bayesian statistics. Philosophical Transactions of the Royal Society A 381(2247), 20220156 (2023)
  * Kuhfeld et al. [1994] Kuhfeld, W.F., Tobias, R.D., Garratt, M.: Efficient experimental design with marketing research applications. Journal of Marketing Research 31(4), 545–557 (1994)
  * Zhang et al. [2018] Zhang, J.F., Papanikolaou, N.E., Kypraios, T., Drovandi, C.C.: Optimal experimental design for predator–prey functional response experiments. Journal of The Royal Society Interface 15(144), 20180186 (2018)
  * Myung et al. [2013] Myung, J.I., Cavagnaro, D.R., Pitt, M.A.: A tutorial on adaptive design optimization. Journal of Mathematical Psychology 57(3-4), 53–67 (2013)
  * Watson [2017] Watson, A.B.: Quest+: A general multidimensional Bayesian adaptive psychometric method. Journal of Vision 17(3), 10–10 (2017)
  * Huan and Marzouk [2013] Huan, X., Marzouk, Y.M.: Simulation-based optimal Bayesian experimental design for nonlinear systems. Journal of Computational Physics 232(1), 288–317 (2013)
  * Loredo [2004] Loredo, T.J.: Bayesian adaptive exploration. In: AIP Conference Proceedings, vol. 707, pp. 330–346 (2004). American Institute of Physics
  * Cheng and Shen [2005] Cheng, Y., Shen, Y.: Bayesian adaptive designs for clinical trials. Biometrika 92(3), 633–646 (2005)
  * Bolstad et al. [2004] Bolstad, B., Collin, F., Simpson, K., Irizarry, R., Speed, T.: Experimental design and low-level analysis of microarray data. International Review of Neurobiology 60, 25–58 (2004)
  * Tsou [2010] Tsou, J.-C.: Production system with process quality control: modelling and application. International Journal of Systems Science 41(7), 865–874 (2010)
  * Edwards [2010] Edwards, P.: Questionnaires in clinical trials: guidelines for optimal design and administration. Trials 11, 1–8 (2010)
  * Tourangeau et al. [2004] Tourangeau, R., Couper, M.P., Conrad, F.: Spacing, position, and order: Interpretive heuristics for visual features of survey questions. Public Opinion Quarterly 68(3), 368–393 (2004)
  * Fuentes et al. [2007] Fuentes, M., Chaudhuri, A., Holland, D.M.: Bayesian entropy for spatial sampling design of environmental data. Environmental and Ecological Statistics 14, 323–340 (2007)
  * Thilan et al. [2023] Thilan, A., Menéndez, P., McGree, J.: Assessing the ability of adaptive designs to capture trends in hard coral cover. Environmetrics, 2802 (2023)
  * Abeysiri Wickrama Liyanaarachchige et al. [2022] Abeysiri Wickrama Liyanaarachchige, P.T., Fisher, R., Thompson, H., Menendez, P., Gilmour, J., McGree, J.M.: Adaptive monitoring of coral health at scott reef where data exhibit nonlinear and disturbed trends over time. Ecology and Evolution 12(9), 9233 (2022)
  * Scarpa et al. [2007] Scarpa, R., Campbell, D., Hutchinson, W.G.: Benefit estimates for landscape improvements: sequential Bayesian design and respondents’ rationality in a choice experiment. Land Economics 83(4), 617–634 (2007)
  * Buxton et al. [1999] Buxton, A., Livingston, A.G., Pistikopoulos, E.N.: Optimal design of solvent blends for environmental impact minimization. AIChE Journal 45(4), 817–843 (1999)
  * Nikolopoulou and Ierapetritou [2012] Nikolopoulou, A., Ierapetritou, M.G.: Optimal design of sustainable chemical processes and supply chains: A review. Computers & Chemical Engineering 44, 94–103 (2012)
  * Rainforth et al. [2023] Rainforth, T., Foster, A., Ivanova, D.R., Smith, F.B.: Modern Bayesian experimental design. arXiv preprint arXiv:2302.14545 (2023)
  * Beck et al. [2018] Beck, J., Dia, B.M., Espath, L.F., Long, Q., Tempone, R.: Fast Bayesian experimental design: Laplace-based importance sampling for the expected information gain. Computer Methods in Applied Mechanics and Engineering 334, 523–553 (2018)
  * Peterson et al. [2013] Peterson, E.E., Ver Hoef, J.M., Isaak, D.J., Falke, J.A., Fortin, M.-J., Jordan, C.E., McNyset, K., Monestiez, P., Ruesch, A.S., Sengupta, A., et al.: Modelling dendritic ecological networks in space: an integrated network perspective. Ecology Letters 16(5), 707–719 (2013)
  * Buchhorn et al. [2023] Buchhorn, K., Mengersen, K., Santos-Fernandez, E., Peterson, E.E., McGree, J.M.: Bayesian design with sampling windows for complex spatial processes. Journal of the Royal Statistical Society Series C: Applied Statistics, 099 (2023) https://doi.org/10.1093/jrsssc/qlad099
  * Leigh et al. [2019] Leigh, C., Alsibai, O., Hyndman, R.J., Kandanaarachchi, S., King, O.C., McGree, J.M., Neelamraju, C., Strauss, J., Talagala, P.D., Turner, R.D., et al.: A framework for automated anomaly detection in high frequency water-quality data from in situ sensors. Science of the Total Environment 664, 885–898 (2019)
  * Rodriguez-Perez et al. [2020] Rodriguez-Perez, J., Leigh, C., Liquet, B., Kermorvant, C., Peterson, E., Sous, D., Mengersen, K.: Detecting technical anomalies in high-frequency water-quality data using artificial neural networks. Environmental Science & Technology 54(21), 13719–13730 (2020)
  * Jones et al. [2022] Jones, A.S., Jones, T.L., Horsburgh, J.S.: Toward automating post processing of aquatic sensor data. Environmental Modelling & Software 151, 105364 (2022)
  * Santos-Fernandez et al. [2023] Santos-Fernandez, E., Ver Hoef, J., Peterson, E., McGree, J., Villa, C., Leigh, C., Turner, R., Roberts, C., Mengersen, K.: Unsupervised anomaly detection in spatio-temporal stream network sensor data (2023) https://doi.org/10.13140/RG.2.2.33200.74241
  * Buchhorn et al. [2023] Buchhorn, K., Santos-Fernandez, E., Mengersen, K., Salomone, R.: Graph neural network-based anomaly detection for river network systems. F1000Research 12(991) (2023) https://doi.org/10.12688/f1000research.136097.1
  * Talagala et al. [2020] Talagala, P.D., Hyndman, R.J., Smith-Miles, K., Kandanaarachchi, S., Munoz, M.A.: Anomaly detection in streaming nonstationary temporal data. Journal of Computational and Graphical Statistics 29(1), 13–27 (2020)
  * Lindley [1956] Lindley, D.V.: On a measure of the information provided by an experiment. The Annals of Mathematical Statistics 27(4), 986–1005 (1956)
  * Stone [1959] Stone, M.: Application of a measure of information to the design and comparison of regression experiments. The Annals of Mathematical Statistics, 55–70 (1959)
  * Ver Hoef and Peterson [2010] Ver Hoef, J.M., Peterson, E.E.: A moving average approach for spatial statistical models of stream networks. Journal of the American Statistical Association 105(489), 6–18 (2010)
  * Santos-Fernandez et al. [2022] Santos-Fernandez, E., Ver Hoef, J.M., Peterson, E.E., McGree, J., Isaak, D.J., Mengersen, K.: Bayesian spatio-temporal models for stream networks. Computational Statistics & Data Analysis 170, 107446 (2022)
  * Overstall et al. [2017] Overstall, A., Woods, D., Adamou, M.: acebayes: An R package for Bayesian optimal design of experiments via approximate coordinate exchange. arXiv preprint arXiv:1705.08096 (2017)
  * Ahmed et al. [2016] Ahmed, M., Mahmood, A.N., Hu, J.: A survey of network anomaly detection techniques. Journal of Network and Computer Applications 60, 19–31 (2016)
  * Nassif et al. [2021] Nassif, A.B., Talib, M.A., Nasir, Q., Dakalbab, F.M.: Machine learning for anomaly detection: A systematic review. Ieee Access 9, 78658–78700 (2021)
  * Pang et al. [2021] Pang, G., Shen, C., Cao, L., Hengel, A.V.D.: Deep learning for anomaly detection: A review. ACM computing surveys (CSUR) 54(2), 1–38 (2021)
  * Chandola et al. [2009] Chandola, V., Banerjee, A., Kumar, V.: Anomaly detection: A survey. ACM computing surveys (CSUR) 41(3), 1–58 (2009)
  * Talagala et al. [2019] Talagala, P.D., Hyndman, R.J., Leigh, C., Mengersen, K., Smith-Miles, K.: A feature-based procedure for detecting technical outliers in water-quality data from in situ sensors. Water Resources Research 55(11), 8547–8568 (2019)
  * Matthews [1975] Matthews, B.W.: Comparison of the predicted and observed secondary structure of t4 phage lysozyme. Biochimica et Biophysica Acta (BBA)-Protein Structure 405(2), 442–451 (1975)
  * Ver Hoef et al. [2014] Ver Hoef, J., Peterson, E., Clifford, D., Shah, R.: SSN: An R package for spatial statistical modeling on stream networks. Journal of Statistical Software 56, 1–45 (2014)
  * Kang et al. [2016] Kang, S.Y., McGree, J.M., Drovandi, C.C., Caley, M.J., Mengersen, K.L.: Bayesian adaptive design: improving the effectiveness of monitoring of the Great Barrier Reef. Ecological Applications 26(8), 2637–2648 (2016)
  * Armour et al. [2009] Armour, J., Hateley, L., Pitt, G.: Catchment modelling of sediment, nitrogen and phosphorus nutrient loads with SedNet/ANNEX in the Tully–Murray basin. Marine and Freshwater Research 60(11), 1091–1096 (2009)