SEPIA model math

There are \(n\) physical (observational) experiments. From the \(i\) inputs \({\bf x}^{obs}_i=(x^{obs}_{i1}, \ldots, x^{obs}_{ip})\), the observation \({\bf y}^{obs}_{i}({\bf x}^{obs}_i)\) (an \(n_{y^{obs}_{i}} \times 1\) vector) is modeled by

\[{\bf y}^{obs}_{i}({\bf x}^{obs}_{i})= \boldsymbol \eta({\bf x}^{obs}_{i},\boldsymbol \theta)+ \boldsymbol \delta({\bf x}^{obs}_{i}) + {\bf e}^{obs}_{i},\]

where the observation error vector \({\bf e}^{obs}_{i}\) is modeled by

\[{\bf e}^{obs}_{i} \sim MVN\left({\bf 0}_{n_{y^{obs}_{i}}}, \, \frac{1}{\lambda_{y^{obs}}^{\tt Os}} \Sigma^{obs}_i \right).\]

\(\boldsymbol \eta(\cdot)\) is an emulator from a simulation code, \(\boldsymbol \theta\) corresponds to inputs of the parameter, and \(\boldsymbol \delta(\cdot)\) is a discrepancy from reality.

There are \(m\) simulation experiments. From the \(i\) inputs \({\bf x}^{sim}_i=(x^{sim}_{i1}, \ldots, x^{sim}_{ip})\) and \({\bf t}^{sim}_i=(t^{sim}_{i1}, \ldots, t^{sim}_{iq})\), the observation \({\bf y}^{sim}_{i}({\bf x}^{sim}_i,{\bf t}^{sim}_i)\) (an \(n_{y^{sim}_{i}} \times 1\) vector) is modeled by

\[{\bf y}^{sim}_{i}({\bf x}^{sim}_i,{\bf t}^{sim}_i)= \boldsymbol \eta({\bf x}^{sim}_i,{\bf t}^{sim}_i)+ {\bf e}^{sim}_{i},\]

where the error vector \({\bf e}^{sim}_{i}\) is modeled by \(MVN\left({\bf 0}_{n_{y^{sim}_{i}}}, \, \frac{1}{\lambda^{\tt WOs}_{y^{sim}}} {\bf I} \right)\) and \(I_m\) is the \(m \times m\) identity matrix.

We re-express \(\boldsymbol \eta({\bf x}_i,{\bf t}_i)\) and \(\boldsymbol \delta({\bf x}_i)\) by linear combinations of basis functions and approximate them using a subset of the complete set of basis functions. Consequently,

\[\boldsymbol \eta({\bf x}^{obs}_i,\boldsymbol \theta) \approx \sum_{j=1}^{p_u} {\bf K}^{obs}_j u_j({\bf x}^{obs}_i, \boldsymbol \theta)\]

for \(p_{u}\) basis functions \({\bf K}^{obs}_j\). So the matrix \({\bf K}^{obs}=({\bf K}^{obs}_1 \cdots {\bf K}^{obs}_{p_u})\). Similarly,

\[\boldsymbol \delta({\bf x}^{obs}_i) \approx \sum_{j=1}^{p_v} {\bf D}^{obs}_j v_j({\bf x}^{obs}_i)\]

for \(p_{v}\) basis functions \({\bf D}^{obs}_j\). So the matrix \({\bf D}^{obs}=({\bf D}^{obs}_1 \cdots {\bf D}^{obs}_{p_v})\).

For the simulations,

\[\boldsymbol \eta({\bf x}^{sim}_i,{\bf t}^{sim}_i) \approx \sum_{j=1}^{p_u} {\bf K}^{sim}_j w_j({\bf x}^{sim}_i,{\bf t}^{sim}_i)\]

for \(p_{u}\) basis functions \({\bf K}^{sim}_j\), where \(w_j({\bf x}^{sim}_i,{\bf t}^{sim}_i)=u_j({\bf x}^{sim}_i,{\bf t}^{sim}_i)+\epsilon^{sim,nug}_j\). So the matrix \({\bf K}^{sim}=({\bf K}^{sim}_1 \cdots {\bf K}^{sim}_{p_u})\).

Note that \(\frac{1}{\lambda^{\tt Ws}_{\epsilon^{sim,nug}_j}}\) is the variance of an i.i.d. Normally distributed nugget \(\epsilon^{sim,nug}_j\) with mean 0 to account for small numerical fluctuations in the simulator. The nugget is used only in fitting, not in prediction.

In the above equations, any error from the truncated basis approximations is assumed to be part of \({\bf e}^{obs}_{i}\) or \({\bf e}^{sim}_{i}\).

The \(u_j({\bf x},{\bf t}), \, j=1, \ldots, p_u\) are modeled as a GP with mean \({\bf 0}_{n}\) and variance covariance matrix \(\frac{1}{\lambda^{\tt Uz}_{u_j}}R^{u}_j\), where

\[R^{u}_j(({\bf x}_i,{\bf t}_i)),({\bf x}_l,{\bf t}_l))=\prod_{k=1}^p \left({\rho^{u}_{jk}}\right)^{4|x_{ik}-x_{lk}|^2} \prod_{k=1}^q \left({\rho^{u}_{(j+p)k}}\right)^{4 |t_{ik}-t_{lk}|^2}.\]

A more familiar form (revealing the squared exponential covariance function form) is

\[R^{u}_j(({\bf x}_i,{\bf t}_i)),({\bf x}_l,{\bf t}_l))=\prod_{k=1}^p \exp(-{{\beta^{u}_{jk}}}|x_{ik}-x_{lk}|^2) \prod_{k=1}^q \exp(-{{\beta^{u}_{(j+p)k}}}|x_{ik}-x_{lk}|^2),\]

so that \(\beta^{u}_{jk}= -{4}\log\left(\rho^u_{jk}\right)\) or \(\rho^u_{jk}=\exp\left(-\frac{\beta^u_{jk}}{4}\right)\).

Similarly, the \(v_j({\bf x}^{obs}_i), \, j=1, \ldots, n\) are modeled as a GP with mean \({\bf 0}_{n}\) and variance covariance matrix \(\frac{1}{\lambda^{\tt Vz}_{v^{obs}_j}}R^{v}_j\), where

\[R^{v}_j({\bf x}^{obs}_i,{\bf x}^{obs}_l)=\prod_{k=1}^p ({\rho^{v}_{jk}})^{4|x_{ik}-x_{lk}|^2},\]

whose more familiar form is

\[R^{v}_j({\bf x}^{obs}_i,{\bf x}^{obs}_l)=\prod_{k=1}^q \exp(-{{\beta^{v}_{jk}}}|x_{ik}-x_{lk}|^2),\]

so that \(\beta^{v}_{jk}= -4\log\left(\rho^v_{jk}\right)\) or \(\rho^v_{jk}=\exp\left(-\frac{\beta^v_{jk}}{4}\right)\).

Note that the \({\bf x}, {\bf t}, \boldsymbol \theta\) are transformed to [0, 1] and the \({\bf y}^{obs}_{i}({\bf x}^{obs}_{i})\) and \({\bf y}^{sim}_{i}({\bf x}^{sim}_i,{\bf t}^{sim}_i)\) are normalized to have sample mean \({\bf 0}\) and covariance matrices equal to identity matrices. Consequently, the \(\Sigma^{obs}_i\) for \({\bf y}^{obs}_{i}({\bf x}^{obs}_{i})\) has to be normalized in the same way that the \({\bf y}^{obs}_{i}({\bf x}^{obs}_{i})\) are.