# Bayes Linear Analysis for Statistical Modelling with Uncertain Inputs

Samuel E<EMAIL_ADDRESS>David C. Woods

Durham University, Durham, UK

###### Abstract

Statistical models typically capture uncertainties in our knowledge of the
corresponding real-world processes, however, it is less common for this
uncertainty specification to capture uncertainty surrounding the values of the
inputs to the model, which are often assumed known. We develop general
modelling methodology with uncertain inputs in the context of the Bayes linear
paradigm, which involves adjustment of second-order belief specifications over
all quantities of interest only, without the requirement for probabilistic
specifications. In particular, we propose an extension of commonly-employed
second-order modelling assumptions to the case of uncertain inputs, with
explicit implementation in the context of regression analysis, stochastic
process modelling, and statistical emulation. We apply the methodology to a
regression model for extracting aluminium by electrolysis, and emulation of
the motivating epidemiological simulator chain to model the impact of an
airborne infectious disease.

##  1 Introduction

Most often, it is assumed that the input, or independent, variables of a
statistical model are known, and that uncertainties largely result from
measurement error in the response and an imperfect description of the link
between the independent variables and the responses. However, in reality, the
inputs may have uncertainty surrounding their values, that is, they may
themselves be random variables. This may arise as a result of imperfect
measurements of the corresponding quantities in the physical experiment, or
when those input variables may have been themselves modelled as an independent
process. Such modelling is sometimes referred to as Error in Variables models
(Durbin, 1954; Figueroa-Zuniga et al., 2022), although such models are
primarily concerned with estimation at specified independent variable values
when the independent variables for the training data have been measured with
error (often a parameterised homoscedastic error structure to be estimated).

In this paper, we develop general Bayes linear modelling methods for
experiments where there is specified uncertainty in the values of the inputs
for both the training data used to estimate the model and use-case data for
which predictions are required. Following de Finetti (de Finetti, 1974, 1975;
Whittle, 1992), Bayes linear methods (Goldstein & Wooff, 2007) involve prior
second-order belief specification over all quantities of interest; the beliefs
about unobserved quantities then being adjusted in light of those which have
been observed. A general advantage of the Bayes linear approach over the fully
Bayesian approach is the lack of requirement to specify full probabilistic
distributions over all quantities of interest, these often being difficult to
specify meaningfully and thus often chosen for computational convenience.
Bayes linear methods have seen diverse application, including in the
petrochemical (Craig et al., 1996), medical (Gosling et al., 2013) and climate
(Astfalck et al., 2021) sciences.

We propose an extension of second-order modelling assumptions to the case of
uncertain inputs. Whilst the methodology is generally applicable, two
important scenarios will be studied in detail. Firstly, uncertain input
modelling for regression analyses, applied to a specific example taken from
Goldstein & Wooff (1998). In this example, input uncertainty is assumed to
arise from inaccurate physical measurements training input values. Secondly,
we address emulation of computer models, or simulators, using a Bayes linear
framework, often termed Bayes linear emulation (Cumming & Goldstein, 2009;
Goldstein & Huntley, 2016). We assume uncertain inputs arising, for example,
as the outputs of a stochastic process (perhaps another emulated simulator, as
is discussed in Section 5.3) about which we are uncertain. Such application
first requires consideration of stochastic process modelling.

The article is structured as follows. In Section 2, we formally introduce
concepts and notation concerning Bayes linear methods, before discussing our
developed general framework to statistical modelling with uncertain inputs
within a Bayes linear paradigm. In Section 3, we demonstrate the Bayes linear
modelling approach in the context of (exchangeable) regression. In Section 4,
we develop the methodology to model stochastic processes. In Section 5, we
combine the scenarios presented in Sections 3 and 4 to develop Uncertain Input
Bayes Linear Emulation and demonstrate its application to emulating a chain of
linked simulators. The methods in this section are implemented in R, with the
developed R packages available at https://github.com/Jackson-SE/UIBLE. Section
6 contains a brief discussion and some directions for future research.

##  2 Bayes Linear Analysis and Uncertain Inputs

In this article, we focus on the Bayes Linear approach (Hartigan, 1969;
O’Hagan, 1987; Goldstein & Wooff, 2007) to statistical inference, which deals
with second-order belief specifications (that is, expectations, variances and
covariances) of observable quantities. Probabilities can be represented as the
expectation of the corresponding indicator functions when required.

More precisely, suppose that there are two collections of random quantities,
$\mathcal{B}=(B_{1},...,B_{r})$ and $\mathcal{D}=(D_{1},...,D_{s})$. Bayes
linear analysis involves updating subjective beliefs about $\mathcal{B}$ given
observation of $\mathcal{D}$. In order to do so, prior mean vectors and
covariance matrices for $\mathcal{B}$ and $\mathcal{D}$ (that is,
$\textrm{E}[\mathcal{B}]$, $\textrm{E}[\mathcal{D}]$,
$\textrm{Var}[\mathcal{B}]$ and $\textrm{Var}[\mathcal{D}]$), along with a
covariance matrix between $\mathcal{B}$ and $\mathcal{D}$ (that is,
$\textrm{Cov}[\mathcal{B},\mathcal{D}]$), must be specified. Second-order
beliefs about $\mathcal{B}$ can be adjusted in the light of $\mathcal{D}$
using the Bayes linear update formulae:

$\displaystyle\textrm{E}_{\mathcal{D}}[\mathcal{B}]$ $\displaystyle=$
$\displaystyle\textrm{E}[\mathcal{B}]+\textrm{Cov}[\mathcal{B},\mathcal{D}]\textrm{Var}[\mathcal{D}]^{-1}(\mathcal{D}-\textrm{E}[\mathcal{D}]),$
(1) $\displaystyle\textrm{Var}_{\mathcal{D}}[\mathcal{B}]$ $\displaystyle=$
$\displaystyle\textrm{Var}[\mathcal{B}]-\textrm{Cov}[\mathcal{B},\mathcal{D}]\textrm{Var}[\mathcal{D}]^{-1}\textrm{Cov}[\mathcal{D},\mathcal{B}],$
(2) $\displaystyle\textrm{Cov}_{\mathcal{D}}[\mathcal{B}_{1},\mathcal{B}_{2}]$
$\displaystyle=$
$\displaystyle\textrm{Cov}[\mathcal{B}_{1},\mathcal{B}_{2}]-\textrm{Cov}[\mathcal{B}_{1},\mathcal{D}]\textrm{Var}[\mathcal{D}]^{-1}\textrm{Cov}[\mathcal{D},\mathcal{B}_{2}].$
(3)

Equations (1)-(3) are the backbone of a Bayes Linear Analysis, and permit
construction of Bayes linear statistical models, such as presented in the
later sections of this article. $\textrm{E}_{\mathcal{D}}[\mathcal{B}]$ and
$\textrm{Var}_{\mathcal{D}}[\mathcal{B}]$ are termed the adjusted expectation
and variance of $\mathcal{B}$ given $\mathcal{D}$.
$\textrm{Cov}_{\mathcal{D}}[\mathcal{B}_{1},\mathcal{B}_{2}]$ is termed the
adjusted covariance of $\mathcal{B}_{1}$ and $\mathcal{B}_{2}$ given
$\mathcal{D}$, where $\mathcal{B}_{1}$ and $\mathcal{B}_{2}$ are
subcollections of $\mathcal{B}$.

The Bayes linear approach differs both philosophically and practically to the
fully Bayesian approach to statistical inference (Goldstein, 1999). For
example, in practical problems, there are often many relevant sources of
uncertainty. A coherent fully Bayesian analysis requires specification of a
full joint prior probability distribution and likelihood to reflect beliefs
about the high-dimensional structure of these uncertainties (Garthwaite et
al., 2005). Such specification can be very difficult, hence approximations are
frequently made for mathematical convenience which causes the specification to
reflect some, but not all, aspects of a person’s beliefs. The theoretical
coherence of the full Bayesian analysis can then get lost due to the required
practical simplifications and assumptions. Furthermore, the resulting Bayesian
analysis is often too computationally intensive to carry out in reasonable
time. In contrast, by only requiring belief specification up to the second
order (Hartigan, 1969), uncertainty in model assumptions, along with any other
uncertainties, can be more easily incorporated into a Bayes linear analysis.
Since linear fitting is generally computationally simpler than full
conditioning, it can make for a more straightforward approach to the analysis
of complex problems. For a more detailed overview and thorough treatment of
Bayes linear methods, see Goldstein (1999) and Goldstein & Wooff (2007). For a
comparison of Bayes linear methods with the full Bayesian approach, see, for
example, Vernon et al. (2010).

We now present a general statistical model for
$\mathbf{y}(\mathbf{x})=(y_{1}(\mathbf{x}),\dots,y_{q}(\mathbf{x}))\in\mathbb{R}^{q}$
being a vector of quantities of interest assuming known input setting
$\mathbf{x}\in\mathbb{X}\subseteq\mathbb{R}^{p}$. The model is given by

$y_{i}(\mathbf{x})=f_{i}(\mathbf{x})+\epsilon_{i}(\mathbf{x}),\qquad\qquad
i=1,\dots,q.$ (4)

Here, $\mathbf{f}(\cdot)=(f_{1}(\cdot),\dots,f_{q}(\cdot))$ is a statistical
model of input $\mathbf{x}$, and
$\bm{\epsilon}(\cdot)=(\epsilon_{1}(\mathbf{x}),\dots,\epsilon_{q}(\mathbf{x}))$
is a residual process attempting to capture the discrepancy between
$\mathbf{f}$ and $\mathbf{y}$ under $\mathbf{x}$.

In order to fit the model, we assume a collection of observed training data
given by the $nq$-vector
$\mathcal{Y}=\mathcal{Y}(\mathcal{X})=(\mathbf{y}_{1}(\mathcal{X})^{T},\dots,\mathbf{y}_{q}(\mathcal{X})^{T})^{T}$
corresponding to each row of the $n\times p$ matrix
$\mathcal{X}=(\mathbf{x}^{(1)},\dots,\mathbf{x}^{(n)})^{T}$, with
$\mathbf{y}_{k}(\mathcal{X}),k=1,\dots,q$, being $n$-vectors of observations
corresponding to the $k$th quantity $y_{k}$.

We aim to adjust our second-order prior belief specification about
$\mathbf{y}(\mathbf{x})$ across $\mathcal{X}$ by $\mathcal{Y}$ using the Bayes
linear update equations to obtain posterior quantities:

${\rm E}_{\mathcal{Y}}[\mathbf{y}(\mathbf{x})]={\rm
E}[\mathbf{y}(\mathbf{x})]+{\rm Cov}[\mathbf{y}(\mathbf{x}),\mathcal{Y}]{\rm
Var}[\mathcal{Y}]^{-1}(\mathcal{Y}-{\rm E}[\mathcal{Y}]),$ (5) ${\rm
Cov}_{\mathcal{Y}}\left[\mathbf{y}(\mathbf{x}),\mathbf{y}(\mathbf{x}^{\prime})\right]={\rm
Cov}[\mathbf{y}(\mathbf{x}),\mathbf{y}(\mathbf{x}^{\prime})]-{\rm
Cov}[\mathbf{y}(\mathbf{x}),\mathcal{Y}]{\rm Var}[\mathcal{Y}]^{-1}{\rm
Cov}[\mathcal{Y},\mathbf{y}(\mathbf{x}^{\prime})]\,.$ (6)

The magnitude of the adjustment of our beliefs for $\mathbf{y}(\mathbf{x})$ in
light of $\mathcal{Y}$ is largely governed by the prior belief specification
${\rm E}[\mathbf{y}(\mathbf{x})]$ and ${\rm
Cov}[\mathbf{y}(\mathbf{x}),\mathbf{y}(\mathbf{x}^{\prime})]$ across
$\mathbb{X}$. For most statistical models, such prior belief specifications
can be expressed as smooth functions of $\mathbf{x}$ and
$\mathbf{x}^{\prime}$, this function being determined by the model structure
for $\mathbf{f}$ and $\bm{\epsilon}$ and the induced belief specification over
them. For example, we often consider a separable covariance structure of the
form

${\rm
Cov}[\mathbf{y}(\mathbf{x}),\mathbf{y}(\mathbf{x}^{\prime})]=c(\mathbf{x},\mathbf{x}^{\prime})\,\bm{\mathit{\Sigma}},$
(7)

where $c(\mathbf{x},\mathbf{x}^{\prime})$ is a stationary correlation function
of $\mathbf{x}$ and $\mathbf{x}^{\prime}$, and $\bm{\mathit{\Sigma}}$ is an
output covariance matrix. It is the structure of ${\rm
Cov}[\mathbf{y}(\mathbf{x}),\mathbf{y}(\mathbf{x}^{\prime})]$, and in
particular correlation function $c(\mathbf{x},\mathbf{x}^{\prime})$, that we
generalise to extend Bayes linear modelling to include uncertain inputs,
firstly in general terms here, and then in specific (but broad) contexts and
examples in the subsequent sections.

Consider the generalisation of the statistical model given by Equation (4) to
the following form:

$\mathbf{y}(\mathbf{X})=\mathbf{f}(\mathbf{X})+\bm{\epsilon}(\mathbf{X}),$ (8)

where now $\mathbf{X}$ is a random variable with second order belief
specification $\\{{\rm E}[\mathbf{X}],{\rm
Cov}[\mathbf{X},\mathbf{X}^{\prime}]\\}$ available for any
$\mathbf{X},\mathbf{X}^{\prime}$. Note that $y(\mathbf{X})$ is therefore the
random quantity in which interest lies.

In a Bayes linear setting, we wish to obtain

${\rm E}_{\mathcal{Y}}[\mathbf{y}(\mathbf{X})]={\rm
E}[\mathbf{y}(\mathbf{X})]+{\rm Cov}[\mathbf{y}(\mathbf{X}),\mathcal{Y}]{\rm
Var}[\mathcal{Y}]^{-1}(\mathcal{Y}-{\rm E}[\mathcal{Y}]),$ (9) ${\rm
Cov}_{\mathcal{Y}}\left[\mathbf{y}(\mathbf{X}),\mathbf{y}(\mathbf{X}^{\prime})\right]={\rm
Cov}[\mathbf{y}(\mathbf{X}),\mathbf{y}(\mathbf{X}^{\prime})]-{\rm
Cov}[\mathbf{y}(\mathbf{X}),\mathcal{Y}]{\rm Var}[\mathcal{Y}]^{-1}{\rm
Cov}[\mathcal{Y},\mathbf{y}(\mathbf{X}^{\prime})]\,,$ (10)

where training data
$\mathcal{Y}=\mathcal{Y}(\mathcal{X})=(\mathbf{y}_{1}(\mathcal{X})^{T},\dots,\mathbf{y}_{q}(\mathcal{X})^{T})^{T}$
are now observations corresponding to each row of
$\mathcal{X}=(\mathbf{X}^{(1)},\dots,\mathbf{X}^{(n)})^{T}$, each element of
which is specified up to second order. In other words, we have ${\rm
E}[\mathbf{X}^{(i)}]$ and ${\rm Cov}[\mathbf{X}^{(i)},\mathbf{X}^{(j)}]$ for
all $i,j=1,\dots,n$. Prior belief specification for $\mathcal{Y}$ is given by

$\displaystyle{\rm E}[\mathcal{Y}]$ $\displaystyle=$ $\displaystyle({\rm
E}[y(\mathbf{X}^{(1)})],\dots,{\rm E}[y(\mathbf{X}^{(n)})]),$ (11)
$\displaystyle{\rm Var}[\mathcal{Y}]$ $\displaystyle=$ $\displaystyle\\{{\rm
Cov}[y(\mathbf{X}^{(i)}),y(\mathbf{X}^{(j)})]\\}_{i,j=1,\dots,n}.$ (12)

Appropriate methodology to obtain ${\rm
E}_{\mathcal{Y}}[\mathbf{y}(\mathbf{X})]$ and ${\rm
Cov}_{\mathcal{Y}}\left[\mathbf{y}(\mathbf{X}),\mathbf{y}(\mathbf{X}^{\prime})\right]$
by Equations (9) and (10) is therefore achieved by specification of
appropriate prior belief structures ${\rm E}[\mathbf{y}(\mathbf{X})]$ and
${\rm Cov}[\mathbf{y}(\mathbf{X}),\mathbf{y}(\mathbf{X}^{\prime})]$ over any
possible specification of the set
$\\{{\rm E}[\mathbf{X}],{\rm E}[\mathbf{X}^{\prime}],{\rm
Cov}[\mathbf{X},\mathbf{X}^{\prime}]\\}$. Such belief structures often exist
as generalisations of those for known $\mathbf{x},\mathbf{x}^{\prime}$ alluded
to above. For example, the separable covariance structure given by Equation
(7) can be generalised to

${\rm
Cov}[\mathbf{y}(\mathbf{X}),\mathbf{y}(\mathbf{X}^{\prime})]=c(\mathbf{X},\mathbf{X}^{\prime})\,\bm{\mathit{\Sigma}},$
(13)

with

$c(\mathbf{X},\mathbf{X}^{\prime})=c({\rm E}[\mathbf{X}],{\rm
E}[\mathbf{X}^{\prime}],{\rm Cov}[\mathbf{X},\mathbf{X}^{\prime}]),$ (14)

implying that the correlation between $\mathbf{y}(\mathbf{X})$ and
$\mathbf{y}(\mathbf{X}^{\prime})$ for two uncertain inputs $\mathbf{X}$ and
$\mathbf{X}^{\prime}$ is a function of the second order belief specification
about and between the two input variables.

In the following sections, we proceed to discuss Bayes linear statistical
modelling with uncertain inputs in two specific, but broad, contexts; firstly,
that of regression analyses, and secondly, that of stochastic processes and
emulation.

##  3 Regression Analysis

Consider a regression model, presented here with scalar output $y$ for
simplicity of notation and exposition:

$y(\mathbf{X})=\mathbf{X}^{T}\bm{\beta}+\epsilon(\mathbf{X}),$ (15)

where $\bm{\beta}\in\mathbb{R}^{p}$ is a vector of exchangeable regression
parameters, and $\mathbf{X}\in\mathbb{R}^{p}$ is a vector random variable.

We consider a generic prior specification over $\bm{\beta}$ of the form ${\rm
E}[\bm{\beta}]=\bm{\Gamma}$ and ${\rm Var}[\bm{\beta}]=\bm{\mathit{\Delta}}$.
Under the usual scenario of known $\mathbf{x}$, we usually specify ${\rm
E}[\epsilon(\mathbf{x})]=0$ for all $\mathbf{x}$, thus we view it a reasonable
extension to specify ${\rm E}[\epsilon(\mathbf{X})]=0$ in the random variable
scenario. It is common for $\epsilon$ to be viewed as uncorrelated with
$\bm{\beta}$ _a priori_ , that is to specify ${\rm
Cov}[\bm{\beta},\epsilon(\mathbf{x})]=0$, so again we view it as reasonable to
extend this in the random variable context to ${\rm
Cov}[\bm{\beta},\epsilon(\mathbf{X})]=0$.

Following Equations (9) and (10), the required prior specifications are ${\rm
E}[y(\mathbf{X})]$ and ${\rm Cov}[y(\mathbf{X}),y(\mathbf{X}^{\prime})]$ given
any possible specification of the set $\\{{\rm E}[\mathbf{X}],{\rm
E}[\mathbf{X}^{\prime}],{\rm Cov}[\mathbf{X},\mathbf{X}^{\prime}]\\}$ for
random variables $\mathbf{X},\mathbf{X}^{\prime}$. In the remainder of this
section, we proceed to state results related to appropriate expression of
prior beliefs of this form, both at a general level and for two specific
examples of possible residual structure. Extended derivation of these results
is presented in the supplementary material.

We have that:

$\displaystyle{\rm E}[y(\mathbf{X})]$ $\displaystyle=$ $\displaystyle{\rm
E}[\mathbf{X}^{T}]\,\bm{\Gamma},$ (16) $\displaystyle{\rm
Cov}[y(\mathbf{X}),y(\mathbf{X}^{\prime})]$ $\displaystyle=$
$\displaystyle{\rm Cov}[\mathbf{X}^{T}\bm{\beta},\mathbf{X}^{\prime
T}\bm{\beta}]+{\rm Cov}[\epsilon(\mathbf{X}),\epsilon(\mathbf{X}^{\prime})],$
(17)

noting that ${\rm
Cov}[\mathbf{X}^{T}\bm{\beta},\epsilon(\mathbf{X}^{\prime})]=0$. Further,
using the law of total covariance,

$\displaystyle{\rm Cov}[\mathbf{X}^{T}\bm{\beta},\mathbf{X}^{\prime
T}\bm{\beta}]$ $\displaystyle=$ $\displaystyle{\rm
E}[\mathbf{X}^{T}\bm{\mathit{\Delta}}\mathbf{X}^{\prime}]+\bm{\Gamma}^{T}{\rm
Cov}[\mathbf{X},\mathbf{X}^{\prime}]\bm{\Gamma}.$ (18)

The second term of Equation (79) is a covariance specification over the
residual process $\epsilon(\mathbf{X})$. In Sections A.1 and A.2, we present
two example model belief structures in the context of statistical modelling
with known inputs, demonstrating reasonable generalisation to the random
variable input scenario using the ideas presented in Section 2.

###  Linear Regression with Uncorrelated Random Error

In the common known input case, the standard linear regression model with
uncorrelated and homoscedastic random error can be obtained by specifying a
prior covariance structure over $\epsilon(\mathbf{x})$ of

${\rm
Cov}[\epsilon(\mathbf{x}),\epsilon(\mathbf{x}^{\prime})]=\mathbb{I}_{\mathbf{x}=\mathbf{x}^{\prime}}\sigma^{2},$
(19)

for all $\mathbf{x},\mathbf{x}^{\prime}$, where $\mathbb{I}$ is the indicator
function taking value $1$ if the statement is true and $0$ otherwise, and
$\sigma^{2}$ is a common scalar variance parameter.

This specification can be simply extended to the random variable case as
follows:

${\rm
Cov}[\epsilon(\mathbf{X}),\epsilon(\mathbf{X}^{\prime})]=\mathbb{I}_{\mathbf{X}=\mathbf{X}^{\prime}}\sigma^{2},$
(20)

where the covariance between two random variable inputs is zero unless they
are known to be the same. Note that two random variables being the same is not
equivalent to two different random variables having the same second-order
belief specification.

###  Linear Regression with Correlated Error - An Example

In this section, we develop an uncertain input Bayes linear approach for
exchangeable regressions with correlated errors, building on the example of
extracting aluminium by electrolysis over time presented in Goldstein & Wooff
(1998). Using the covariance structure presented in that paper for known
inputs, we present a generalisation to the random variable input scenario
using the ideas presented in Section 2.

The model is of the form

$y(t)=\beta_{0}+\beta_{1}\,t+\epsilon(t),$ (21)

and is exchangeable over $\bm{\beta}$. Note that this regression model is
consistent with the general form presented in Equation (73), with $p=2$
variables being an intercept and single controllable parameter $t$. In
Goldstein & Wooff (1998), it is assumed that $t=1,\dots,13$, with a structured
error model as follows:

$\displaystyle\epsilon(t)$ $\displaystyle=$ $\displaystyle A(t)+J(t)+H(t),$
(22) $\displaystyle J(t)$ $\displaystyle=$ $\displaystyle J(t-1)+Q(t),$
$\displaystyle H(t)$ $\displaystyle=$ $\displaystyle\psi\,H(t-1)+R(t),t\geq
2,$

with $\psi\in(0,1)$, and where, for generality, we denote ${\rm
Var}[A(t)]=\sigma^{2}_{A}$, ${\rm Var}[Q(t)]={\rm Var}[J(1)]=\sigma^{2}_{Q}$,
${\rm Var}[R(t)]=\sigma^{2}_{R}$, and ${\rm Var}[H(1)]=\sigma^{2}_{1}$. These
terms express discrepancies from the linear trend as the sum of a pure
measurement error $A(t)$, a stochastic development of the discrepancy as a
random walk with drift $J(t)$, and an autoregressive term expressing the
measurement of the suspended particles in the chemical analysis $H(t)$. For
further details, see Goldstein & Wooff (1998). Whilst the original model
assumes discrete timesteps $t$, we will generalise the regression structure to
continuous $t$-values, subject to $t\geq 3$ for consistency issues, before
addressing the uncertain input $T$ scenario.

We assume that the three terms in Equation (22) are uncorrelated, leading to

$\displaystyle{\rm Cov}[\epsilon(t),\epsilon(t)]$ $\displaystyle=$
$\displaystyle{\rm Cov}[A(t),A(t^{\prime})]+{\rm Cov}[J(t),J(t^{\prime})]+{\rm
Cov}[H(t),H(t^{\prime})].$ (23)

The three terms on the right hand side of Equation 88 have the following
structure, where we define $d=|t-t^{\prime}|$, and $t_{m}=\min(t,t^{\prime})$,
with greater exposition presented in the supplementary material:

$\displaystyle{\rm Cov}[A(t),A(t^{\prime})]$ $\displaystyle=$
$\displaystyle\mathbb{I}_{t=t^{\prime}}\sigma_{A}^{2},$ (24)
$\displaystyle{\rm Cov}[J(t),J(t^{\prime})]$ $\displaystyle=$ $\displaystyle
t_{m}\,\sigma^{2}_{Q},$ (25) $\displaystyle{\rm Cov}[H(t),H(t^{\prime})]$
$\displaystyle=$
$\displaystyle\psi^{d}\,\left(\psi^{2(t_{m}-1)}\,\sigma^{2}_{1}+\frac{1-\psi^{2(t_{m}-2)}}{1-\psi^{2}}\,\sigma_{R}^{2}\right).$
(26)

Hence we have a prior residual covariance belief structure across all $t$ of

$\displaystyle{\rm Cov}[\epsilon(t),\epsilon(t)]$ $\displaystyle=$
$\displaystyle\mathbb{I}_{t=t^{\prime}}\sigma_{A}^{2}+t_{m}\,\sigma^{2}_{Q}+\psi^{d}\,\left(\psi^{2(t_{m}-1)}\,\sigma^{2}_{1}+\frac{1-\psi^{2(t_{m}-2)}}{1-\psi^{2}}\,\sigma_{R}^{2}\right).$
(27)

We now proceed to consider suitable extensions to the uncertain input
scenario. In this case, we need an expectation and covariance structure for
any possible specification of the set $\\{{\rm E}[T],{\rm E}[T^{\prime}],{\rm
Cov}[T,T^{\prime}]\\}$. We consider that the statements of expectation in the
error structure simply extend to being ${\rm E}[A(T)]=0$, ${\rm E}[J(T)]=0$
and ${\rm E}[H(T)]=0$; logical derivation of these statements can be shown
using the law of total expectation. Regarding covariance structure, Equation
(89) simply extends to the following:

${\rm Cov}[A(T),A(T^{\prime})]=\mathbb{I}_{T=T^{\prime}}\sigma_{A}^{2},$ (28)

this being similar to the uncorrelated error term presented in Section A.1.

Equation (90) extends as follows, where the random variables
$T_{m}=\min(T,T^{\prime})$ and $D=|T-T^{\prime}|$ correspond to known
quantities $t_{m}$ and $d$ above:

$\displaystyle{\rm Cov}[J(T),J(T^{\prime})]$ $\displaystyle=$
$\displaystyle\sigma^{2}_{Q}\,{\rm E}[T_{m}].$ (29)

The tricky specification here is eliciting ${\rm E}[T_{m}]$ from ${\rm
E}[T],{\rm E}[T^{\prime}]$ and ${\rm Cov}[T,T^{\prime}]$, since a function
which takes the minimum of two random quantities is non-linear. If we know
which of our random variables is bigger (without loss of generality, that
$T<T^{\prime}$, say), then ${\rm E}[T_{m}]={\rm E}[T]$. This is not an
unrealistic situation; for example, we may have measurements at two times
about which we are uncertain, but know that one of them certainly happened
after the other.

More generally, we propose a reasonable covariance structure as follows. Let
$\mathbb{I}=\mathbb{I}_{T<T^{\prime}}$ and $p={\rm E}[\mathbb{I}]$. We then
condition on this event using the law of total expectation to give:

$\displaystyle{\rm E}[T_{m}]$ $\displaystyle=$ $\displaystyle p\,{\rm
E}_{\mathbb{I}=1}[T_{m}]+(1-p)\,{\rm E}_{\mathbb{I}=0}[T_{m}],$ (30)

where

$\displaystyle{\rm E}_{\mathbb{I}=1}[T_{m}]$ $\displaystyle=$
$\displaystyle{\rm E}[T]+\frac{{\rm Cov}[T,\mathbb{I}]}{p},$ (31)
$\displaystyle{\rm E}_{\mathbb{I}=0}[T_{m}]$ $\displaystyle=$
$\displaystyle{\rm E}[T^{\prime}]-\frac{{\rm
Cov}[T^{\prime},\mathbb{I}]}{1-p},$ (32)

and

$\displaystyle{\rm Cov}[T,\mathbb{I}]$ $\displaystyle=$ $\displaystyle{\rm
Corr}\left[T,\mathbb{I}\right]\,\sqrt{{\rm Var}[T]\,p(1-p)},$ (33)
$\displaystyle{\rm Cov}[T^{\prime},\mathbb{I}]$ $\displaystyle=$
$\displaystyle{\rm Corr}\left[T^{\prime},\mathbb{I}\right]\,\sqrt{{\rm
Var}[T^{\prime}]\,p(1-p)},$ (34)

so that:

$\displaystyle{\rm E}[T_{m}]$ $\displaystyle=$ $\displaystyle p\,{\rm
E}[T]+(1-p)\,{\rm E}[T^{\prime}]$ (35)
$\displaystyle\,\,\,\,\,\,+\sqrt{p(1-p)}\left({\rm
Corr}\left[T,\mathbb{I}\right]\,\sqrt{{\rm Var}[T]}-{\rm
Corr}\left[T^{\prime},\mathbb{I}\right]\,\sqrt{{\rm Var}[T^{\prime}]}\right).$

We may be able to specify values for all quantities required in Equation
(103), in which case we have our expression for ${\rm E}[T_{m}]$. If
specification of ${\rm Corr}\left[T,\mathbb{I}\right]$ and ${\rm
Corr}\left[T^{\prime},\mathbb{I}\right]$ is proving difficult, then let us
first note that it is logical to assume that $-1<{\rm
Corr}\left[T,\mathbb{I}\right]<0$ (since $\mathbb{I}=1$ necessarily
corresponds, on average in some sense, to smaller $T$, as in this case
$T<T^{\prime}$), and $0<{\rm Corr}\left[T^{\prime},\mathbb{I}\right]<1$, so
that

${\rm E}[T_{m}]>p\,{\rm E}[T]+(1-p)\,{\rm
E}[T^{\prime}]-\sqrt{p(1-p)}\left(\sqrt{{\rm Var}[T]}+\sqrt{{\rm
Var}[T^{\prime}]}\right).$ (36)

Again, if we feel able to specify a value for $p$, then we can use it in
Equation (104). Alternatively, in order to be conservative with respect to
variance resolution, we can use elementary calculus to minimise the expression
on the right had side of Equation (104) over $p\in[0,1]$ to get that

${\rm E}[T_{m}]\geq\frac{1}{2}\left({\rm E}[T]+{\rm
E}[T^{\prime}]-\sqrt{\left({\rm E}[T^{\prime}]-{\rm
E}[T]\right)^{2}+\left(\sqrt{{\rm Var}[T]}+\sqrt{{\rm
Var}[T^{\prime}]}\right)^{2}}\right).$ (37)

We thus propose, for $T\neq T^{\prime}$, the covariance structure:

${\rm Cov}[J(T),J(T^{\prime})]=\frac{\sigma^{2}_{Q}}{2}\left({\rm E}[T]+{\rm
E}[T^{\prime}]-\sqrt{\left({\rm E}[T^{\prime}]-{\rm
E}[T]\right)^{2}+\left(\sqrt{{\rm Var}[T]}+\sqrt{{\rm
Var}[T^{\prime}]}\right)^{2}}\right).$ (38)

From Equation (90) with $t=t^{\prime}$, it follows that ${\rm
Var}[J(T)]=\sigma^{2}_{Q}\,{\rm E}[T],$ since $T_{m}=T$. We thus make use of
an indicator function in the expression for ${\rm Cov}[J(T),J(T^{\prime})]$ to
ensure it is appropriate when it is known that $T=T^{\prime}$:

$\displaystyle{\rm Cov}[J(T),J(T^{\prime})]$ $\displaystyle=$
$\displaystyle\frac{\sigma^{2}_{Q}}{2}\left({\rm E}[T]+{\rm
E}[T^{\prime}]-\sqrt{\left({\rm E}[T^{\prime}]-{\rm
E}[T]\right)^{2}+\left(1-\mathbb{I}_{T=T^{\prime}}\right)\left(\sqrt{{\rm
Var}[T]}+\sqrt{{\rm Var}[T^{\prime}]}\right)^{2}}\right)$

Note that:

  1. 1.

When ${\rm Var}[T]={\rm Var}[T^{\prime}]=0$, the expression on the right hand
side of Inequality (112) yields

$\frac{1}{2}\left(t+t^{\prime}-\sqrt{(t-t^{\prime})^{2}}\right)=\frac{1}{2}\left(t+t^{\prime}-|t-t^{\prime}|\right)=\min(t,t^{\prime}),$

thus resolving to the known input scenario.

  2. 2.

If it is known that $T<T^{\prime}$, then, as mentioned earlier, it would be
logical to have that ${\rm E}[T_{m}]={\rm E}[T]$. Since

${\rm E}[T]=\frac{1}{2}\left({\rm E}[T]+{\rm E}[T^{\prime}]-\sqrt{\left({\rm
E}[T^{\prime}]-{\rm E}[T]\right)^{2}}\right),$ (40)

it is clear that ${\rm E}[T_{m}]$ satisfies Inequality (112) in this case.

We extend Equation (92) by noticing that

$\displaystyle{\rm Cov}[H(T),H(T^{\prime})]$ $\displaystyle\geq$
$\displaystyle\sigma_{1}^{2}\exp(\log\psi\,{\rm
E}[T_{M}+T_{m}-2])+\sigma_{R}^{2}\exp(\log\psi\,{\rm E}[T_{M}-T_{m}]),$ (41)

where $t_{M}=\max(t,t^{\prime})$, and correspondingly
$T_{M}=\max(T,T^{\prime})$.

We note that ${\rm E}[T_{M}+T_{m}]={\rm E}[T]+{\rm E}[T^{\prime}]$, and thus
seek an upper bound for ${\rm E}[T_{M}-T_{m}]$ in order to find a lower bound
for Expression (117) (since $\log\psi<0$ if $\psi\in(0,1)$). Similar to
Equation (98), we use the law of total expectation on ${\rm E}[T_{M}]$ to get
that

$\displaystyle{\rm E}[T_{M}]$ $\displaystyle=$ $\displaystyle p\,{\rm
E}_{\mathbb{I}=1}[T_{M}]+(1-p)\,{\rm E}_{\mathbb{I}=0}[T_{M}],$ (42)

where:

$\displaystyle{\rm E}_{\mathbb{I}=1}[T_{M}]$ $\displaystyle=$
$\displaystyle{\rm E}[T^{\prime}]+\frac{{\rm Cov}[T^{\prime},\mathbb{I}]}{p},$
(43) $\displaystyle{\rm E}_{\mathbb{I}=0}[T_{M}]$ $\displaystyle=$
$\displaystyle{\rm E}[T]-\frac{{\rm Cov}[T,\mathbb{I}]}{1-p}.$ (44)

Combining Equations (98) and (118), we get that

$\displaystyle{\rm E}[T_{M}-T_{m}]$ $\displaystyle=$
$\displaystyle(1-2p)\left({\rm E}[T]-{\rm E}[T^{\prime}]\right)$ (45)
$\displaystyle\,\,\,\,\,\,+2\sqrt{p(1-p)}\left({\rm
Corr}\left[T^{\prime},\mathbb{I}\right]\sqrt{{\rm Var}[T^{\prime}]}-{\rm
Corr}\left[T,\mathbb{I}\right]\sqrt{{\rm Var}[T]}\right).$

Again, if we don’t feel able to specify ${\rm
Corr}\left[T^{\prime},\mathbb{I}\right]$ or ${\rm
Corr}\left[T,\mathbb{I}\right]$, then using the assumptions that $-1<{\rm
Corr}\left[T,\mathbb{I}\right]<0$ and $0<{\rm
Corr}\left[T^{\prime},\mathbb{I}\right]<1$ we have that

${\rm E}[T_{M}-T_{m}]<(1-2p)\left({\rm E}[T]-{\rm
E}[T^{\prime}]\right)+2\sqrt{p(1-p)}\left(\sqrt{{\rm
Var}[T^{\prime}]}+\sqrt{{\rm Var}[T]}\right).$ (46)

If we have specified a value for $p$, again we can use this in Expression
(122), otherwise we can maximise over $p\in[0,1]$ to obtain

$\displaystyle{\rm E}[T_{M}-T_{m}]$ $\displaystyle\leq$
$\displaystyle\sqrt{\left({\rm E}[T]-{\rm
E}[T^{\prime}]\right)^{2}+\left(\sqrt{{\rm Var}[T]}+\sqrt{{\rm
Var}[T^{\prime}]}\right)^{2}}.$ (47)

For $T\neq T^{\prime}$, we therefore propose

$\displaystyle{\rm Cov}[H(T),H(T^{\prime})]$ (48) $\displaystyle=$
$\displaystyle\sigma_{1}^{2}\exp(\log\psi\,({\rm E}[T]+{\rm
E}[T^{\prime}]-2))$
$\displaystyle\,\,\,\,\,\,+\,\sigma_{R}^{2}\exp\left(\log\psi\,\sqrt{\left({\rm
E}[T]-{\rm E}[T^{\prime}]\right)^{2}+\left(\sqrt{{\rm Var}[T]}+\sqrt{{\rm
Var}[T^{\prime}]}\right)^{2}}\right).$

When $T=T^{\prime}$, the above expression needs to be commensurate with a
reasonable expression for ${\rm Var}[H(T)]$. Unlike ${\rm Var}[J(T)]$, we find
that ${\rm Var}[H(T)]$ itself needs approximating. However, for a variance
term we need to find an upper bound in order to be conservative with regard to
prior variance specification (i.e. overestimate rather than underestimate). We
have that

$\displaystyle{\rm Var}[H(T)]$ $\displaystyle\leq$
$\displaystyle\mathbb{I}_{\sigma_{1}^{2}\leq
K}K+(1-\mathbb{I}_{\sigma_{1}^{2}\leq
K})(\psi^{4}\sigma_{1}^{2}+\sigma_{R}^{2}),$ (49)

where $K=\frac{\sigma_{R}^{2}}{1-\psi^{2}}$ (see the supplementary material
for justification).

We can now set

$\displaystyle{\rm Cov}[H(T),H(T^{\prime})]$ (50) $\displaystyle=$
$\displaystyle\mathbb{I}_{T=T^{\prime}}L+(1-\mathbb{I}_{T=T^{\prime}})\bigg{(}\sigma_{1}^{2}\exp(\log\psi\,({\rm
E}[T]+{\rm E}[T^{\prime}]-2))$
$\displaystyle\,\,\,\,\,\,+\,\sigma_{R}^{2}\exp\left(\log\psi\,\sqrt{\left({\rm
E}[T]-{\rm E}[T^{\prime}]\right)^{2}+\left(\sqrt{{\rm Var}[T]}+\sqrt{{\rm
Var}[T^{\prime}]}\right)^{2}}\right)\bigg{)},$

where

$L=\mathbb{I}_{\sigma_{1}^{2}\leq K}K+(1-\mathbb{I}_{\sigma_{1}^{2}\leq
K})(\psi^{4}\sigma_{1}^{2}+\sigma_{R}^{2}).$

Note that Expression (50) underestimates ${\rm Cov}[H(t),H(t^{\prime})]$, as
given by Equation (92), when ${\rm Var}[T]={\rm Var}[T^{\prime}]=0$. We do not
view this as a problem, as we consider interest to lie in the case of
statistical modelling with random variable input. If it is assumed that all
$t,t^{\prime}$ are known, we would use Equation (92) throughout.

Finally, to conclude our prior specification in the case of uncertain inputs
for the model presented in this section, note that, with
$\mathbf{T}=(1,T)^{T}$, we have

$\displaystyle{\rm E}[\mathbf{T}^{T}\bm{\mathit{\Delta}}\mathbf{T}^{\prime}]$
$\displaystyle=$ $\displaystyle\delta_{00}+\delta_{01}({\rm
E}[T^{\prime}]+{\rm E}[T])+\delta_{11}{\rm E}[T]{\rm E}[T^{\prime}]+{\rm
Cov}[T,T^{\prime}],$ (51) $\displaystyle\bm{\Gamma}^{T}{\rm
Cov}[\mathbf{T},\mathbf{T}^{\prime}]\bm{\Gamma}$ $\displaystyle=$
$\displaystyle\gamma_{1}^{2}{\rm Cov}[T,T^{\prime}],$ (52)

so that full prior specification for $y(T)$ can therefore be given by

${\rm E}[y(T)]=\gamma_{0}+{\rm E}[T]\gamma_{1},$ (53)

and

$\displaystyle{\rm Cov}[y(T),y(T^{\prime})]$ (54) $\displaystyle=$
$\displaystyle\delta_{00}+\delta_{01}({\rm E}[T^{\prime}]+{\rm
E}[T])+\delta_{11}{\rm E}[T]{\rm E}[T^{\prime}]+{\rm Cov}[T,T^{\prime}]$
$\displaystyle\,\,\,\,\,\,+\,\,\gamma_{1}^{2}{\rm Cov}[T,T^{\prime}]$
$\displaystyle\,\,\,\,\,\,+\,\,\mathbb{I}_{T=T^{\prime}}\sigma_{A}^{2}$
$\displaystyle\,\,\,\,\,\,+\,\,\frac{\sigma^{2}_{Q}}{2}\left({\rm E}[T]+{\rm
E}[T^{\prime}]-\sqrt{\left({\rm E}[T^{\prime}]-{\rm
E}[T]\right)^{2}+\left(1-\mathbb{I}_{T=T^{\prime}}\right)\left(\sqrt{{\rm
Var}[T]}+\sqrt{{\rm Var}[T^{\prime}]}\right)^{2}}\right)$
$\displaystyle\,\,\,\,\,\,+\,\,\mathbb{I}_{T=T^{\prime}}L+(1-\mathbb{I}_{T=T^{\prime}})\bigg{(}\sigma_{1}^{2}\exp(\log\psi\,({\rm
E}[T]+{\rm E}[T^{\prime}]-2))$
$\displaystyle\,\,\,\,\,\,\qquad+\sigma_{R}^{2}\exp\left(\log\psi\,\sqrt{\left({\rm
E}[T]-{\rm E}[T^{\prime}]\right)^{2}+\left(\sqrt{{\rm Var}[T]}+\sqrt{{\rm
Var}[T^{\prime}]}\right)^{2}}\right)\bigg{)}.$

Given these prior specifications, we can use the Bayes linear update Equations
(9) and (10) to adjust belief specifications for any second-order
specification $\\{{\rm E}[T],{\rm Var}[T]\\}$ for $T$, given knowledge of
system behaviour $y(\mathcal{T})$ for training set
$\mathcal{T}=(T^{(1)},\dots,T^{(n)})^{T}$. Here, beliefs across $\mathcal{T}$
are specified up to second order.

##  4 Bayes Linear Stochastic Processes

We here consider modelling the relationship between input
$\mathbf{x}\in\mathbb{R}^{p}$ and output
$\mathbf{y}(\mathbf{x})\in\mathbb{R}^{q}$ using a second-order stochastic
process, which provides an intuitive setting for the discussion of many of our
novel extensions in Uncertain Input modelling required for Bayes linear
emulation in Section 5. More precisely, we assume

$\mathbf{y}(\mathbf{x})=\mathbf{u}(\mathbf{x}),$ (55)

with a prior assumption that ${\rm E}[\mathbf{u}(\mathbf{x})]=0$, and assuming
that the between-outputs covariance and spatial correlation across the input
space is separable, thus given by

${\rm
Cov}[\mathbf{u}(\mathbf{x}),\mathbf{u}(\mathbf{x}^{\prime})]=c(\mathbf{x},\mathbf{x}^{\prime})\,\bm{\mathit{\Sigma}},$
(56)

for two inputs $\mathbf{x}$ and $\mathbf{x}^{\prime}$. Here,
$\bm{\mathit{\Sigma}}$ is a $q\times q$ output covariance matrix and
$c(\mathbf{x},\mathbf{x}^{\prime})$ is a stationary correlation function of
$\mathbf{x}$ and $\mathbf{x}^{\prime}$ (Koehler & Owen, 1996; Kennedy &
O’Hagan, 2001); for example, the Gaussian correlation function for modelling a
deterministic relationship (for example, that of a computer experiment) is
given by

$c(\mathbf{x},\mathbf{x}^{\prime})=\exp\left\\{-\sum_{r=1}^{p}\left(\frac{x_{(r)}-x^{\prime}_{(r)}}{\theta_{r}}\right)^{2}\right\\},$
(57)

which depends on the specification of the correlation length parameters
$\theta_{r},r=1,...,p$.

The second-order stochastic process in Equation (55) can be viewed similar to
a Gaussian process in the fully Bayesian paradigm with corresponding second-
order specification. Gaussian process modelling with uncertain inputs has been
considered (Dallaire et al., 2009; McHutchon & Rasmussen, 2011; Ye et al.,
2022; Wang et al., 2022), however, the uncertain input location noise is often
assumed to be homoscedastic and unknown, whereas we assume specification of
$\\{{\rm E}[\mathbf{X}],{\rm E}[\mathbf{X}^{\prime}],{\rm
Cov}[\mathbf{X},\mathbf{X}^{\prime}]\\}$, but for the modelling structure to
hold for any such specification.

In a Bayes linear setting with known inputs, the assumed covariance structure,
for example as presented as separable in Equation (56) and for a specific
example form in Equation (57), makes it straightforward to adjust our second-
order prior belief specification about $\mathbf{y}(\mathbf{x})$ across
$\mathcal{X}$ by $\mathbf{Y}$ to obtain posterior quantities, using the Bayes
linear update Equations (5) and (6).

Extension of Bayes linear stochastic process modelling to the random variable
input case with second-order specification requires extension of the prior
belief specifications for $\mathbf{u}$, such as given by Equations (56) and
(57). We propose ${\rm E}[\mathbf{u}(\mathbf{X})]=0$ and

${\rm
Cov}[\mathbf{u}(\mathbf{X}),\mathbf{u}(\mathbf{X}^{\prime})]=c(\mathbf{X},\mathbf{X}^{\prime})\,\bm{\mathit{\Sigma}}$
(58)

where

$c(\mathbf{X},\mathbf{X}^{\prime})=c({\rm E}[\mathbf{X}],{\rm
E}[\mathbf{X}^{\prime}],{\rm Cov}[\mathbf{X},\mathbf{X}^{\prime}])$ (59)

which implies that the correlation between $\mathbf{u}(\mathbf{X})$ and
$\mathbf{u}(\mathbf{X}^{\prime})$ for two uncertain inputs $\mathbf{X}$ and
$\mathbf{X}^{\prime}$ is a function of the second order belief specification
about and between the two input variables. As an example, we propose the
following extension to the Gaussian correlation function, given by Equation
(57):

$\begin{split}c(\mathbf{X},\mathbf{X}^{\prime})&=\exp\left\\{-{\rm
E}[(\mathbf{X}-\mathbf{X}^{\prime})^{T}\bm{\mathit{\Theta}}^{-2}(\mathbf{X}-\mathbf{X}^{\prime})]\right\\}\\\
&=\exp\left\\{-\sum_{r=1}^{p}\left(\frac{\textrm{E}[X_{(r)}-X^{\prime}_{(r)}]^{2}+\textrm{Var}[X_{(r)}-X^{\prime}_{(r)}]}{\theta^{2}_{r}}\right)\right\\}\,,\end{split}$
(60)

with positive-definite diagonal matrix $\bm{\mathit{\Theta}}^{-2}$ having
entries $(\Theta^{-2})_{rr}=1/\theta_{r}^{2}$, and the second line obtained
from standard results on the expected value of a quadratic form (Harville,
2018, pp. 200-201). This choice satisfies the desirable property that it
reduces to a standard form of correlation function if
$(\mathbf{X},\mathbf{X}^{\prime})=(\mathbf{x},\mathbf{x}^{\prime})$ are known.
We also derive two further important results for this new correlation function
in the form of two lemmas, proofs of which can be found in the supplementary
material.

###### Lemma 4.0.1

For random variables $\mathbf{X},\mathbf{X}^{\prime}$ with finite first and
second moments, the kernel function $c(\mathbf{X},\mathbf{X}^{\prime})$ from
(60) is positive-definite.

###### Lemma 4.0.2

Covariance function (59), with $c(\mathbf{X},\mathbf{X}^{\prime})$ given by
(60), is a lower bound under the Loewner (partial) ordering on the covariance
obtained by assuming the conditional covariance

${\rm
Cov}[\mathbf{u}(\mathbf{X}),\mathbf{u}(\mathbf{X}^{\prime})\,|\,\mathbf{X}=\mathbf{x},\mathbf{X}^{\prime}=\mathbf{x}^{\prime}]=\exp\left\\{-\sum_{r=1}^{p}\left(\frac{x_{(r)}-x^{\prime}_{(r)}}{\theta_{r}}\right)^{2}\right\\}\,\bm{\mathit{\Sigma}}\,.$
(61)

Moreover, the $kk$th element of (60) is a lower bound on the expected value,
with respect to $\mathbf{X},\mathbf{X}^{\prime}$, of the $kk$th element of
(61) ($k=1,\ldots q$).

Whilst the proposed correlation function of Equation (60) can be viewed simply
as a modelling assumption, Lemma 4.0.2 shows that it can also be derived as an
approximation to the conditional covariance between $\mathbf{u}(\mathbf{X})$
and $\mathbf{u}(\mathbf{X}^{\prime})$ assuming the standard Gaussian
correlation function that one might use for known
$\mathbf{x},\mathbf{x}^{\prime}$. In terms of a Bayes linear adjustment, this
quantity reflects the amount of resolved uncertainty given the training data,
hence an underestimation (lower bound) of this quantity is preferable to an
overestimation. Similar derivations for other correlation function forms
commonly presented in the literature (for example, given by Paulo, 2005), is
an avenue for future research.

##  5 Uncertain Input Bayes Linear Emulation

We extend uncertain input modelling methodology to Bayes linear emulation,
which can in some sense be viewed as combining the applications of Sections 3
and 4 together. Emulators are typically utilised as statistical approximations
of mathematical models, represented as computer code, or simulators, of
scientific processes. Such mathematical models encapsulate the key features of
the system and facilitate prediction and decision making. Emulation is often
required as a result of the computational expense of a typical simulator
leading to substantial output uncertainty across the input space due to the
small number of input combinations for which it is feasible to run the
simulator. The difficulties encountered as a result of the overall modelling
uncertainties are exacerbated in the case of Uncertain Inputs, as discussed in
this section. We begin with a review of Bayes linear emulation.

###  Bayes Linear Emulation

We assume system behaviour of interest $\mathbf{y}$ is modelled by a generic
deterministic simulator $\mathbf{f}$. The simulator has input vector
$\mathbf{x}=(x_{(1)},\dots,x_{(p)})\in\mathbb{X}\subseteq\mathbb{R}^{p}$, and
outputs vector
$\mathbf{f}(\mathbf{x})=(f_{1}(\mathbf{x}),\dots,f_{q}(\mathbf{x}))\in\mathbf{f}(\mathbb{X})\subseteq\mathbb{R}^{q}$.
We represent our beliefs about the behaviour of $\mathbf{f}(\mathbf{x})$ in
the following form (Goldstein & Rougier, 2004):

$\mathbf{f}(\mathbf{x})=\mathbf{g}(\mathbf{x})^{T}\bm{B}+\mathbf{u}(\mathbf{x})\,,$
(62)

where $\mathbf{g}(\mathbf{x})$ is an $m$-vector of known basis regression
functions, $\bm{B}$ is an $m\times q$ matrix of unknown regression
coefficients, and $\mathbf{u}(\mathbf{x})$ is a $q$-dimensional second-order
weakly-stationary stochastic process.

Let $\bm{\beta}=\textrm{vec}{(\bm{B})}$ be an $mq$-vector resulting from
stacking the columns of $\bm{B}$, with a generic prior specification ${\rm
E}[\bm{\beta}]=\bm{\Gamma}$ and ${\rm Var}[\bm{\beta}]=\bm{\mathit{\Delta}}$.
We also make the common assumptions that ${\rm
E}[\mathbf{u}(\mathbf{x})]=\mathbf{0}$, ${\rm
Cov}[\bm{\beta},\mathbf{u}(\mathbf{x})]=\bm{\mathit{0}}$, and covariance
between $\mathbf{u}(\mathbf{x})$ and $\mathbf{u}(\mathbf{x}^{\prime})$ is of
the separable form given by Equation (56).

Consistent with previous notation, suppose
$\mathbf{F}=\mathbf{f}(\mathcal{X})=(\mathbf{f}_{1}(\mathcal{X})^{T},\dots,\mathbf{f}_{q}(\mathcal{X})^{T})^{T}$
is an $nq$-vector with $\mathbf{f}_{k}(\mathcal{X})$ being $n$-vectors of
simulator output $k=1,\dots,q$ run at each row of the $n\times p$ design
matrix $\mathcal{X}=(\mathbf{x}^{(1)},\dots,\mathbf{x}^{(n)})^{T}$. We can
adjust our second-order prior belief specification about
$\mathbf{f}(\mathbf{x})$ across $\mathcal{X}$ by $\mathbf{F}$ using the Bayes
linear update Equations (63) and (64), now with $\mathbf{y}$ replaced with
$\mathbf{f}$, to obtain posterior quantities:

${\rm E}_{\mathbf{F}}[\mathbf{f}(\mathbf{x})]={\rm
E}[\mathbf{f}(\mathbf{x})]+{\rm Cov}[\mathbf{f}(\mathbf{x}),\mathbf{F}]{\rm
Var}[\mathbf{F}]^{-1}(\mathbf{F}-{\rm E}[\mathbf{F}]),$ (63) ${\rm
Var}_{\mathbf{F}}[\mathbf{f}(\mathbf{x})]={\rm
Var}[\mathbf{f}(\mathbf{x})]-{\rm Cov}[\mathbf{f}(\mathbf{x}),\mathbf{F}]{\rm
Var}[\mathbf{F}]^{-1}{\rm Cov}[\mathbf{F},\mathbf{f}(\mathbf{x})]\,.$ (64)

###  Bayes Linear Emulation with Uncertain Inputs

Uncertain Input Bayes Linear Emulation (UIBLE) is a computationally efficient
extension to Bayes linear emulation in the case of uncertain inputs. We
consider an emulator setup for $\mathbf{f}$ similar to that discussed in
Section 5.1. Following Equation (62), we choose to decompose the vector of
training runs $\mathbf{F}$ as follows:

$\mathbf{F}=\textrm{vec}(\bm{G}\bm{B})+\mathbf{U}=\bm{W}\bm{\beta}+\mathbf{U},$

where
$\bm{G}=(\mathbf{g}(\mathbf{x}^{(1)}),\dots,\mathbf{g}(\mathbf{x}^{(n)}))^{T}$
is an $n\times m$ matrix of regressors at the known design points in
$\mathcal{X}$, $\bm{W}=\bm{I}_{q}\otimes\bm{G}$, $\bm{I}_{q}$ is a $q\times q$
identity matrix, $\otimes$ represents the kronecker product,
$\mathbf{U}=\mathbf{u}(\mathcal{X})$ is an $nq$ vector of residuals, and
recall $\bm{\beta}=\textrm{vec}(\bm{B})$ with prior specification ${\rm
E}[\bm{\beta}]=\bm{\Gamma}$ and ${\rm Var}[\bm{\beta}]=\bm{\mathit{\Delta}}$.

We wish to make inference about $\mathbf{f}(\mathbf{X})$, where
$\mathbf{X}\in\mathbb{R}^{p}$ is an uncertain (random variable) input to
$\mathbf{f}$. Following Equation (62), $\mathbf{f}(\mathbf{X})$ can be written
as

$\mathbf{f}(\mathbf{X})=\mathbf{g}(\mathbf{X})^{T}\bm{B}+\mathbf{u}(\mathbf{X})=\bm{w}(\mathbf{X})\bm{\beta}+\mathbf{u}(\mathbf{X}),$

where $\bm{w}(\mathbf{X})=\bm{I}_{q}\otimes\mathbf{g}(\mathbf{X})^{T}$. We
assume ${\rm E}[\mathbf{u}(\mathbf{X})]=\mathbf{0}$ and ${\rm
Cov}[\bm{\beta},\mathbf{u}(\mathbf{X})]=\mathbf{0}$, these assumptions being a
reasonable extension to the prior specification that may be made in the case
of known inputs (Jackson, 2018). We also utilise the covariance structure for
$\mathbf{u}(\mathbf{X})$ as given by Equations (58) and (59). We therefore
have that ${\rm E}[\mathbf{U}]=\mathbf{0}$, ${\rm
Var}[\mathbf{U}]=\bm{\mathit{\Omega}}=\bm{\mathit{\Sigma}}\otimes\bm{C}$ and
${\rm Cov}[\bm{\beta},\mathbf{U}]=\bm{\mathit{0}}$, where we define

$\bm{C}=\left(\begin{array}[]{cccc}c(\mathbf{x}_{1},\mathbf{x}_{1})&c(\mathbf{x}_{1},\mathbf{x}_{2})&\cdots&c(\mathbf{x}_{1},\mathbf{x}_{n})\\\
c(\mathbf{x}_{2},\mathbf{x}_{1})&c(\mathbf{x}_{2},\mathbf{x}_{2})&\cdots&c(\mathbf{x}_{2},\mathbf{x}_{n})\\\
\vdots&\vdots&\ddots&\vdots\\\
c(\mathbf{x}_{n},\mathbf{x}_{1})&c(\mathbf{x}_{n},\mathbf{x}_{2})&\cdots&c(\mathbf{x}_{n},\mathbf{x}_{n})\end{array}\right).$

We also define
$\bm{v}(\mathbf{X})=\bm{\mathit{\Sigma}}\otimes\mathbf{c}(\mathbf{X})$ and
$\mathbf{c}(\mathbf{X})=(c(\mathbf{X},\mathbf{x}_{1}),\dots,c(\mathbf{X},\mathbf{x}_{n}))$.
We now proceed to state the adjusted belief formulae for
$\mathbf{f}(\mathbf{X})$ by $\mathbf{F}$ in the form of two lemmas, proofs of
which can be found in the supplementary material.

###### Lemma 5.2.1

The expected value of $\mathbf{f}(\mathbf{X})$, adjusted by $\mathbf{F}$, is
given by:

${\rm E}_{\mathbf{F}}[\mathbf{f}(\mathbf{X})]={\rm
E}[\bm{w}(\mathbf{X})]\,{\rm
E}_{\mathbf{F}}[\bm{\beta}]+\bm{v}(\mathbf{X})\,\bm{\mathit{\Omega}}^{-1}\,(\mathbf{F}-\bm{W}\,{\rm
E}_{\mathbf{F}}[\bm{\beta}])\,.$ (65)

###### Lemma 5.2.2

The variance of $\mathbf{f}(\mathbf{X})$, adjusted by $\mathbf{F}$, is given
by:

$\displaystyle{\rm Var}_{\mathbf{F}}[\mathbf{f}(\mathbf{X})]$ $\displaystyle=$
$\displaystyle{\rm E}[\bm{w}(\mathbf{X})\,{\rm
Var}_{\mathbf{F}}[\bm{\beta}]\,\bm{w}(\mathbf{X})^{T}]+{\rm
E}_{\mathbf{F}}[\bm{\beta}^{T}]\,{\rm Var}[\bm{w}(\mathbf{X})]\,{\rm
E}_{\mathbf{F}}[\bm{\beta}]+\,\bm{\mathit{\Sigma}}$ (66)
$\displaystyle\,\,\,\,\,\,\,-\,\bm{v}(\mathbf{X})\,\bm{\mathit{\Omega}}^{-1}\,\bm{v}(\mathbf{X})^{T}+\bm{v}(\mathbf{X})\,\bm{\mathit{\Omega}}^{-1}\,\bm{W}\,{\rm
Var}_{\mathbf{F}}[\bm{\beta}]\,\bm{W}^{T}\,\bm{\mathit{\Omega}}^{-1}\,\bm{v}(\mathbf{X})^{T}$
$\displaystyle\,\,\,\,\,\,\,-\,{\rm E}[\bm{w}(\mathbf{X})]\,{\rm
Var}_{\mathbf{F}}[\bm{\beta}]\,\bm{W}\,\bm{\mathit{\Omega}}^{-1}\,\bm{v}(\mathbf{X})^{T}$
$\displaystyle\,\,\,\,\,\,\,-\,({\rm E}[\bm{w}(\mathbf{X})]\,{\rm
Var}_{\mathbf{F}}[\bm{\beta}]\,\bm{W}\,\bm{\mathit{\Omega}}^{-1}\,\bm{v}(\mathbf{X})^{T})^{T}\,.$

Specification of ${\rm E}[\bm{w}(\mathbf{X})]$ and ${\rm
Var}[\bm{w}(\mathbf{X})]$ is straight forward for first-order linear
regression functions. It is also possible for further functions of the input
components, but these transformed input components will require a sensible
second-order specification. Expressions for ${\rm E}_{\mathbf{F}}[\bm{\beta}]$
and ${\rm Var}_{\mathbf{F}}[\bm{\beta}]$ (extended derivations for which are
provided in the supplementary material) are given by

$\displaystyle{\rm E}_{\mathbf{F}}[\bm{\beta}]$ $\displaystyle=$
$\displaystyle(\bm{W}^{T}\,\bm{\mathit{\Omega}}^{-1}\,\bm{W}+\bm{\mathit{\Delta}}^{-1})^{-1}(\bm{W}^{T}\,\bm{\mathit{\Omega}}^{-1}\,\bm{W}\,\hat{\bm{\beta}}_{GLS}+\bm{\mathit{\Delta}}^{-1}\bm{\Gamma}),$
(67) $\displaystyle{\rm Var}_{\mathbf{F}}[\bm{\beta}]$ $\displaystyle=$
$\displaystyle(\bm{W}^{T}\,\bm{\mathit{\Omega}}^{-1}\,\bm{W}+\bm{\mathit{\Delta}}^{-1})^{-1},$
(68)

where

$\hat{\bm{\beta}}_{GLS}=\bm{I}_{q}\otimes(\bm{G}^{T}\,\bm{C}^{-1}\,\bm{G})^{-1}\bm{G}^{T}\,\bm{C}^{-1}\,\mathbf{F}.$
(69)

As for the common known input case, vague priors on $\bm{\beta}$, which we
define to mean that the eigenvectors of $\bm{\mathit{\Delta}}$ tend to
$\infty$ and thus the prior has negligible effect on the posterior, resulting
in ${\rm E}_{\mathbf{F}}[\bm{\beta}]=\hat{\bm{\beta}}_{GLS}$ and ${\rm
Var}_{\mathbf{F}}[\bm{\beta}]=(\bm{W}^{T}\,\bm{\mathit{\Omega}}^{-1}\,\bm{W})^{-1}$.

Assuming an appropriate correlation structure across any possible
specification of the set $\\{{\rm E}[\mathbf{X}],{\rm
E}[\mathbf{X}^{\prime}],{\rm Cov}[\mathbf{X},\mathbf{X}^{\prime}]\\}$, the
results of Lemmas 5.2.1 and 5.2.2 can be used to provide a second-order
approximation of the output (${\rm E}_{\mathbf{F}}[\mathbf{f}(\mathbf{X})]$
and ${\rm Var}_{\mathbf{F}}[\mathbf{f}(\mathbf{X})]$) of any simulator at
random variable input $\mathbf{X}$ for which a second-order belief
specification $\\{{\rm E}[\mathbf{X}],{\rm Var}[\mathbf{X}]\\}$ is itself
provided.

###  Emulation of Simulator Networks

An application of UIBLE is that of the Bayes linear emulation of simulator
networks.

Complex systems can often be most appropriately modelled as a network of
simpler component simulators, that together form a composite simulator of the
entire system of interest. One simple, but important, example from
epidemiology, that will be considered further in Section 5.4, combines an
atmospheric Anthrax dispersion simulator (Legrand et al., 2009), labelled
$d(\cdot)$, with a dose-response (DR) simulator (Groer, 1978), labelled
$\rho(\cdot)$, in a simple chain network; see Figure 1. The composite
dispersion dose-response (DDR) simulator $h=\rho(d(\mathbf{z}))$ models the
overall process, where $\mathbf{z}$ can be viewed as the input to $d$ or $h$.

Utilising simulator networks for uncertainty quantification is challenging,
largely due to the variety of sources of uncertainty for each individual
component simulator (Goldstein et al., 2013) and the necessity of propagating
that uncertainty through the network. The issues surrounding computational
expense, discussed in the introductory paragraph to Section 5, are also
exacerbated, requiring use of emulation. A key question in such scenarios is
whether combining emulators for the component simulators within a network can
result in more powerful approximations than emulating the network as a single
composite simulator.

The output of a Bayes linear emulator, or UIBLE, as given by Equations (63)
and (64), or (65) and (156) respectively, is a second order belief
specification across the output components of $\mathbf{f}(\mathbf{x})$, which
can subsequently be used as the second order belief specification of the input
to a subsequent simulator or emulator. As a result, UIBLE can be used to
approximate arbitrarily large networks of simulators, where the random input
$\mathbf{X}$ to one simulator is taken to have a second-order belief
specification arising from a previous emulator.

Emulation of simulator networks has been addressed elsewhere in the
literature. Kyzyurova et al. (2018) proposed coupling two simulators by
linking independently developed Gaussian process emulators of the simulators.
Their motivation arose from potentially having separate training runs for two
simulators $\mathbf{f}^{1}=$ _bent_ and $\mathbf{f}^{2}=$ _puff_ , where
_bent_ simulates volcanic ash plumes arising from a vent and _puff_ simulates
ash dispersion. As a result, direct emulation (similar to that discussed in
Section 5.1) of the composite simulator defined by the chain is not possible.
For specific Gaussian process forms, they derived closed form expressions for
the overall mean and variance arising from linking the two emulators and
applied these quantities within a normal distribution approximation for the
composite emulator. Subsequently, Ming & Guillas (2021) extended the work of
emulating coupled simulators to much larger networks, also deriving closed-
form mean and variance expressions when using a class of Matérn correlation
functions. The availability of second-order posterior statistics for the
emulation of simulator networks provided in this literature, but the lack of a
closed-form distribution, naturally suggests the use of Bayes linear
approaches in the same context, hence the application of UIBLE.

Emulation, particularly Gaussian process emulation, of simulator networks, can
also be compared to emulation using deep Gaussian processes (Damianou &
Lawrence, 2013; Dunlop et al., 2018). Deep Gaussian processes arise from
belief networks about simulator behaviour based on Gaussian process mappings,
such that layers of Gaussian process latent variables exist between simulator
input and output, these being marginalised variationally (see, for example,
Titisias & Lawrence, 2010). Whilst similar, the intermediate variables in a
simulator network represent physical system properties, which aids the
construction and modelling of the emulators for each of the component
processes. Direct use of a deep Gaussian process for the entire network will
not exploit this additional information. Linking deep Gaussian processes of
component simulators can make use of the advantages of both observable and
latent variables (Ming et al., 2022), however, extension of such strategies to
the Bayes linear paradigm, where one is usually concerned with belief
specification over observable quantities, is an area for future research.

###  Application to a DDR Chain of Simulators

In this section, we apply UIBLE to an important example from epidemiology
which combines an atmospheric Anthrax dispersion simulator (Legrand et al.,
2009), labelled $d(\cdot)$, with a dose-response (DR) simulator (Groer, 1978),
labelled $\rho(\cdot)$, in a simple chain network; see Figure 1.

wind speed wind direction source mass
$\mathbf{z}$$d(\mathbf{z})$$\rho(d(\mathbf{z}))$ casualty proportion $d$$\rho$
Figure 1: Graphical representation of the DDR network of simulators
$h(\mathbf{z})=\rho(d(\mathbf{z}))$.

The dispersion simulator models the spread of a released biological agent
across a given spatial domain, with input parameters of interest corresponding
to physical quantities wind speed ($z_{WS}$), wind direction ($z_{WD}$) and
source mass ($z_{SM}$). Simulator outputs $d(\mathbf{z})$ represent dose at
each location across the domain, however, we consider a single spatial
location output for illustrative purposes. Due to the behaviour of $d(\cdot)$,
we chose to emulate a transformation of the output, namely
$f^{1}(\mathbf{x}_{1})=\log(d(\mathbf{x}_{1})+1)$, treating this transformed
function $f^{1}(\cdot)$ as the first simulator of the network.

The DR simulator takes dose, $x$, as input, and outputs casualties, $\rho(x)$,
as a proportion of the population at a particular spatial location. However,
to be consistent with $f^{1}(\cdot)$, we consider the second simulator to be
$f^{2}(x_{2})=\rho(\exp(x_{2})-1)$, so that
$h(\mathbf{z})=f^{2}(f^{1}(\mathbf{z}))=\rho(d(\mathbf{z}))$. We also note
that whilst $d(\cdot)$ is computationally expensive, dose-response model
$\rho(\cdot)$ is not; however we emulate both simulators to demonstrate the
efficacy of our methods. Our methods are also applicable and effective when
only a subset of the simulators in a network require emulation. As
$f^{2}(\cdot)$ here is straightforward to emulate, our application also
effectively demonstrates the use of the methods for this special case. The
simplicity of this second simulator additionally highlights a key benefit of
the emulation of simulator networks methodology in general, since if the
composite simulator is emulated as a single model, neither the computational
efficiency of running the second simulation nor the simplicity of the
behaviour of this simulator can be exploited.

The composite simulator $h=f^{2}\cdot f^{1}=\rho\cdot d$ takes wind speed,
wind direction and source mass as input $\mathbf{z}$, and directly outputs a
proportion of casualties $h(\mathbf{z})$. An expanded DAG showing Figure 1
shows the links between the original simulators $d$ and $\rho$, their inputs,
output and corresponding physical quantities.

The primary interest of decision makers is the impact of release conditions on
the proportion of casualties. This assessment requires linking the two
component simulators, each of which implements modelling from two different
groups of experts. To address the question of whether emulating individual
simulators in the network is preferable to a single emulator of the composite
simulator, we will compare direct Bayes linear emulation of $h(\cdot)$ with
application of UIBLE to $f^{2}(\cdot)$ given standard Bayes linear emulation
of $f^{1}(\cdot)$.

We proceeded to construct Bayes linear emulators for each of the component
simulators $f^{1}$ and $f^{2}$, as well as the composite simulator $h$. Ranges
of interest of the inputs to simulator $f^{1}$ (and thus $h$) are
$(z_{WD},z_{WS},z_{SM})\in[37,63]^{\circ}\times[1,150]\textrm{ms}^{-1}\times[0.001,1]\textrm{kg}$,
each of which were scaled to $[-1,1]$ for the purposes of our analysis. We
constructed a training point design for $f^{1}$ and $h$ using a maximin Latin
hypercube (McKay et al., 1979; Currin et al., 1991) of size 50 across the
three input dimensions. In contrast, simulator $f^{2}$ is one-dimensional,
thus the need for fewer training points, so we take a random sample of 20
points from a uniform distribution. This idea of using independent space-
filling designs for the construction of each of the two emulators is similar
to that proposed in Kyzyurova et al. (2018). More sophisticated design
strategies are discussed in Section 6, however, yield much scope for future
research. In order to explore the effects of design size, we subsequently
repeated the analysis using 30 training points for each simulator, this
following the ad-hoc $10d$ rule-of-thumb (Loeppky et al., 2009). It is correct
to use 30 (as opposed to 10) training points even for $f^{2}$, since it is the
power of the component-wise emulation of simulator networks that permits
separation of the 1-dimensional emulator from the 3-dimensional composite
simulator. For each of the emulators for $f^{1},f^{2}$ and $h$, we assumed a
Gaussian correlation function, as given by Equation (57), along with a first-
order polynomial mean function. We represent the scalar variance parameter and
correlation length vectors as $\sigma^{2}_{1},\sigma^{2}_{2},\sigma^{2}_{h}$
and $\bm{\theta}_{1},\theta_{2},\bm{\theta}_{h}$ respectively. We fit these
parameters using maximum likelihood for each emulator, this permitting a fair
comparison between the emulation methods presented. It should be noted that
whilst it can be argued that use of maximum likelihood lies outside of the
Bayes linear paradigm in which we present the methodology of this article, it
provides a useful tool for obtaining valid emulators. Alternative approaches,
such as making use of cross-validation, might also be used.

Figure 2: Adjusted expectation (red points) $\pm 3$ standard deviations for
each of the diagnostic points (indexed along the x-axis according to
increasing size of the simulator output) when 50 simulator training points
were used for $f^{1}$, and $h$, and 20 were used for $f^{2}$, with the
simulator output value given as a green point if falling within the $\pm 3$
standard deviation error bar, and black otherwise. Figure 3: Adjusted
expectation (red points) $\pm 3$ standard deviations for each of the
diagnostic points (indexed along the x-axis according to increasing size of
the simulator output) when 30 simulator training points were used for $f^{1}$,
$f^{2}$ and $h$, with the simulator output value given as a green point if
falling within the $\pm 3$ standard deviation error bar, and black otherwise.

Given the component emulators for $f^{1}$ and $f^{2}$, we can then combine
them using UIBLE to yield chained emulators for $h$. Figures 2 and 3 show
plots, for the case of 50 and 30 training points respectively, of adjusted
expectation (red points) $\pm 3$ standard deviations for each of a set of
diagnostic points (indexed along the x-axis according to increasing size of
the simulator output), with the simulator output value given as a green point
if falling within the $\pm 3$ standard deviation error bar, and black
otherwise. The input designs for these diagnostic runs were constructed in the
same manner as the training run designs.

In addition to the plots, Table 1 shows the Mean Absolute Standardised
Prediction Error (MASPE):

$\frac{1}{n}\sum_{k=1}^{n}\frac{|f(\mathbf{x}^{(k)})-\mu_{f}(\mathbf{x}^{(k)})|}{\sqrt{\nu_{f}(\mathbf{x}^{(k)})}},$
(70)

Root Mean Squared Prediction Error (RMSPE) (Bastos & O’Hagan, 2008):

$\sqrt{\frac{1}{n}\sum_{k=1}^{n}(f(\mathbf{x}^{(k)})-\mu_{f}(\mathbf{x}^{(k)}))^{2}},$
(71)

and Mean Generalised Entropy Score (MGES), as defined by Equation (27) of
Gneiting & Raftery (2007):

$-\,\frac{1}{n}\sum_{k=1}^{n}\left\\{\frac{\left[f(\mathbf{x}^{(k)})-\mu_{f}(\mathbf{x}^{(k)})\right]^{2}}{\nu_{f}(\mathbf{x}^{(k)})}+\log(\nu_{f}(\mathbf{x}^{(k)}))\right\\}$
(72)

for the diagnostic runs for each of the four emulators, with $\mu_{f}$,
$\nu_{f}$ representing appropriate mean and variance estimators corresponding
to generic simulator output $f$. MASPE is a measure of emulator validity;
heuristically we expect this value to be roughly 1 (assuming normal errors
this value should be $\sqrt{2/\pi}$). RMSPE permits comparison of emulator
accuracy. MGES is larger (better) for approximations that are both valid and
precise.

Table 1: The MASPE, RMSPE and MGES for each of the four approximations and two
design sizes discussed in Section 5.4. MASPE and RMSPE are smaller-the-better
quantities; MGES is larger-the-better.

DE of $f^{1}$ DE of $f^{2}$ DE UIBLE MASPE 0.850 0.580 1.676 0.485 50-50-20
RMSPE 0.520 0.00031 0.0308 0.0223 MGES 0.343 15.161 3.960 6.616 MASPE 0.798
0.987 2.075 0.754 30-30-30 RMSPE 0.616 0.00025 0.0534 0.0379 MGES 0.111 15.176
-0.019 5.683

We can see from Figure 2a that the emulator for $f^{1}$ is fairly accurate,
with less accuracy where simulator output lies towards the lower end of the
output range. The emulator for $f^{2}$ (Figure 2b) is very accurate,
reflecting the fact that emulator predictions can be taken with almost as much
certainty as running the simulator itself. The direct emulator for $h$ (Figure
2c) yields predictions with underestimated uncertainty. By comparison, the
estimated uncertainty for UIBLE is perhaps slightly overestimated, but in
general yields both more appropriate MASPE values and improved MGES values. In
addition, the accuracy of the predictions for the UIBLE are, on the whole,
improved, this being confirmed by the RMSPE values.

Looking now to Figure 3, we can see greater estimated uncertainty and lower
levels of accuracy for the direct emulators of $f^{1}$ and $h$, as well as the
UIBLE, relative to the corresponding predictions made using the emulators
trained with 50 simulator runs. What is consistent for both design sizes,
however, is that direct emulation underestimates the uncertainty towards the
upper end of the input space, where the predictions are in addition more
inaccurate than those resulting from UIBLE.

It is particularly noticable, in both Figures 2d and 3d, that UIBLE seems to
overestimate uncertainty, particularly where the simulator output lies towards
the lower end of the output range. This is likely to be a consequence of the
uncertainty propagation approach. UIBLE has uncertainty from the regression
part and covariance structure. As with standard Bayes linear emulation that
uses a single correlation structure across $\mathbb{X}$ , there is some
averaging of the uncertainty estimates for simulator prediction across the
input space, even if the behaviour at some points is smoother than others. Use
of uncertain input adaptations to non-stationary covariance structures could
address this, and is an area for future research.

To summarise, we feel that the results presented give evidence for utilising
UIBLE for linking emulators of component simulators in a network over using a
direct emulator of the composite simulator in many cases.

##  6 Conclusion

We have presented novel general Bayes linear methodology for the analysis of
statistical models with uncertain inputs. This has been achieved by extending
commonly-employed second-order modelling assumptions to be applicable in the
context of uncertain inputs.

Section A.2 demonstrated the applicability of the methodology within the
context of a regression model with correlated error structure in the case of
known inputs. By considering adaptations of the constituent components of the
prior correlation structure, we developed an appropriate prior covariance
structure which could be adjusted under observations (at uncertain inputs) of
the modelled electrolysis extraction process over time. An important feature
of the presented adapted covariance structures is the extended homogeneity of
the error structure across possible specifications of the set $\\{{\rm
E}[\mathbf{X}],{\rm E}[\mathbf{X}^{\prime}],{\rm
Cov}[\mathbf{X},\mathbf{X}^{\prime}]\\}$. This homogeneity may result in
conservative approximations in the sense of effectively overestimating
variances ${\rm Var}[\mathbf{y}(\mathbf{X})]$, and underestimating covariances
${\rm Cov}[\mathbf{y}(\mathbf{X}),\mathbf{y}(\mathbf{X}^{\prime})]$ for
$\mathbf{X}\neq\mathbf{X}^{\prime}$, relative to the expected values of the
corresponding covariances, that is ${\rm
E}_{\mathbf{X},\mathbf{X}^{\prime}}[{\rm
Cov}[\mathbf{y}(\mathbf{X}),\mathbf{y}(\mathbf{X}^{\prime})]]$, were we to
assume a probability distribution over $\mathbf{X}$ and $\mathbf{X}^{\prime}$
and integrate out. Such approximations can enter either by conservatively
optimising over introduced parameters (such as $p$ in Equation (103) whilst
establishing ${\rm Cov}[J(T),J(T^{\prime})]$), or by making direct statements
via inequalities (such as in (117) whilst establishing ${\rm
Cov}[H(T),H(T^{\prime})]$). Whilst practically making little difference, it is
perhaps important to bear in mind that whilst the former approximation
resolves to the expression for ${\rm Cov}[J(t),J(t^{\prime})]$ in the case
that ${\rm Var}[T]={\rm Var}[T^{\prime}]=0$, the latter approximation does not
have such resolution.

We justify these covariance structures as useful when we are unwilling or
unable to meaningfully specify a probability distribution over the inputs
$\mathbf{X}$ and $\mathbf{X}^{\prime}$. Whilst perhaps deemed conservative,
there are many statistical models for which it can be viewed as preferable to
underestimate resolved uncertainty in a Bayes linear paradigm than to
overestimate resolved uncertainty by arbitrary specification of too
concentrated a probability distribution over the inputs, leading to an
underestimation of adjusted uncertainty for future inputs, perhaps leading to
real-world decision-making with dire consequences. Having said this, less
conservative covariance structures may be available via exploration of
heterogeneous covariance structures over the space of possible specifications
$\\{{\rm E}[\mathbf{X}],{\rm E}[\mathbf{X}^{\prime}],{\rm
Cov}[\mathbf{X},\mathbf{X}^{\prime}]\\}$ for $\mathbf{X},\mathbf{X}^{\prime}$,
this being an area for future research.

In Section 5, we developed UIBLE, which is a direct extension of Bayes linear
emulation methodology to the case of training emulators with, and predicting
for, simulator output from uncertain inputs. UIBLE is a computationally
efficient approximation of simulator output given a second-order belief
specification regarding the input, with each evaluation being akin to a run of
a standard Bayes linear emulator. The modelling assumptions are pragmatic, but
the resulting emulator can, and should, be assessed using diagnostic summaries
and plots for overall adequacy as anyway necessary for a standard Bayes linear
emulator (or indeed any approximating statistical model).

The specific application example which we presented was for component-wise
emulation of a chain or network of simulators. In this case, the input to the
emulator for the dose-response simulator was uncertain as a result of it being
the emulated output of the dispersion simulator. Such networks of simulators
are common, both explicitly as here, and given the fact that many complex
simulators, for example climate simulators, are usually made up of smaller
components which arguably form a network. The potential benefits in many
scenarios of linking emulators of component-simulators in contrast to direct
emulation of the composite simulator is exemplified through the results of
Section 5.4, along with the results of Kyzyurova et al. (2018) and Ming &
Guillas (2021), albeit with these works focusing on Gaussian process emulation
rather than Bayes linear emulation. In some cases, the links between the
simulators in a network, and particularly those between the outputs of an
earlier simulator and the inputs of a later simulator, may not be directly
compatible. In such cases, these differences are a potential source of model
discrepancy, which should be taken into account, and an area for future
research, however, we have assumed a direct correspondence within this article
to aid the illustration of our methodology.

There are future directions for UIBLE in the context of the emulation of a
network of simulators. One example would be the consideration of suitable
adaptations to additional covariance structures, such as the Matérn
correlation function (Matern, 1947; Paulo, 2005), or additional adaptations
for correlation functions used for the emulation of stochastic simulators with
uncertain inputs.

In this work, the two emulators were constructed independently, before being
linked together to form an approximation to the composite simulators. However,
both design and parameter estimation for all component emulators could be
considered simultaneously, and provide much scope for future research. The
issue regarding design has been partially addressed in Ming & Guillas (2021),
where it is correctly highlighted that ignorance of structural dependence
cause unnecessary refinements of component emulators that are insignificant to
the global output. They propose a variance-based adaptive design, which
selects at each iteration first the component emulator with greatest average
variance contribution to the linked emulator, before selecting the single
design point for this emulator based on individual contribution to the linked
emulator. Whilst shown to be effective in comparison to the non-sequential
design strategy of Kyzyurova et al. (2018), it should be noted that the
emulator with greatest average variance contribution may not contain the
optimal simulator run across all component simulators for reducing the
targeted aspects of linked emulator variance, for example, average or maximal
variance.

A similar design approach based on predictive means in the Bayes linear
paradigm, which may get around this issue, is outlined as follows. For each
simulator $f^{i}$, locate input $\mathbf{x}_{i}$ with (near-)optimal predicted
composite emulator variance reduction by constructing the relevant component
emulator under the assumption that $f^{i}(\mathbf{x}_{i})={\rm
E}[f^{i}(\mathbf{x}_{i})]$, and hence ${\rm Var}[f^{i}(\mathbf{x}_{i})]=0$,
for a large range of possible $\mathbf{x}_{i}$. Exploration over possible
$\mathbf{x}_{i}$ may involve a large space-filling design, or more likely an
efficient optimisation algorithm, for each component emulator. Once the
optimal point for each simulator is found, select the simulator for which the
component optimal composite emulator variance reduction is greatest. Whilst we
postulate this design algorithm, we do not apply it here because it is not the
main focus of this work, and thorough exploration into this topic, which makes
advances on this suggestion and the algorithm of Ming & Guillas (2021), is a
separate research area requiring thorough attention. It should be noted in
particular that both these algorithms neglect the computational efficiency of
the component simulators, which are likely to vary significantly. A realistic
and applicable design algorithm, assuming it to be restricted by overall
computational budget, must take this into account.

In addition to situations involving networks of simulators, UIBLEs have more
general application. For example, the inputs of a simulator may be uncertain
as a result of wishing to make predictions about the corresponding system
under a specific scenario, but not knowing the direct relationship between the
real-world quantities representing the scenario and the corresponding
parameter values in the model. Alternatively, UIBLEs would permit efficient
sensitivity analyses; several evaluations of an emulator with constant ${\rm
E}[\mathbf{X}]$ and varying ${\rm Var}[\mathbf{X}]$ could quickly provide an
indication of the influence of individual inputs on simulator output
behaviour.

More generally, the methodology presented in this article may fit smoothly
within the context of Bayes linear Bayes graphical models and Bayes linear
kinematics (Goldstein & Shaw, 2004). In particular, both approaches involve
specification and adjustment of our beliefs about uncertain quantities within
a Bayes linear structure based on revised, or updated, beliefs about other
quantities within the structure. We view this connection as an area for future
research, since the inputs to the statistical model over which we are
uncertain may enter into the modelling process in a non-linear way, for
example, through use of correlation functions. In addition, we wish the
specification, of model output, to hold for all possible specifications
$\\{{\rm E}[\mathbf{X}],{\rm Var}[\mathbf{X}]\\}$ of the input variable
$\mathbf{X}$ that we may have.

To conclude, this work opens up many doors for the analysis of statistical
modelling with uncertain inputs, particular within a Bayes linear paradigm.
The most appropriate belief structures for different practical applications
will be situation-specific, and scope for future research themselves, but will
benefit from the foundations laid out within this article.

### Acknowledgements

This work was supported by Chemical and Biological Technologies Department
(contract HDTRA1-17-C-0028). We are grateful to Crystalcast project members
for invaluable discussions, comments, and provision of the simulators for the
dispersion dose-response application. Particular thanks are due to Professor
Veronica Bowman and Dr Daniel Silk (Defence Science and Technology Laboratory,
UK), and Dr Daria Semochkina (University of Southampton, UK).

## References

  * (1)
  * Astfalck et al. (2021) Astfalck, L., Williamson, D., Gandy, N., Gregoire, L. & Ivanovic, R. (2021), Coexchangeable process modelling for uncertainty quantification in joint climate reconstruction.
  * Bastos & O’Hagan (2008) Bastos, T. S. & O’Hagan, A. (2008), ‘Diagnostics for Gaussian process emulators.’, Technometrics 51, 425–438.
  * Craig et al. (1996) Craig, P. S., Goldstein, M., Seheult, A. H. & Smith, J. A. (1996), Bayes linear strategies for matching hydrocarbon reservoir history, in J. M. Bernardo, J. O. Berger, A. P. Dawid & A. F. M. Smith, eds, ‘Bayesian Statistics 5’, Clarendon Press, Oxford, pp. 69–95.
  * Cumming & Goldstein (2009) Cumming, J. & Goldstein, M. (2009), ‘Bayes linear uncertainty analysis for oil reservoirs based on multiscale computer experiments’, Technometrics 51(4), 377–388.
  * Currin et al. (1991) Currin, C., Mitchell, T., Morris, M. & Ylvisaker, D. (1991), ‘Bayesian prediction of deterministic functions with applications to the design and analysis of computer experiments’, Journal of the American Statistical Association 86(416), 953–963.
  * Dallaire et al. (2009) Dallaire, P., Besse, C. & Chaib-Draa, B. (2009), ‘Learning gaussian process models from uncertain data’, International Conference on Neural Information Processing pp. 433–440.
  * Damianou & Lawrence (2013) Damianou, A. C. & Lawrence, N. D. (2013), ‘Deep gaussian processes’, Proceedings of the 16th International Conference on Artificial Intelligence and Statistics 31.
  * de Finetti (1974) de Finetti, B. (1974), Theory of Probability, Vol. 1, Wiley.
  * de Finetti (1975) de Finetti, B. (1975), Theory of Probability, Vol. 2, Wiley.
  * Dunlop et al. (2018) Dunlop, M. M., Girolami, M. A., Stuart, A. M. & Teckentrup, A. L. (2018), ‘How deep are deep gaussian processes’, Journal of Machine Learning Research 19, 1–46.
  * Durbin (1954) Durbin, J. (1954), ‘Errors in variables’, Review of the International Statistical Institute 22(1/3), 23–32.
  * Figueroa-Zuniga et al. (2022) Figueroa-Zuniga, J. I., Bayes, C. L., Leiva, V. & Liu, S. (2022), ‘Robust beta regression modeling with errors-in-variables: a bayesian approach and numerical applications’, Statistical Papers 63, 919–942.
  * Garthwaite et al. (2005) Garthwaite, P., Kadane, J. & O’Hagan, A. (2005), ‘Statistical methods for eliciting probability distributions’, Journal of the American Statistical Association 100, 680–701.
  * Gneiting & Raftery (2007) Gneiting, T. & Raftery, A. E. (2007), ‘Stricitly proper scoring rules, prediction and estimation’, Journal of the American Statistical Association 102(477), 359–378.
  * Goldstein (1999) Goldstein, M. (1999), Bayes linear analysis, in S. Kotz, C. B. Read, N. Balakrishnan & B. Vidakovic, eds, ‘Encyclopedia of statistical Sciences’, Wiley, New York, chapter Bayes Linear Analysis, pp. 29–34.
  * Goldstein & Huntley (2016) Goldstein, M. & Huntley, N. (2016), Bayes Linear Emulation, History Matching, and Forecasting for Complex Computer Simulators, Springer International Publishing, Cham, pp. 1–24.
https://doi.org/10.1007/978-3-319-11259-6_14-1

  * Goldstein & Rougier (2004) Goldstein, M. & Rougier, J. C. (2004), ‘Probabilistic formulations for transferring inferences from mathematical models to physical systems’, SIAM Journal on Scientific Computing 26(2), 467–487.
  * Goldstein et al. (2013) Goldstein, M., Seheult, A. & Vernon, I. (2013), Assessing model adequacy, in J. Wainwright & M. Mulligan, eds, ‘Environmental Modelling: Finding Simplicity in Complexity’, John Wiley and Sons, Chichester.
  * Goldstein & Shaw (2004) Goldstein, M. & Shaw, S. (2004), ‘Bayes linear kinematics and Bayes linear Bayes graphical models’, Biometrika 91(2), 425–446.
  * Goldstein & Wooff (2007) Goldstein, M. & Wooff, D. (2007), Bayes Linear Statistics, Wiley, Chichester.
  * Goldstein & Wooff (1998) Goldstein, M. & Wooff, D. A. (1998), ‘Adjusting exchangeable beliefs’, Biometrika 85(1), 39–54.
  * Gosling et al. (2013) Gosling, J. P., Hart, A., Owen, H., Davies, M., Li, J. & MacKay, C. (2013), ‘A Bayes linear approach to weight-of-evidence risk assessment for skin allergy’, Bayesian Analysis 8(1), 169–186.
  * Groer (1978) Groer, P. G. (1978), ‘Dose-response curves and competing risks’, Proceedings of the National Academy of Sciences of the United States of America 75(9), 4087–4091.
  * Hartigan (1969) Hartigan, J. A. (1969), ‘Linear Bayesian methods’, Journal of the Royal Statistical Society 31, 446–454.
  * Harville (2018) Harville, D. A. (2018), Linear Models and the relevant distributions and matrix algebra, CRC press, Boca Raton.
  * Jackson (2018) Jackson, S. E. (2018), Design of Physical System Experiments Using Bayes Linear Emulation and History Matching Methodology with Application to Arabidopsis Thaliana, PhD thesis, Durham University.
  * Kennedy & O’Hagan (2001) Kennedy, M. C. & O’Hagan, A. (2001), ‘Bayesian calibration of computer models’, Journal of the Royal Statistical Society 63(3), 425–464.
  * Koehler & Owen (1996) Koehler, J. R. & Owen, A. B. (1996), ‘Computer experiments’.
  * Kyzyurova et al. (2018) Kyzyurova, K. N., Berger, J. O. & Wolpert, R. L. (2018), ‘Coupling computer models through linking their statistical emulators’, Journal on Uncertainty Quantification 6(3), 1151–1171.
  * Legrand et al. (2009) Legrand, J., Egan, J. R., Hall, I. M., Cauchemez, S., Leach, S. & Ferguson, N. M. (2009), ‘Estimating the location and spatial extent of a covert anthrax release’, PLoS Computational Biology 5(1).
  * Loeppky et al. (2009) Loeppky, J. L., Sacks, J. & Welch, W. J. (2009), ‘Choosing the sample size of a computer experiment: A practical guide’, Technometrics 51(4), 366–376.
  * Mardia et al. (1979) Mardia, K. V., Kent, J. T. & Bibby, J. M. (1979), Multivariate Analysis, London.
  * Matern (1947) Matern, B. (1947), ‘Methods of estimating the accuracy of line and sample plot surveys’, Meddelelser fra Statens Skogsforskningsinstitut 36.
  * McHutchon & Rasmussen (2011) McHutchon, A. & Rasmussen, C. E. (2011), ‘Gaussian process training with input noise’, Advances in Neural Information Processing Systems 24.
  * McKay et al. (1979) McKay, M. D., Beckman, R. J. & Conover, W. J. (1979), ‘A comparison of three methods for selecting values of input variables in the analysis of output from a computer code’, Technometrics 21(2), 239–245.
  * Ming & Guillas (2021) Ming, D. & Guillas, S. (2021), ‘Linked gaussian process emulation for systems of computer models using matern kernels and adaptive design’, Journal of Uncertainty Quantification 9(4), 1615–1642.
  * Ming et al. (2022) Ming, D., Williamson, D. & Guillas, S. (2022), Deep gaussian process emulation using stochastic imputation.
  * O’Hagan (1987) O’Hagan, A. (1987), ‘Bayes linear estimators for randomized response models’, Journal of the American Statistical Association 82, 580–585.
  * Paulo (2005) Paulo, R. (2005), ‘Default priors for Gaussian processes’, The Annals of Statistics 33(2), 556–582.
  * Shawe-Taylor & Cristianini (2011) Shawe-Taylor, J. & Cristianini, N. (2011), Kernel Methods for Pattern Analysis, Cambridge University Press, Cambridge.
  * Titisias & Lawrence (2010) Titisias, M. K. & Lawrence, N. D. (2010), ‘Bayesian gaussian process latent variable model’, Proceedings of the 13th International Conference on Artificial Intelligence and Statistics 9, 844–851.
  * Vernon et al. (2010) Vernon, I., Goldstein, M. & Bower, R. G. (2010), ‘Galaxy formation: A Bayesian uncertainty analysis’, Bayesian Analysis 5(4), 619–669.
  * Wang et al. (2022) Wang, W., Yue, X., Haaland, B. & Wu, C. F. J. (2022), ‘Gaussian processes with input location error and applications to the composite parts assembly process’, Journal on Uncertainty Quantification 10(2).
  * Whittle (1992) Whittle, P. (1992), Probability Via Expectation, Springer.
  * Ye et al. (2022) Ye, D., Zun, P., Krzhizhanovskaya, V. & Hoekstra, A. (2022), ‘Uncertainty quantification of a three-dimensional in-stent restenosis model with surrogate modelling’, Royal Society Interface 18.

##  Appendix A Extended Calculations for Section 3

We here present the material covered in Section 3 of the main text again, but
with extended detail in calculation for the interested reader.

Consider a regression model, presented here with scalar output $y$ for
simplicity of notation and exposition:

$y(\mathbf{X})=\mathbf{X}^{T}\bm{\beta}+\epsilon(\mathbf{X})$ (73)

where $\bm{\beta}\in\mathbb{R}^{p}$ is a vector of regression parameters, and
$\mathbf{X}\in\mathbb{R}^{p}$ is a vector random variable.

We consider a generic prior specification over $\bm{\beta}$ of the form ${\rm
E}[\bm{\beta}]=\bm{\Gamma}$ and ${\rm Var}[\bm{\beta}]=\bm{\mathit{\Delta}}$.
Under the usual scenario of known $\mathbf{x}$, we are likely to specify ${\rm
E}[\epsilon(\mathbf{x})]=0$, thus we view it a reasonable extension to specify
${\rm E}[\epsilon(\mathbf{X})]=0$ in the random variable scenario. It is
common for $\epsilon$ to be viewed as uncorrelated with $\bm{\beta}$ _a
priori_ , that is to specify ${\rm Cov}[\bm{\beta},\epsilon(\mathbf{x})]=0$,
so again we view it as reasonable to extend this in the random variable
context to ${\rm Cov}[\bm{\beta},\mathbf{X}]=0$.

Following the above prior specifications, adjusted belief specification for
$y(\mathbf{X})$ can be obtained as follows:

$\displaystyle{\rm E}_{\mathcal{Y}}[y(\mathbf{X})]$ $\displaystyle=$
$\displaystyle{\rm E}[y(\mathbf{X})]+{\rm
Cov}[y(\mathbf{X}),\mathcal{Y},\,]{\rm
Var}[\mathcal{Y}]^{-1}\,(\mathcal{Y}-{\rm E}[\mathcal{Y}]),$ (74)
$\displaystyle{\rm
Cov}_{\mathcal{Y}}\left[y(\mathbf{X}),y(\mathbf{X}^{\prime})\right]$
$\displaystyle=$ $\displaystyle{\rm
Cov}[y(\mathbf{X}),y(\mathbf{X}^{\prime})]-{\rm
Cov}[y(\mathbf{X}),\mathcal{Y}]\,{\rm Var}[\mathcal{Y}]^{-1}\,{\rm
Cov}[\mathcal{Y},y(\mathbf{X}^{\prime})],$ (75)

where

$\displaystyle{\rm E}[\mathcal{Y}]$ $\displaystyle=$ $\displaystyle({\rm
E}[y(\mathbf{X}^{(1)})],\dots,{\rm E}[y(\mathbf{X}^{(n)})]),$ (76)
$\displaystyle{\rm Var}[\mathcal{Y}]$ $\displaystyle=$ $\displaystyle\\{{\rm
Cov}[y(\mathbf{X}^{(i)}),y(\mathbf{X}^{(j)})]\\}_{i,j=1,\dots,n}.$ (77)

Therefore the required prior specifications are ${\rm E}[y(\mathbf{X})]$ and
${\rm Cov}[y(\mathbf{X}),y(\mathbf{X}^{\prime})]$ for any
$\mathbf{X},\mathbf{X}^{\prime}\in\mathbb{X}$.

We have that:

$\displaystyle{\rm E}[y(\mathbf{X})]$ $\displaystyle=$ $\displaystyle{\rm
E}[\mathbf{X}^{T}\bm{\beta}]+{\rm E}[\epsilon(\mathbf{X})]$ (78)
$\displaystyle=$ $\displaystyle{\rm E}[\mathbf{X}^{T}]\,\bm{\Gamma},$

and

$\displaystyle{\rm Cov}[y(\mathbf{X}),y(\mathbf{X}^{\prime})]$
$\displaystyle=$ $\displaystyle{\rm
Cov}[\mathbf{X}^{T}\bm{\beta}+\epsilon(\mathbf{X}),\mathbf{X}^{\prime
T}\bm{\beta}+\epsilon(\mathbf{X}^{\prime})]$ (79) $\displaystyle=$
$\displaystyle{\rm Cov}[\mathbf{X}^{T}\bm{\beta},\mathbf{X}^{\prime
T}\bm{\beta}]+{\rm Cov}[\epsilon(\mathbf{X}),\epsilon(\mathbf{X}^{\prime})]$
$\displaystyle\,\,+\,{\rm
Cov}[\mathbf{X}^{T}\bm{\beta},\epsilon(\mathbf{X}^{\prime})]+{\rm
Cov}[\epsilon(\mathbf{X}),\mathbf{X}^{\prime T}\bm{\beta}]$ $\displaystyle=$
$\displaystyle{\rm Cov}[\mathbf{X}^{T}\bm{\beta},\mathbf{X}^{\prime
T}\bm{\beta}]+{\rm Cov}[\epsilon(\mathbf{X}),\epsilon(\mathbf{X}^{\prime})],$

with the final line of Equation (79) resulting from the law of total
covariance:

$\displaystyle{\rm
Cov}[\mathbf{X}^{T}\bm{\beta},\epsilon(\mathbf{X}^{\prime})]$ (80)
$\displaystyle=$ $\displaystyle{\rm E}[{\rm
Cov}[\mathbf{X}^{T}\bm{\beta},\epsilon(\mathbf{X}^{\prime})|\mathbf{X}=\mathbf{x},\mathbf{X}^{\prime}=\mathbf{x}^{\prime}]]+{\rm
Cov}[{\rm E}[\mathbf{X}^{T}\bm{\beta}|\mathbf{X}=\mathbf{x}],{\rm
E}[\epsilon(\mathbf{X}^{\prime})|\mathbf{X}^{\prime}=\mathbf{x}^{\prime}]]$
$\displaystyle=$ $\displaystyle{\rm E}[\mathbf{X}^{T}\,{\rm
Cov}[\bm{\beta},\epsilon(\mathbf{X}^{\prime})]]+{\rm Cov}[{\rm
E}[\mathbf{X}^{T}\bm{\beta}|\mathbf{X}=\mathbf{x}],0]$ $\displaystyle=$
$\displaystyle 0+0=0.$

For the first term of Equation (79), we have, again using the law of total
covariance, that

$\displaystyle{\rm Cov}[\mathbf{X}^{T}\bm{\beta},\mathbf{X}^{\prime
T}\bm{\beta}]$ $\displaystyle=$ $\displaystyle{\rm E}[{\rm
Cov}[\mathbf{X}^{T}\bm{\beta},\mathbf{X}^{\prime
T}\bm{\beta}|\mathbf{X}=\mathbf{x},\mathbf{X}^{\prime}=\mathbf{x}^{\prime}]]+{\rm
Cov}[{\rm E}[\mathbf{X}^{T}\bm{\beta}|\mathbf{X}=\mathbf{x}],{\rm
E}[\mathbf{X}^{\prime T}\bm{\beta}|\mathbf{X}^{\prime}=\mathbf{x}^{\prime}]]$
$\displaystyle=$ $\displaystyle{\rm
E}[\mathbf{X}^{T}\bm{\mathit{\Delta}}\mathbf{X}^{\prime}]+\bm{\Gamma}^{T}{\rm
Cov}[\mathbf{X},\mathbf{X}^{\prime}]\bm{\Gamma},$ (82)

with

$\displaystyle{\rm E}[\mathbf{X}^{T}\bm{\mathit{\Delta}}\mathbf{X}^{\prime}]$
$\displaystyle=$ $\displaystyle{\rm
E}[\textrm{trace}(\mathbf{X}^{T}\bm{\mathit{\Delta}}\mathbf{X}^{\prime})]$
(83) $\displaystyle=$ $\displaystyle{\rm
E}[\textrm{trace}(\bm{\mathit{\Delta}}\mathbf{X}^{\prime}\mathbf{X}^{T})]$
$\displaystyle=$ $\displaystyle\textrm{trace}({\rm
E}[\bm{\mathit{\Delta}}\mathbf{X}^{\prime}\mathbf{X}^{T}])$ $\displaystyle=$
$\displaystyle\textrm{trace}(\bm{\mathit{\Delta}}{\rm
E}[\mathbf{X}^{\prime}\mathbf{X}^{T}])$ $\displaystyle=$
$\displaystyle\textrm{trace}(\bm{\mathit{\Delta}}({\rm
E}[\mathbf{X}^{\prime}]{\rm E}[\mathbf{X}^{T}]+{\rm
Cov}[\mathbf{X}^{\prime},\mathbf{X}])).$

The second term of Equation (79) is a covariance specification over the
residual process $\epsilon(\mathbf{X})$. In Sections A.1 and A.2, we present
two example model belief structures in the context of statistical modelling
with known inputs, demonstrating reasonable generalisation to the random
variable input scenario using the ideas presented in Section 2 of the main
text.

###  Linear Regression with Uncorrelated Random Error

In the common known input case, the standard linear regression model with
uncorrelated and homoscedastic random error can be obtained by specifying a
prior covariance structure over $\epsilon(\mathbf{x})$ of

${\rm
Cov}[\epsilon(\mathbf{x}),\epsilon(\mathbf{x}^{\prime})]=\mathbb{I}_{\mathbf{x}=\mathbf{x}^{\prime}}\sigma^{2},$
(84)

where $\mathbb{I}$ is the indicator function taking value $1$ if the statement
is true and $0$ otherwise, and $\sigma^{2}$ is a common scalar variance
parameter.

This specification can be simply extended to the random variable case as
follows:

${\rm
Cov}[\epsilon(\mathbf{X}),\epsilon(\mathbf{X}^{\prime})]=\mathbb{I}_{\mathbf{X}=\mathbf{X}^{\prime}}\sigma^{2},$
(85)

where the covariance across two random variable inputs is zero unless they are
known to be the same. Note that we could have two non-identical random
variables with the identical second order belief specification.

###  Example: Extracting Aluminium by Electrolysis

State that $t\geq 3$ in order to generalise to continuous $t$ values.

In this section, we consider an adaptation of the exchangeable regressions
example for extracting aluminium by electrolysis over time presented in
Goldstein & Wooff (1998). We first consider the covariance structure as
presented in that paper with known input, before presenting a reasonable
generalisation to the random variable input scenario using the ideas presented
in Section 2 of the main text.

The model is of the form

$y(t)=\beta_{0}+\beta_{1}\,t+\epsilon(t),$ (86)

and is exchangeable over $\bm{\beta}$. Note that this regression model is
consistent with the general form presented in Equation (73), with $p=2$
variables being an intercept and single controllable parameter $t$. In
Goldstein & Wooff (1998), it is assumed that $t=1,\dots,13$, as well as a
structured error model as follows:

$\displaystyle\epsilon(t)$ $\displaystyle=$ $\displaystyle A(t)+J(t)+H(t),$
(87) $\displaystyle J(t)$ $\displaystyle=$ $\displaystyle J(t-1)+Q(t),$
$\displaystyle H(t)$ $\displaystyle=$ $\displaystyle\psi\,H(t-1)+R(t),$

with $\psi\in(0,1)$, and where, for generality, we denote ${\rm
Var}[A(t)]=\sigma^{2}_{A}$, ${\rm Var}[Q(t)]={\rm Var}[J(1)]=\sigma^{2}_{Q}$,
${\rm Var}[R(t)]=\sigma^{2}_{R}$, and ${\rm Var}[H(1)]=\sigma^{2}_{1}$. These
terms express discrepancies from the linear trend as the sum of a pure
measurement error $A(t)$, a stochastic development of the discrepancy as a
random walk with drift $J(t)$, and an autoregressive term expressing the
measurement of the suspended particles in the chemical analysis $H(t)$. For
further details see Goldstein & Wooff (1998). Whilst the original model
assumes discrete timesteps $t$, we will generalise the regression structure to
continuous $t$-values, subject to $t\geq 3$ for consistency issues, before
addressing the uncertain input $T$ scenario.

In this example, we have that:

$\displaystyle{\rm Cov}[\epsilon(t),\epsilon(t)]$ $\displaystyle=$
$\displaystyle{\rm
Cov}[A(t)+J(t)+H(t),A(t^{\prime})+J(t^{\prime})+H(t^{\prime})]$ (88)
$\displaystyle=$ $\displaystyle{\rm Cov}[A(t),A(t^{\prime})]+{\rm
Cov}[J(t),J(t^{\prime})]+{\rm Cov}[H(t),H(t^{\prime})].$

We consider each term of Equation (88) separately. To begin with, we have

${\rm Cov}[A(t),A(t^{\prime})]=\mathbb{I}_{t=t^{\prime}}\sigma_{A}^{2}.$ (89)

Secondly, defining $d=|t-t^{\prime}|$, and $t_{m}=\min(t,t^{\prime})$, we have

$\displaystyle{\rm Cov}[J(t),J(t^{\prime})]$ $\displaystyle=$
$\displaystyle{\rm Cov}[J(t_{m}),J(t_{m}+d)]$ (90) $\displaystyle=$
$\displaystyle{\rm Cov}[J(t_{m}),J(t_{m})+Q_{t_{m}+1}+...+Q_{t_{m}+d}]$
$\displaystyle=$ $\displaystyle{\rm Var}[J(t_{m})]$ $\displaystyle=$
$\displaystyle t_{m}\,\sigma^{2}_{Q},$

resulting in a symmetric correlation structure as follows:

$\displaystyle{\rm Corr}\left[J(t),J(t^{\prime})\right]$ $\displaystyle=$
$\displaystyle\frac{{\rm Cov}[J(t),J(t^{\prime})]}{\sqrt{{\rm Var}[J(t)]{\rm
Var}[J(t^{\prime})]}}$ (91) $\displaystyle=$ $\displaystyle\frac{{\rm
Cov}[J(t_{m}),J(t_{m}+d)]}{\sqrt{{\rm Var}[J(t_{m})]{\rm Var}[J(t_{m}+d)]}}$
$\displaystyle=$
$\displaystyle\frac{t_{m}\,\sigma^{2}_{Q}}{\sqrt{t_{m}\,\sigma^{2}_{Q}\,(t_{m}+d)\,\sigma^{2}_{Q}}}$
$\displaystyle=$ $\displaystyle\sqrt{\frac{t_{m}}{t_{m}+d}}.$

Similarly,

$\displaystyle{\rm Cov}[H(t),H(t^{\prime})]$ $\displaystyle=$
$\displaystyle{\rm Cov}[H(t_{m}),H(t_{m}+d)]$ (92) $\displaystyle=$
$\displaystyle{\rm
Cov}[H(t_{m}),\psi^{d}\,H(t_{m})+\psi^{d-1}\,R_{t_{m}+1}+...+\psi\,R_{t_{m}+d}]$
$\displaystyle=$ $\displaystyle\psi^{d}\,{\rm Var}[H(t_{m})]$ $\displaystyle=$
$\displaystyle\psi^{d}\,\left(\psi^{2(t_{m}-1)}\,\sigma^{2}_{1}+\frac{1-\psi^{2(t_{m}-2)}}{1-\psi^{2}}\,\sigma_{R}^{2}\right),$

since

$\displaystyle{\rm Var}[H(t)]$ $\displaystyle=$ $\displaystyle{\rm
Var}[\psi^{t-1}\,H_{1}+\psi^{t-2}\,R_{2}+...+\psi\,R_{t-1}+R_{t}]$ (93)
$\displaystyle=$
$\displaystyle\psi^{2(t-1)}\,\sigma^{2}_{1}\,+\,\psi^{2(t-2)}\,{\rm
Var}[R_{2}]+...+\psi^{2}\,{\rm Var}[R_{t-1}]+{\rm Var}[R_{t}]$
$\displaystyle=$
$\displaystyle\psi^{2(t-1)}\,\sigma^{2}_{1}+(\psi^{2(t-2)}+\psi^{2(t-3)}+...+\psi^{2}+1)\,\sigma_{R}^{2}$
$\displaystyle=$
$\displaystyle\psi^{2(t-1)}\,\sigma^{2}_{1}+\frac{1-\psi^{2(t-2)}}{1-\psi^{2}}\,\sigma_{R}^{2},$

resulting in a symmetric correlation structure as follows:

$\displaystyle{\rm Corr}\left[H(t),H(t^{\prime})\right]$ $\displaystyle=$
$\displaystyle\frac{{\rm Cov}[H(t),H(t^{\prime})]}{\sqrt{{\rm Var}[H(t)]{\rm
Var}[H(t^{\prime})]}}$ (94) $\displaystyle=$ $\displaystyle\frac{{\rm
Cov}[H(t_{m}),H(t_{m}+d)]}{\sqrt{{\rm Var}[H(t_{m})]{\rm Var}[H(t_{m}+d)]}}$
$\displaystyle=$ $\displaystyle\frac{\psi^{d}\,{\rm Var}[H(t_{m})]}{\sqrt{{\rm
Var}[H(t_{m})]{\rm Var}[H(t_{m}+d)]}}$ $\displaystyle=$
$\displaystyle\psi^{d}\,\sqrt{\frac{{\rm Var}[H(t_{m})]}{{\rm
Var}[H(t_{m}+d)]}}$ $\displaystyle=$
$\displaystyle\psi^{d}\,\sqrt{\frac{\psi^{2(t_{m}-1)}\,\sigma^{2}_{1}+\frac{1-\psi^{2(t_{m}-2)}}{1-\psi^{2}}\,\sigma_{R}^{2}}{\psi^{2(t_{m}+d-1)}\,\sigma^{2}_{1}+\frac{1-\psi^{2(t_{m}+d-2)}}{1-\psi^{2}}\,\sigma_{R}^{2}}}$
$\displaystyle=$
$\displaystyle\psi^{d}\,\sqrt{\frac{(1-\psi^{2})\psi^{2(t_{m}-1)}\,\sigma^{2}_{1}+(1-\psi^{2(t_{m}-2)})\,\sigma_{R}^{2}}{(1-\psi^{2})\psi^{2(t_{m}+d-1)}\,\sigma^{2}_{1}+(1-\psi^{2(t_{m}+d-2)})\,\sigma_{R}^{2}}}.$

Combining the results of Equations (89), (90) and (92) with Equation (88)
leads to a prior belief structure across all $t$ of

$\displaystyle{\rm Cov}[\epsilon(t),\epsilon(t)]$ $\displaystyle=$
$\displaystyle\mathbb{I}_{t=t^{\prime}}\sigma_{A}^{2}+t_{m}\,\sigma^{2}_{Q}+\psi^{d}\,\left(\psi^{2(t_{m}-1)}\,\sigma^{2}_{1}+\frac{1-\psi^{2(t_{m}-2)}}{1-\psi^{2}}\,\sigma_{R}^{2}\right).$
(95)

We now proceed to consider suitable extensions to the uncertain input
scenario. In this case, we need an expectation and covariance structure
spanning ${\rm E}[T],{\rm E}[T^{\prime}],{\rm Cov}[T,T^{\prime}]$-space. We
consider that the statements of expectation in the error structure simply
extend to being ${\rm E}[A(T)]=0$, ${\rm E}[J(T)]=0$ and ${\rm E}[H(T)]=0$,
with logical derivation of these statements being shown using the law of total
expectation. Regarding covariance structure, Equation (89) simply extends to
the following:

${\rm Cov}[A(T),A(T^{\prime})]=\mathbb{I}_{T=T^{\prime}}\sigma_{A}^{2},$ (96)

this being similar to the uncorrelated error term shown in Section A.1 of this
supplementary material.

Equation (90) extends as follows, where the random variables
$T_{m}=\min(T,T^{\prime})$ and $D=|T-T^{\prime}|$ correspond to known
quantities $t_{m}$ and $d$ above:

$\displaystyle{\rm Cov}[J(T),J(T^{\prime})]$ (97) $\displaystyle=$
$\displaystyle{\rm Cov}[J(T_{m}),J(T_{m}+D)]$ $\displaystyle=$
$\displaystyle{\rm E}[{\rm Cov}[J(T_{m}),J(T_{m}+D)|T_{m}=t_{m},D=d]]+{\rm
Cov}[{\rm E}[J(T_{m})|T_{m}],{\rm E}[J(T_{m}+D)|T_{m},D]]$ $\displaystyle=$
$\displaystyle{\rm E}[T_{m}\,\sigma^{2}_{Q}]+0$ $\displaystyle=$
$\displaystyle\sigma^{2}_{Q}\,{\rm E}[T_{m}].$

The tricky specification here is eliciting ${\rm E}[T_{m}]$ from ${\rm
E}[T],{\rm E}[T^{\prime}]$ and ${\rm Cov}[T,T^{\prime}]$, since a function
which takes the minimum of two random quantities is a non-linear function.

If we know which of our random variables is bigger (without loss of
generality, that $T<T^{\prime}$, say), then ${\rm E}[T_{m}]={\rm E}[T]$. This
is not an unrealistic situation; for example, we may have measurements at two
times about which we are uncertain, but know that one of them certainly
happened after the other.

More generally, we propose a reasonable covariance structure as follows. Let
$\mathbb{I}=\mathbb{I}_{T<T^{\prime}}$ and $p={\rm E}[\mathbb{I}]$. We then
condition on this event using the law of total expectation to give:

$\displaystyle{\rm E}[T_{m}]$ $\displaystyle=$ $\displaystyle{\rm E}[{\rm
E}_{\mathbb{I}}[T_{m}]]$ (98) $\displaystyle=$ $\displaystyle{\rm
E}_{\mathbb{I}=1}[T_{m}]\,{\rm E}[\mathbb{I}]+{\rm
E}_{\mathbb{I}=0}[T_{m}]\,(1-{\rm E}[\mathbb{I}])$ $\displaystyle=$
$\displaystyle p\,{\rm E}_{\mathbb{I}=1}[T_{m}]+(1-p)\,{\rm
E}_{\mathbb{I}=0}[T_{m}],$

where

$\displaystyle{\rm E}_{\mathbb{I}=1}[T_{m}]$ $\displaystyle=$
$\displaystyle{\rm E}_{\mathbb{I}=1}[T]$ (99) $\displaystyle=$
$\displaystyle{\rm E}[T]+{\rm Cov}[T,\mathbb{I}]\,{\rm
Var}[\mathbb{I}]^{-1}\,(\mathbb{I}-{\rm E}[\mathbb{I}])$ $\displaystyle=$
$\displaystyle{\rm E}[T]+\frac{{\rm Cov}[T,\mathbb{I}]}{p(1-p)}\times(1-p)\,$
$\displaystyle=$ $\displaystyle{\rm E}[T]+\frac{{\rm Cov}[T,\mathbb{I}]}{p},$
$\displaystyle{\rm E}_{\mathbb{I}=0}[T_{m}]$ $\displaystyle=$
$\displaystyle{\rm E}_{\mathbb{I}=0}[T^{\prime}]$ (100) $\displaystyle=$
$\displaystyle{\rm E}[T^{\prime}]+{\rm Cov}[T^{\prime},\mathbb{I}]\,{\rm
Var}[\mathbb{I}]^{-1}\,(\mathbb{I}-{\rm E}[\mathbb{I}])$ $\displaystyle=$
$\displaystyle{\rm E}[T^{\prime}]+\frac{{\rm
Cov}[T^{\prime},\mathbb{I}]}{p(1-p)}\times(-p)\,$ $\displaystyle=$
$\displaystyle{\rm E}[T^{\prime}]-\frac{{\rm
Cov}[T^{\prime},\mathbb{I}]}{1-p},$

and

$\displaystyle{\rm Cov}[T,\mathbb{I}]$ $\displaystyle=$ $\displaystyle{\rm
Corr}\left[T,\mathbb{I}\right]\,\sqrt{{\rm Var}[T]\,{\rm Var}[\mathbb{I}]}$
(101) $\displaystyle=$ $\displaystyle{\rm
Corr}\left[T,\mathbb{I}\right]\,\sqrt{{\rm Var}[T]\,p(1-p)},$
$\displaystyle{\rm Cov}[T^{\prime},\mathbb{I}]$ $\displaystyle=$
$\displaystyle{\rm Corr}\left[T^{\prime},\mathbb{I}\right]\,\sqrt{{\rm
Var}[T^{\prime}]\,{\rm Var}[\mathbb{I}]}$ (102) $\displaystyle=$
$\displaystyle{\rm Corr}\left[T^{\prime},\mathbb{I}\right]\,\sqrt{{\rm
Var}[T^{\prime}]\,p(1-p)},$

so that:

$\displaystyle{\rm E}[T_{m}]$ $\displaystyle=$ $\displaystyle p\left({\rm
E}[T]+\frac{1}{p}\,{\rm Corr}\left[T,\mathbb{I}\right]\,\sqrt{{\rm
Var}[T]\,p(1-p)}\right)$ (103) $\displaystyle\,\,\,\,\,\,+\,(1-p)\left({\rm
E}[T^{\prime}]-\frac{1}{1-p}\,{\rm
Corr}\left[T^{\prime},\mathbb{I}\right]\,\sqrt{{\rm
Var}[T^{\prime}]\,p(1-p)}\right)$ $\displaystyle=$ $\displaystyle p\,{\rm
E}[T]+(1-p)\,{\rm E}[T^{\prime}]$
$\displaystyle\,\,\,\,\,\,+\sqrt{p(1-p)}\left({\rm
Corr}\left[T,\mathbb{I}\right]\,\sqrt{{\rm Var}[T]}-{\rm
Corr}\left[T^{\prime},\mathbb{I}\right]\,\sqrt{{\rm Var}[T^{\prime}]}\right).$

We may be able to specify values for all quantities required in Equation
(103), in which case we have our expression for ${\rm E}[T_{m}]$. If
specification of ${\rm Corr}\left[T,\mathbb{I}\right]$ and ${\rm
Corr}\left[T^{\prime},\mathbb{I}\right]$ is proving difficult, then let us
first note that it is logical to assume that $-1<{\rm
Corr}\left[T,\mathbb{I}\right]<0$ and $0<{\rm
Corr}\left[T^{\prime},\mathbb{I}\right]<1$, so that

${\rm E}[T_{m}]>p\,{\rm E}[T]+(1-p)\,{\rm
E}[T^{\prime}]-\sqrt{p(1-p)}\left(\sqrt{{\rm Var}[T]}+\sqrt{{\rm
Var}[T^{\prime}]}\right).$ (104)

Again, if we feel able to specify a value for $p$, then we can use it in
Equation (104). Alternatively, in order to be conservative with respect to
variance resolution, we can minimise the exspression on the right had side of
Equation (104) over $p\in[0,1]$.

First, we note that both Equations (103) and (104) are of the form

$h(p)=ap+b(1-p)-c\sqrt{p(1-p)},$ (105)

with $a={\rm E}[T]$, $b={\rm E}[T^{\prime}]$, and $c=\sqrt{{\rm
Var}[T]}+\sqrt{{\rm Var}[T^{\prime}]}$ or $c={\rm
Corr}\left[T,\mathbb{I}\right]\,\sqrt{{\rm Var}[T]}-{\rm
Corr}\left[T^{\prime},\mathbb{I}\right]\,\sqrt{{\rm Var}[T^{\prime}]}$.

Use elementary calculus to find the minimum, by first taking the derivative of
$h(p)$ with respect to $p$:

$\frac{\mathrm{d}h(p)}{\mathrm{d}p}=a-b+\frac{c(2p-1)}{2\sqrt{p(1-p)}}.$ (106)

We then set $\frac{\mathrm{d}h(p)}{\mathrm{d}p}=0$ to find $p$ which minimises
$h(p)$:

$\displaystyle a-b+\frac{c(2p-1)}{2\sqrt{p(1-p)}}=0$ (107)
$\displaystyle\implies$ $\displaystyle c^{2}(2p-1)^{2}=4p(1-p)(b-a)^{2}$
$\displaystyle\implies$ $\displaystyle
p^{2}(4c^{2}+4(b-a)^{2})-p(4c^{2}+4(b-a)^{2})+c^{2}=0,$

which is of the form

$kp^{2}-kp+c^{2}=0,$ (108)

with $k=4(c^{2}+(b-a)^{2})$. Equation (108) has the solution

$\displaystyle\hat{p}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(1+\sqrt{\frac{m^{2}}{n}}\right),$ (109)

where $m=b-a$ and $n=c^{2}+m^{2}$.

Inserting this into Expression (105), we have

$\displaystyle h(\hat{p})$ $\displaystyle=$ $\displaystyle
a\hat{p}+b(1-\hat{p})-c\sqrt{\hat{p}(1-\hat{p})}.$ (110) $\displaystyle=$
$\displaystyle\frac{1}{2}\left(a+b-\sqrt{n}\left(\frac{c^{2}}{n}\pm\frac{m^{2}}{n}\right)\right)$

The minimum of Expression (110) occurs when $\pm$ is $+$, hence

$\displaystyle h(\hat{p})$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(a+b-\sqrt{n}\left(\frac{n-m^{2}}{n}+\frac{m^{2}}{n}\right)\right)$
(111) $\displaystyle=$
$\displaystyle\frac{1}{2}\left(a+b-\sqrt{(b-a)^{2}+c^{2}}\right).$

Therefore we have that

${\rm E}[T_{m}]\geq\frac{1}{2}\left({\rm E}[T]+{\rm
E}[T^{\prime}]-\sqrt{\left({\rm E}[T^{\prime}]-{\rm
E}[T]\right)^{2}+\left(\sqrt{{\rm Var}[T]}+\sqrt{{\rm
Var}[T^{\prime}]}\right)^{2}}\right),$ (112)

and we might set

${\rm Cov}[J(T),J(T^{\prime})]=\frac{\sigma^{2}_{Q}}{2}\left({\rm E}[T]+{\rm
E}[T^{\prime}]-\sqrt{\left({\rm E}[T^{\prime}]-{\rm
E}[T]\right)^{2}+\left(\sqrt{{\rm Var}[T]}+\sqrt{{\rm
Var}[T^{\prime}]}\right)^{2}}\right).$ (113)

From Equation (90) with $t=t^{\prime}$, it follows that

$\displaystyle{\rm Var}[J(T)]$ $\displaystyle=$ $\displaystyle{\rm E}[{\rm
Var}[J(T)|T=t]]+{\rm Var}[{\rm E}[J(T)|T]]$ (114) $\displaystyle=$
$\displaystyle{\rm E}[\sigma_{Q}^{2}\,T]+0$ $\displaystyle=$
$\displaystyle\sigma^{2}_{Q}\,{\rm E}[T],$

since $T_{m}=T$. We thus make use of an indicator function in the expression
for ${\rm Cov}[J(T),J(T^{\prime})]$ to ensure it is appropriate when it is
known that $T=T^{\prime}$:

$\displaystyle{\rm Cov}[J(T),J(T^{\prime})]$ $\displaystyle=$
$\displaystyle\frac{\sigma^{2}_{Q}}{2}\left({\rm E}[T]+{\rm
E}[T^{\prime}]-\sqrt{\left({\rm E}[T^{\prime}]-{\rm
E}[T]\right)^{2}+\left(1-\mathbb{I}_{T=T^{\prime}}\right)\left(\sqrt{{\rm
Var}[T]}+\sqrt{{\rm Var}[T^{\prime}]}\right)^{2}}\right).$

Note that:

  1. 1.

When ${\rm Var}[T]={\rm Var}[T^{\prime}]=0$, the expression on the right hand
side of Inequality (112) yields

$\frac{1}{2}\left(t+t^{\prime}-\sqrt{(t-t^{\prime})^{2}}\right)=\frac{1}{2}\left(t+t^{\prime}-|t-t^{\prime}|\right)=\min(t,t^{\prime}),$

thus resolving to the known input scenario.

  2. 2.

If it is known that $T<T^{\prime}$, then, as mentioned earlier, it would be
logical to have that ${\rm E}[T_{m}]={\rm E}[T]$. Since

${\rm E}[T]=\frac{1}{2}\left({\rm E}[T]+{\rm E}[T^{\prime}]-\sqrt{\left({\rm
E}[T^{\prime}]-{\rm E}[T]\right)^{2}}\right),$ (116)

it is clear that ${\rm E}[T_{m}]$ satisfies Inequality (112) in this case.

We extend Equation (92) by noticing that

$\displaystyle{\rm Cov}[H(T),H(T^{\prime})]$ $\displaystyle=$
$\displaystyle{\rm Cov}[H(T_{m}),H(T_{m}+D)]$ (117) $\displaystyle=$
$\displaystyle{\rm E}[{\rm Cov}[H(T_{m}),H(T_{m}+D)|T_{m}=t_{m},D=d]]$
$\displaystyle\,\,\,\,\,\,+{\rm Cov}[{\rm E}[H(T_{m})|T_{m}=t_{m}],{\rm
E}[H(T_{m}+D)|T_{m}=t_{m},D=d]]$ $\displaystyle=$ $\displaystyle{\rm
E}[\psi^{D}\,{\rm Var}[H(T_{m})]]+0$ $\displaystyle=$ $\displaystyle{\rm
E}[\psi^{D}\,\left(\psi^{2(T_{m}-1)}\,\sigma^{2}_{1}+\frac{1-\psi^{2(T_{m}-2)}}{1-\psi^{2}}\,\sigma_{R}^{2}\right)]$
$\displaystyle\geq$ $\displaystyle{\rm
E}[\psi^{D}\,\left(\psi^{2(T_{m}-1)}\,\sigma^{2}_{1}+\sigma_{R}^{2}\right)]$
$\displaystyle=$ $\displaystyle\sigma^{2}_{1}\,{\rm
E}[\psi^{D+2T_{m}-2}]+\sigma^{2}_{R}\,{\rm E}[\psi^{D}]$ $\displaystyle=$
$\displaystyle\sigma^{2}_{1}\,\exp(\log({\rm
E}[\psi^{D+2T_{m}-2}]))+\sigma^{2}_{R}\,\exp(\log({\rm E}[\psi^{D}]))$
$\displaystyle\geq$ $\displaystyle\sigma^{2}_{1}\,\exp({\rm
E}[\log(\psi^{D+2T_{m}-2})])+\sigma^{2}_{R}\,\exp({\rm E}[\log(\psi^{D}]))$
$\displaystyle=$ $\displaystyle\sigma_{1}^{2}\exp(\log\psi\,{\rm
E}[D+2T_{m}-2])+\sigma_{R}^{2}\exp(\log\psi\,{\rm E}[D])$ $\displaystyle=$
$\displaystyle\sigma_{1}^{2}\exp(\log\psi\,{\rm
E}[T_{M}+T_{m}-2])+\sigma_{R}^{2}\exp(\log\psi\,{\rm E}[T_{M}-T_{m}]).$

where the first $\geq$ is a result of the fact that, if $T_{m}\geq 3$, then
$(1-\psi^{2(T_{m}-2)})/(1-\psi^{2})\geq 1$, the second $\geq$ is a result of
the fact that $\log$ is a concave function, i.e. that ${\rm E}[f(X)]\leq
f({\rm E}[X])$ if $f(X)=\log(X)$., and the final line makes use of the fact
that $d=t_{M}-t_{m}$, where $t_{M}=\max(t,t^{\prime})$, and correspondingly
$D=T_{M}-T_{m}$, where $T_{M}=\max(T,T^{\prime})$.

We note that ${\rm E}[T_{M}+T_{m}]={\rm E}[T]+{\rm E}[T^{\prime}]$, and thus
seek an upper bound for ${\rm E}[T_{M}-T_{m}]$ in order to find a lower bound
for Expression (117) (since $\log\psi$ is negative if $\psi\in(0,1)$).

Similar to Equation (98), we use the law of total expectation on ${\rm
E}[T_{M}]$ to get that

$\displaystyle{\rm E}[T_{M}]$ $\displaystyle=$ $\displaystyle{\rm E}[{\rm
E}_{\mathbb{I}}[T_{M}]]$ (118) $\displaystyle=$ $\displaystyle{\rm
E}_{\mathbb{I}=1}[T_{M}]\,{\rm E}[\mathbb{I}]+{\rm
E}_{\mathbb{I}=0}[T_{M}]\,(1-{\rm E}[\mathbb{I}])$ $\displaystyle=$
$\displaystyle p\,{\rm E}_{\mathbb{I}=1}[T_{M}]+(1-p)\,{\rm
E}_{\mathbb{I}=0}[T_{M}],$

where:

$\displaystyle{\rm E}_{\mathbb{I}=1}[T_{M}]$ $\displaystyle=$
$\displaystyle{\rm E}_{\mathbb{I}=1}[T^{\prime}]$ (119) $\displaystyle=$
$\displaystyle{\rm E}[T^{\prime}]+{\rm Cov}[T^{\prime},\mathbb{I}]\,{\rm
Var}[\mathbb{I}]^{-1}\,(\mathbb{I}-{\rm E}[\mathbb{I}])$ $\displaystyle=$
$\displaystyle{\rm E}[T^{\prime}]+\frac{{\rm
Cov}[T^{\prime},\mathbb{I}]}{p(1-p)}\times(1-p)\,$ $\displaystyle=$
$\displaystyle{\rm E}[T^{\prime}]+\frac{{\rm Cov}[T^{\prime},\mathbb{I}]}{p},$
$\displaystyle{\rm E}_{\mathbb{I}=0}[T_{M}]$ $\displaystyle=$
$\displaystyle{\rm E}_{\mathbb{I}=0}[T]$ (120) $\displaystyle=$
$\displaystyle{\rm E}[T]+{\rm Cov}[T,\mathbb{I}]\,{\rm
Var}[\mathbb{I}]^{-1}\,(\mathbb{I}-{\rm E}[\mathbb{I}])$ $\displaystyle=$
$\displaystyle{\rm E}[T]+\frac{{\rm Cov}[T,\mathbb{I}]}{p(1-p)}\times(-p)\,$
$\displaystyle=$ $\displaystyle{\rm E}[T]-\frac{{\rm
Cov}[T,\mathbb{I}]}{1-p}.$

Combining Equations (98) and (118), we get that

$\displaystyle{\rm E}[T_{M}-T_{m}]$ $\displaystyle=$ $\displaystyle{\rm
E}[T_{M}]-{\rm E}[T_{m}]$ (121) $\displaystyle=$ $\displaystyle p\,{\rm
E}[T^{\prime}]+{\rm Cov}[T^{\prime},\mathbb{I}]+(1-p)\,{\rm E}[T]-{\rm
Cov}[T,\mathbb{I}]$ $\displaystyle\,\,\,\,\,\,-\left(p\,{\rm E}[T]+{\rm
Cov}[T,\mathbb{I}]+(1-p)\,{\rm E}[T^{\prime}]-{\rm
Cov}[T^{\prime},\mathbb{I}]\right)$ $\displaystyle=$
$\displaystyle(1-2p)\left({\rm E}[T]-{\rm E}[T^{\prime}]\right)+2\left({\rm
Cov}[T^{\prime},\mathbb{I}]-{\rm Cov}[T,\mathbb{I}]\right)$ $\displaystyle=$
$\displaystyle(1-2p)\left({\rm E}[T]-{\rm E}[T^{\prime}]\right)$
$\displaystyle\,\,\,\,\,\,+2\left({\rm
Corr}\left[T^{\prime},\mathbb{I}\right]\sqrt{{\rm Var}[T^{\prime}]p(1-p)}-{\rm
Corr}\left[T,\mathbb{I}\right]\sqrt{{\rm Var}[T]p(1-p)}\right)$
$\displaystyle=$ $\displaystyle(1-2p)\left({\rm E}[T]-{\rm
E}[T^{\prime}]\right)$ $\displaystyle\,\,\,\,\,\,+2\sqrt{p(1-p)}\left({\rm
Corr}\left[T^{\prime},\mathbb{I}\right]\sqrt{{\rm Var}[T^{\prime}]}-{\rm
Corr}\left[T,\mathbb{I}\right]\sqrt{{\rm Var}[T]}\right).$

Again, if we don’t feel able to specify ${\rm
Corr}\left[T^{\prime},\mathbb{I}\right]$ or ${\rm
Corr}\left[T,\mathbb{I}\right]$, then using the assumptions that $-1<{\rm
Corr}\left[T,\mathbb{I}\right]<0$ and $0<{\rm
Corr}\left[T^{\prime},\mathbb{I}\right]<1$, we have that

${\rm E}[T_{M}-T_{m}]<(1-2p)\left({\rm E}[T]-{\rm
E}[T^{\prime}]\right)+2\sqrt{p(1-p)}\left(\sqrt{{\rm
Var}[T^{\prime}]}+\sqrt{{\rm Var}[T]}\right).$ (122)

If we have specified a value for $p$, again we can use this in Expression
(122), otherwise we can maximise over $p\in[0,1]$. This follows from the
calculations to bound ${\rm E}[T_{m}]$ above, since we note again that we have
an expression of the form

$h(p)=ap+b(1-p)-c\sqrt{p(1-p)},$ (123)

this time with $a=-\left({\rm E}[T]-{\rm E}[T^{\prime}]\right)={\rm
E}[T^{\prime}]-{\rm E}[T]$, $b={\rm E}[T]-{\rm E}[T^{\prime}]$, and
$c=-2(\sqrt{{\rm Var}[T]}+\sqrt{{\rm Var}[T^{\prime}]})$, which is maximised
by Equation (110), this time with $-$ instead of $+$, hence

$\displaystyle h(\hat{p})$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(a+b-\sqrt{n}\left(\frac{c^{2}}{n}-\frac{m^{2}}{n}\right)\right)$
(124) $\displaystyle=$
$\displaystyle\frac{1}{2}\left(a+b+\sqrt{(b-a)^{2}+c^{2}}\right).$

Therefore we have that

$\displaystyle{\rm E}[T_{M}-T_{m}]$ $\displaystyle\leq$
$\displaystyle\frac{1}{2}\bigg{(}{\rm E}[T^{\prime}]-{\rm E}[T]+{\rm
E}[T]-{\rm E}[T^{\prime}]$ (125) $\displaystyle\,\,\,\,\,\,+\,\sqrt{\left({\rm
E}[T]-{\rm E}[T^{\prime}]+{\rm E}[T]-{\rm
E}[T^{\prime}]\right)^{2}+4\left(\sqrt{{\rm Var}[T]}+\sqrt{{\rm
Var}[T^{\prime}]}\right)^{2}}\bigg{)}$ $\displaystyle=$
$\displaystyle\sqrt{\left({\rm E}[T]-{\rm
E}[T^{\prime}]\right)^{2}+\left(\sqrt{{\rm Var}[T]}+\sqrt{{\rm
Var}[T^{\prime}]}\right)^{2}}.$

For $T\neq T^{\prime}$, we therefore propose

$\displaystyle{\rm Cov}[H(T),H(T^{\prime})]$ (126) $\displaystyle=$
$\displaystyle\sigma_{1}^{2}\exp(\log\psi\,({\rm E}[T]+{\rm
E}[T^{\prime}]-2))$
$\displaystyle\,\,\,\,\,\,+\,\sigma_{R}^{2}\exp(\log\psi\,\sqrt{\left({\rm
E}[T]-{\rm E}[T^{\prime}]\right)^{2}+\left(\sqrt{{\rm Var}[T]}+\sqrt{{\rm
Var}[T^{\prime}]}\right)^{2}}).$

When $T=T^{\prime}$, the above expression needs to be commensurate with a
reasonable expression for ${\rm Var}[H(T)]$. Unlike ${\rm Var}[J(T)]$, we find
that ${\rm Var}[H(T)]$ itself needs approximating. However, for a variance
term we need to find an upper bound in order to be conservative with regard to
prior variance specification (i.e. overestimate rather than underestimate).

We have that

$\displaystyle{\rm Var}[H(T)]$ $\displaystyle=$ $\displaystyle{\rm E}[{\rm
Var}[H(T)|T=t]]+{\rm Var}[{\rm E}[H(T)|T]]$ (127) $\displaystyle=$
$\displaystyle{\rm
E}[\psi^{2(T-1)}\,\sigma^{2}_{1}+\frac{1-\psi^{2(T-2)}}{1-\psi^{2}}\,\sigma_{R}^{2}]$
$\displaystyle=$ $\displaystyle{\rm
E}[\psi^{2(T-1)}\,(\sigma_{1}^{2}-\sigma_{R}^{2})+\frac{1-\psi^{2(T-1)}}{1-\psi^{2}}\,\sigma_{R}^{2}]$
$\displaystyle\leq$ $\displaystyle\mathbb{I}_{\sigma_{1}^{2}\leq
K}K+(1-\mathbb{I}_{\sigma_{1}^{2}\leq
K})(\psi^{4}\sigma_{1}^{2}+\sigma_{R}^{2}),$

where $K=\frac{\sigma_{R}^{2}}{1-\psi^{2}}$. The inequality holds because as
$t\rightarrow\infty$, the expression tends to $K$, but this is from above or
below depending on whether the variance gets larger or smaller for larger
values of $t$. We can assess this by considering the difference between the
variance at time $t$ and $t+1$.

$\displaystyle{\rm Var}[H(t+1)]-{\rm Var}[H(t)]$ $\displaystyle=$
$\displaystyle\psi^{2(t-1)}\sigma_{R}^{2}-(\psi^{2(t-1)}-\psi^{2t})\sigma_{1}^{2}$
(128) $\displaystyle=$
$\displaystyle\psi^{2(t-1)}\sigma_{R}^{2}+(\psi^{2t}-\psi^{2(t-1)})\sigma_{1}^{2}$
$\displaystyle=$
$\displaystyle\psi^{2(t-1)}(\sigma_{R}^{2}+\psi^{2}\sigma_{1}^{2}-\sigma_{1}^{2})$

thus if $\sigma_{R}^{2}+\psi^{2}\sigma_{1}^{2}-\sigma_{1}^{2}\geq 0$, or in
other words $\sigma_{1}^{2}\leq K$, then the difference is positive, thus we
approach the limit at infinity from below, and the expression is maximised at
this limit. On the other hand, if
$\sigma_{R}^{2}+\psi^{2}\sigma_{1}^{2}-\sigma_{1}^{2}\leq 0$, or in other
words $\sigma_{1}^{2}\geq K$, then the difference is negative, and we approach
the limit at infinity from above, in which case the expression is maximised
for $t=3$ (since we enforced $t\geq 3$), which yields
$\psi^{4}\sigma_{1}^{2}+\sigma_{R}^{2}$.

We can now set

$\displaystyle{\rm Cov}[H(T),H(T^{\prime})]$ (129) $\displaystyle=$
$\displaystyle\mathbb{I}_{T=T^{\prime}}L+(1-\mathbb{I}_{T=T^{\prime}})\bigg{(}\sigma_{1}^{2}\exp(\log\psi\,({\rm
E}[T]+{\rm E}[T^{\prime}]-2))$
$\displaystyle\,\,\,\,\,\,+\,\sigma_{R}^{2}\exp\left(\log\psi\,\sqrt{\left({\rm
E}[T]-{\rm E}[T^{\prime}]\right)^{2}+\left(\sqrt{{\rm Var}[T]}+\sqrt{{\rm
Var}[T^{\prime}]}\right)^{2}}\right)\bigg{)}$

where

$L=\mathbb{I}_{\sigma_{1}^{2}\leq K}K+(1-\mathbb{I}_{\sigma_{1}^{2}\leq
K})(\psi^{4}\sigma_{1}^{2}+\sigma_{R}^{2})$

Finally, to conclude our prior specification in the case of uncertain inputs
for the model presented in this section, note that, with
$\mathbf{T}=(1,T)^{T}$, we have

$\displaystyle{\rm E}[\mathbf{T}^{T}\bm{\mathit{\Delta}}\mathbf{T}^{\prime}]$
(136) $\displaystyle=$ $\displaystyle\textrm{trace}(\bm{\mathit{\Delta}}({\rm
E}[\mathbf{T}^{\prime}]{\rm E}[\mathbf{T}^{T}]+{\rm
Cov}[\mathbf{T},\mathbf{T}^{\prime}]))$ $\displaystyle=$
$\displaystyle\textrm{trace}\left(\left(\begin{array}[]{cc}\delta_{00}&\delta_{01}\\\
\delta_{10}&\delta_{11}\end{array}\right)\left(\begin{array}[]{cc}1&{\rm
E}[T]\\\ {\rm E}[T^{\prime}]&{\rm E}[T]{\rm
E}[T^{\prime}]\end{array}\right)+\left(\begin{array}[]{cc}0&0\\\ 0&{\rm
Cov}[T,T^{\prime}]\end{array}\right)\right)$ $\displaystyle=$
$\displaystyle\textrm{trace}\left(\begin{array}[]{cc}\delta_{00}+\delta_{01}{\rm
E}[T^{\prime}]&\delta_{00}{\rm E}[T]+\delta_{01}{\rm E}[T]{\rm
E}[T^{\prime}]\\\ \delta_{10}+\delta_{11}{\rm E}[T^{\prime}]&\delta_{10}{\rm
E}[T]+\delta_{11}{\rm E}[T]{\rm E}[T^{\prime}]+{\rm
Cov}[T,T^{\prime}]\end{array}\right)$ (139) $\displaystyle=$
$\displaystyle\delta_{00}+\delta_{01}({\rm E}[T^{\prime}]+{\rm
E}[T])+\delta_{11}{\rm E}[T]{\rm E}[T^{\prime}]+{\rm Cov}[T,T^{\prime}],$
(140)

and

$\bm{\Gamma}^{T}{\rm
Cov}[\mathbf{T},\mathbf{T}^{\prime}]\bm{\Gamma}=\left(\begin{array}[]{cc}\gamma_{0}&\gamma_{1}\end{array}\right)\left(\begin{array}[]{cc}0&0\\\
0&{\rm
Cov}[T,T^{\prime}]\end{array}\right)\left(\begin{array}[]{c}\gamma_{0}\\\
\gamma_{1}\end{array}\right)=\gamma_{1}^{2}{\rm Cov}[T,T^{\prime}],$ (141)

so that full prior specification for $y(T)$ can therefore be given by

${\rm E}[y(T)]=\gamma_{0}+{\rm E}[T]\gamma_{1},$ (142)

and

$\displaystyle{\rm Cov}[y(T),y(T^{\prime})]$ (143) $\displaystyle=$
$\displaystyle\delta_{00}+\delta_{01}({\rm E}[T^{\prime}]+{\rm
E}[T])+\delta_{11}{\rm E}[T]{\rm E}[T^{\prime}]+{\rm Cov}[T,T^{\prime}]$
$\displaystyle\,\,\,\,\,\,+\,\,\gamma_{1}^{2}{\rm Cov}[T,T^{\prime}]$
$\displaystyle\,\,\,\,\,\,+\,\,\mathbb{I}_{T=T^{\prime}}\sigma_{A}^{2}$
$\displaystyle\,\,\,\,\,\,+\,\,\frac{\sigma^{2}_{Q}}{2}\left({\rm E}[T]+{\rm
E}[T^{\prime}]-\sqrt{\left({\rm E}[T^{\prime}]-{\rm
E}[T]\right)^{2}+\left(1-\mathbb{I}_{T=T^{\prime}}\right)\left(\sqrt{{\rm
Var}[T]}+\sqrt{{\rm Var}[T^{\prime}]}\right)^{2}}\right)$
$\displaystyle\,\,\,\,\,\,+\,\,\mathbb{I}_{T=T^{\prime}}L+(1-\mathbb{I}_{T=T^{\prime}})\bigg{(}\sigma_{1}^{2}\exp(\log\psi\,({\rm
E}[T]+{\rm E}[T^{\prime}]-2))$
$\displaystyle\,\,\,\,\,\,\qquad+\sigma_{R}^{2}\exp\left(\log\psi\,\sqrt{\left({\rm
E}[T]-{\rm E}[T^{\prime}]\right)^{2}+\left(\sqrt{{\rm Var}[T]}+\sqrt{{\rm
Var}[T^{\prime}]}\right)^{2}}\right)\bigg{)}.$

##  Appendix B Proof of Lemmas 1 - 4

In this section, we prove Lemmas 1 - 4 of Sections 4 and 5 of the main text.

### Proof of Lemma 1

Rewrite Equation (60) as

$\begin{split}c(\mathbf{X},\mathbf{X}^{\prime})&=\exp\left\\{-{\rm
E}[(\mathbf{X}-\mathbf{X}^{\prime})^{T}\Theta^{-2}(\mathbf{X}-\mathbf{X}^{\prime})]\right\\}\\\
&=\exp\left\\{-{\rm
E}[\mathbf{X}^{T}\Theta^{-2}\mathbf{X}]\right\\}\exp\left\\{-{\rm
E}[(\mathbf{X}^{\prime})^{T}\Theta^{-2}\mathbf{X}^{\prime}]\right\\}\exp\left\\{2{\rm
E}[\mathbf{X}^{T}\Theta^{-2}\mathbf{X}^{\prime}]\right\\}\,,\\\ \end{split}$

with kernel

$k(\mathbf{X},\mathbf{X}^{\prime})=\exp\left\\{2{\rm
E}[\mathbf{X}^{T}\Theta^{-2}\mathbf{X}^{\prime}]\right\\}\,.$

The kernel $k(\mathbf{X},\mathbf{X}^{\prime})$ is positive definite since

$\mathrm{tr}(\Theta^{-2}{\rm Var}[\mathbf{X}])+{\rm
E}[\mathbf{X}]^{T}\Theta^{-2}{\rm E}[\mathbf{X}]\geq 0$

with equality if and only if $\mathbf{X}$ has a degenerate distribution at
zero. Postive definiteness of the function $c(\mathbf{X},\mathbf{X}^{\prime})$
then follows from standard properties of kernels (see Shawe-Taylor &
Cristianini 2011, ch. 3)

### Proof of Lemma 2

Covariance can be derived from conditional quantities using the law of total
covariance, hence:

${\rm Cov}[\mathbf{u}(\mathbf{X}),\mathbf{u}(\mathbf{X}^{\prime})]={\rm
E}[{\rm
Cov}[\mathbf{u}(\mathbf{X}),\mathbf{u}(\mathbf{X}^{\prime})\,|\,\mathbf{X}=\mathbf{x},\mathbf{X}^{\prime}=\mathbf{x}^{\prime}]]\\\
+{\rm Cov}[{\rm E}[\mathbf{u}(\mathbf{X})\,|\,\mathbf{X}],{\rm
E}[\mathbf{u}(\mathbf{X}^{\prime})\,|\,\mathbf{X}^{\prime}]]\,.$

Under the assumption that ${\rm E}[\mathbf{u}(\mathbf{X})]=0$, it follows that

${\rm Cov}[{\rm E}[\mathbf{u}(\mathbf{X})\,|\,\mathbf{X}],{\rm
E}[\mathbf{u}(\mathbf{X}^{\prime})\,|\,\mathbf{X}^{\prime}]]=0$

. Hence for conditional covariance (61), we have

$\begin{split}{\rm
Cov}[\mathbf{u}(\mathbf{X}),\mathbf{u}(\mathbf{X}^{\prime})]&=\mathrm{E}\left[\exp\left(-\sum_{r=1}^{p}\left\\{\frac{X_{(r)}-X^{\prime}_{(r)}}{\theta_{r}}\right\\}^{2}\right)\right]\,\bm{\mathit{\Sigma}}\\\
&\succeq\exp\left(-\sum_{r=1}^{p}\mathrm{E}\left[\left\\{\frac{X_{(r)}-X^{\prime}_{(r)}}{\theta_{r}}\right\\}^{2}\,\right]\right)\,\bm{\mathit{\Sigma}}\,,\\\
\end{split}$

under the Loewner (partial) ordering, following from an application of
Jensen’s inequality and the positive-definiteness of $\bm{\mathit{\Sigma}}$.
As the diagonal entries of $\bm{\mathit{\Sigma}}$ are all non-negative,
elementwise inequality for the $kk$th entry follows directly ($k=1,\ldots,q$).

### Proof of Lemma 3

We begin by expanding the terms of the Bayes linear update as follows:

$\displaystyle{\rm E}_{\mathbf{F}}[\mathbf{f}(\mathbf{X})]$ $\displaystyle=$
$\displaystyle{\rm E}_{\mathbf{F}}[\bm{w}(\mathbf{X})\,\bm{\beta}]+{\rm
E}_{\mathbf{F}}[\mathbf{u}(\mathbf{X})]\,.$

Taking the two parts of the right-hand side of this equation separately, we
first have that

${\rm E}_{\mathbf{F}}[\bm{w}(\mathbf{X})\,\bm{\beta}]={\rm
E}_{\mathbf{F}}[\bm{w}(\mathbf{X})]\,{\rm E}_{\mathbf{F}}[\bm{\beta}]={\rm
E}[\bm{w}(\mathbf{X})]\,{\rm E}_{\mathbf{F}}[\bm{\beta}]\,,$ (144)

where we have used the facts that ${\rm
E}_{\mathbf{F}}[\bm{w}(\mathbf{X})]={\rm E}[\bm{w}(\mathbf{X})]$ since ${\rm
Cov}[\mathbf{X},\mathbf{F}]=\mathbf{0}$, and ${\rm
Cov}_{\mathbf{F}}\left[\bm{w}(\mathbf{X}),\bm{\beta}\right]=\mathbf{0}$ since
${\rm
Cov}_{\mathbf{F}}\left[\mathbf{g}(\mathbf{X}),\bm{\beta}\right]=\mathbf{0}$.
We then have, using basic rules of linear algebra (Mardia et al. 1979), that

$\displaystyle{\rm E}_{\mathbf{F}}[\mathbf{u}(\mathbf{X})]$ $\displaystyle=$
$\displaystyle{\rm E}[\mathbf{u}(\mathbf{X})]+{\rm
Cov}[\mathbf{u}(\mathbf{X}),\mathbf{F}]\,{\rm
Var}[\mathbf{F}]^{-1}\,(\mathbf{F}-{\rm E}[\mathbf{F}])$ (145)
$\displaystyle=$ $\displaystyle{\rm
Cov}[\mathbf{u}(\mathbf{X}),\mathbf{U}]\,{\rm
Var}[\mathbf{F}]^{-1}\,(\mathbf{F}-\bm{W}\,\bm{\Gamma})$ $\displaystyle=$
$\displaystyle\bm{v}(\mathbf{X})\,(\bm{W}\,\bm{\mathit{\Delta}}\,\bm{W}^{T}+\bm{\mathit{\Omega}})^{-1}\,(\mathbf{F}-\bm{W}\,\bm{\Gamma})$
$\displaystyle=$
$\displaystyle\bm{v}(\mathbf{X})\,\big{(}(\bm{W}\,\bm{\mathit{\Delta}}\,\bm{W}^{T}+\bm{\mathit{\Omega}})^{-1}\,\mathbf{F}-(\bm{W}\,\bm{\mathit{\Delta}}\,\bm{W}^{T}+\bm{\mathit{\Omega}})^{-1}\,\bm{W}\,\bm{\Gamma}\big{)}$
$\displaystyle=$
$\displaystyle\bm{v}(\mathbf{X})\bigg{(}\big{(}\bm{\mathit{\Omega}}^{-1}-\bm{\mathit{\Omega}}^{-1}\,\bm{W}\,(\bm{\mathit{\Delta}}^{-1}+\bm{W}^{T}\,\bm{\mathit{\Omega}}^{-1}\,\bm{W})^{-1}\,\bm{W}^{T}\,\bm{\mathit{\Omega}}^{-1}\big{)}\,\mathbf{F}$
$\displaystyle\,\,\,\,\,\,-\bm{\mathit{\Omega}}^{-1}\,\bm{W}\,(\bm{W}^{T}\,\bm{\mathit{\Omega}}^{-1}\,\bm{W}+\bm{\mathit{\Delta}}^{-1})^{-1}\,\bm{\mathit{\Delta}}^{-1}\,\bm{W}\bigg{)}$
$\displaystyle=$
$\displaystyle\bm{v}(\mathbf{X})\,\bm{\mathit{\Omega}}^{-1}\,\big{(}\mathbf{F}-\bm{W}\,(\bm{\mathit{\Delta}}^{-1}+\bm{W}^{T}\,\bm{\mathit{\Omega}}^{-1}\,\bm{W})^{-1}\,(\bm{\mathit{\Delta}}^{-1}\,\bm{W}+\bm{W}\,\bm{\mathit{\Omega}}^{-1}\,\mathbf{F})\big{)}$
$\displaystyle=$
$\displaystyle\bm{v}(\mathbf{X})\,\bm{\mathit{\Omega}}^{-1}\,(\mathbf{F}-\bm{W}\,{\rm
E}_{\mathbf{F}}[\bm{\beta}])\,,$

where we have that

$\displaystyle{\rm E}_{\mathbf{F}}[\bm{\beta}]$ $\displaystyle=$
$\displaystyle{\rm E}[\bm{\beta}]+{\rm Cov}[\bm{\beta},\mathbf{F}]{\rm
Var}[\mathbf{F}]^{-1}(\mathbf{F}-{\rm E}[\mathbf{F}])$ (146) $\displaystyle=$
$\displaystyle\bm{\Gamma}+\bm{\mathit{\Delta}}\,\bm{W}^{T}\,(\bm{W}\,\bm{\mathit{\Delta}}\,\bm{W}^{T}+\bm{\mathit{\Omega}})^{-1}(\mathbf{F}-\bm{W}\bm{\Gamma})$
$\displaystyle=$
$\displaystyle\bm{\Gamma}+(\bm{W}^{T}\,\bm{\mathit{\Omega}}^{-1}\,\bm{W}+\bm{\mathit{\Delta}}^{-1})^{-1}\bm{W}^{T}\,\bm{\mathit{\Omega}}^{-1}(\mathbf{F}-\bm{W}\bm{\Gamma})$
$\displaystyle=$
$\displaystyle(\bm{W}^{T}\,\bm{\mathit{\Omega}}^{-1}\,\bm{W}+\bm{\mathit{\Delta}}^{-1})^{-1}\big{(}(\bm{W}^{T}\,\bm{\mathit{\Omega}}^{-1}\,\bm{W}+\bm{\mathit{\Delta}}^{-1})\bm{\Gamma}+\bm{W}^{T}\,\bm{\mathit{\Omega}}^{-1}(\mathbf{F}-\bm{W}\bm{\Gamma})\big{)}$
$\displaystyle=$
$\displaystyle(\bm{W}^{T}\,\bm{\mathit{\Omega}}^{-1}\,\bm{W}+\bm{\mathit{\Delta}}^{-1})^{-1}(\bm{\mathit{\Delta}}^{-1}\bm{\Gamma}+\bm{W}^{T}\,\bm{\mathit{\Omega}}^{-1}\,\mathbf{F})$
$\displaystyle=$
$\displaystyle(\bm{W}^{T}\,\bm{\mathit{\Omega}}^{-1}\,\bm{W}+\bm{\mathit{\Delta}}^{-1})^{-1}(\bm{W}^{T}\,\bm{\mathit{\Omega}}^{-1}\,\bm{W}\,\hat{\bm{\beta}}_{GLS}+\bm{\mathit{\Delta}}^{-1}\bm{\Gamma})$

where line 3 holds as a result of the identity (see, for example, Mardia et
al. (1979))

$AB(DAB+C)^{-1}=(BC^{-1}D+A^{-1})^{-1}BC^{-1}$ (147)

where here, $A,B,C,D$ are generic matrices, and

$\displaystyle\hat{\bm{\beta}}_{GLS}$ $\displaystyle=$
$\displaystyle(\bm{W}^{T}\,\bm{\mathit{\Omega}}^{-1}\,\bm{W})^{-1}\bm{W}^{T}\,\bm{\mathit{\Omega}}^{-1}\,\mathbf{F}$
(148) $\displaystyle=$
$\displaystyle\bm{I}_{q}\otimes(\bm{G}^{T}\,\bm{C}^{-1}\,\bm{G})^{-1}\bm{G}^{T}\,\bm{C}^{-1}\,\mathbf{F}$

is the generalised least squares estimate for $\bm{\beta}$. Combining
Equations (144) and (145) we get:

$\displaystyle{\rm E}_{\mathbf{F}}[\mathbf{f}(\mathbf{X})]$ $\displaystyle=$
$\displaystyle{\rm E}[\bm{w}(\mathbf{X})]\,{\rm
E}_{\mathbf{F}}[\bm{\beta}]+\bm{v}(\mathbf{X})\,\bm{\mathit{\Omega}}^{-1}\,(\mathbf{F}-\bm{W}\,{\rm
E}_{\mathbf{F}}[\bm{\beta}])\,.$

$\Box$

### Proof of Lemma 4

We begin by expanding the terms of the Bayes linear update as follows:

$\displaystyle{\rm Var}_{\mathbf{F}}[\mathbf{f}(\mathbf{X})]$ $\displaystyle=$
$\displaystyle{\rm
Var}_{\mathbf{F}}[\bm{w}(\mathbf{X})\,\bm{\beta}+\mathbf{u}(\mathbf{X})]$
(149) $\displaystyle=$ $\displaystyle{\rm
Var}_{\mathbf{F}}[\bm{w}(\mathbf{X})\,\bm{\beta}]+{\rm
Var}_{\mathbf{F}}[\mathbf{u}(\mathbf{X})]$ $\displaystyle\,\,\,\,\,\,\,+\,{\rm
Cov}_{\mathbf{F}}\left[\bm{w}(\mathbf{X})\,\bm{\beta},\mathbf{u}(\mathbf{X})\right]+{\rm
Cov}_{\mathbf{F}}\left[\mathbf{u}(\mathbf{X}),\bm{w}(\mathbf{X})\,\bm{\beta}\right]\,.$

Starting with the first term on the right-hand side of Equation (149), we
have, using linear algebra (Mardia et al. 1979):

$\displaystyle{\rm Var}_{\mathbf{F}}[\bm{w}(\mathbf{X})\,\bm{\beta}]$
$\displaystyle=$ $\displaystyle{\rm E}[{\rm
Var}_{\mathbf{F}\cup\mathbf{X}}[\bm{w}(\mathbf{X})\,\bm{\beta}]]+{\rm
Var}[{\rm E}_{\mathbf{F}\cup\mathbf{X}}[\bm{w}(\mathbf{X})\,\bm{\beta}]]$
(150) $\displaystyle=$ $\displaystyle{\rm E}[\bm{w}(\mathbf{X})\,{\rm
Var}_{\mathbf{F}}[\bm{\beta}]\,\bm{w}(\mathbf{X})^{T}]+{\rm
E}_{\mathbf{F}}[\bm{\beta}^{T}]\,{\rm Var}[\bm{w}(\mathbf{X})]{\rm
E}_{\mathbf{F}}[\bm{\beta}]\,,$

where the first term of Equation 150 can be calculated component-wise:

$\displaystyle{\rm E}_{ij}[\bm{w}(\mathbf{X}){\rm
Var}_{\mathbf{F}}[\bm{\beta}]\bm{w}(\mathbf{X})^{T}]$ $\displaystyle=$
$\displaystyle{\rm E}[\textrm{trace}(\bm{w}_{i}(\mathbf{X})\,{\rm
Var}_{\mathbf{F}}[\bm{\beta}]\,\bm{w}_{j}(\mathbf{X})^{T})]$ $\displaystyle=$
$\displaystyle{\rm E}[\textrm{trace}({\rm
Var}_{\mathbf{F}}[\bm{\beta}]\,\bm{w}_{j}(\mathbf{X})^{T}\,\bm{w}_{i}(\mathbf{X}))]$
$\displaystyle=$ $\displaystyle\textrm{trace}({\rm E}[{\rm
Var}_{\mathbf{F}}[\bm{\beta}]\,\bm{w}_{j}(\mathbf{X})^{T}\,\bm{w}_{i}(\mathbf{X})])$
$\displaystyle=$ $\displaystyle\textrm{trace}({\rm
Var}_{\mathbf{F}}[\bm{\beta}]{\rm
E}[\bm{w}_{j}(\mathbf{X})^{T}\,\bm{w}_{i}(\mathbf{X})])$ $\displaystyle=$
$\displaystyle\textrm{trace}\big{(}{\rm Var}_{\mathbf{F}}[\bm{\beta}]({\rm
E}[\bm{w}_{j}(\mathbf{X})^{T}]{\rm E}[\bm{w}_{i}(\mathbf{X})])+{\rm
Cov}[\bm{w}_{j}(\mathbf{X}),\bm{w}_{i}(\mathbf{X})]\big{)}$

which is now in terms of specified quantities. Note that ${\rm
Var}_{\mathbf{F}}[\bm{\beta}]$ is given by

$\displaystyle{\rm Var}_{\mathbf{F}}[\bm{\beta}]$ $\displaystyle=$
$\displaystyle{\rm Var}[\bm{\beta}]-{\rm Cov}[\bm{\beta},\mathbf{F}]{\rm
Var}[\mathbf{F}]^{-1}{\rm Cov}[\mathbf{F},\bm{\beta}]$ (152) $\displaystyle=$
$\displaystyle\bm{\mathit{\Delta}}-\bm{\mathit{\Delta}}\,\bm{W}^{T}(\bm{W}\,\bm{\mathit{\Delta}}\,\bm{W}^{T}+\bm{\mathit{\Omega}})^{-1}\bm{W}\,\bm{\mathit{\Delta}}$
$\displaystyle=$
$\displaystyle\bm{\mathit{\Delta}}-(\bm{W}^{T}\,\bm{\mathit{\Omega}}^{-1}\,\bm{W}+\bm{\mathit{\Delta}}^{-1})^{-1}\bm{W}^{T}\,\bm{\mathit{\Omega}}^{-1}\,\bm{W}\,\bm{\mathit{\Delta}}$
$\displaystyle=$
$\displaystyle(\bm{W}^{T}\,\bm{\mathit{\Omega}}^{-1}\,\bm{W}+\bm{\mathit{\Delta}}^{-1})^{-1}\big{(}(\bm{W}^{T}\,\bm{\mathit{\Omega}}^{-1}\,\bm{W}+\bm{\mathit{\Delta}}^{-1})\bm{\mathit{\Delta}}-\bm{W}^{T}\,\bm{\mathit{\Omega}}^{-1}\,\bm{W}\,\bm{\mathit{\Delta}}\big{)}$
$\displaystyle=$
$\displaystyle(\bm{W}^{T}\,\bm{\mathit{\Omega}}^{-1}\,\bm{W}+\bm{\mathit{\Delta}}^{-1})^{-1}$

For the second term of Equation (149), we have that

$\displaystyle{\rm Var}_{\mathbf{F}}[\mathbf{u}(\mathbf{X})]$ (153)
$\displaystyle=$ $\displaystyle{\rm Var}[\mathbf{u}(\mathbf{X})]-{\rm
Cov}[\mathbf{u}(\mathbf{X}),\mathbf{F}]\,{\rm Var}[\mathbf{F}]^{-1}\,{\rm
Cov}[\mathbf{F},\mathbf{u}(\mathbf{X})]$ $\displaystyle=$
$\displaystyle\bm{\mathit{\Sigma}}-\bm{v}(\mathbf{X})\,(\bm{W}\,\bm{\mathit{\Delta}}\,\bm{W}^{T}+\bm{\mathit{\Omega}})^{-1}\,\bm{v}(\mathbf{X})^{T}$
$\displaystyle=$