## 2 The Surrogate Model for Computer Experiments
### The Kriging Model  
The kriging model was originally introduced in Geostatistics to estimate the properties of the most likely distribution of gold based on a given set of sampled sites [17]. The kriging model is also known as the Gaussian process model. The kriging model is an interpolation method for which the unknown predictions are modeled by the Gaussian process, which is expressed by a prior covariance. Stated differently, the unknown predictions (\(\hat{g}(\mathbf{x}_{\mathbf{new}})\)) are modeled as a realization of the Gaussian process with the covariance conditional on the training samples. The kriging model can be expressed as  
\[\hat{g}(\mathbf{x})=\mathbf{f}(\mathbf{x})+\mathbf{Z}(\mathbf{x})\ \text{with}\ \mathbf{Z}(\mathbf{x})\ \sim\ GP \big{(}\mathbf{0},\ \mathbf{C}\big{(}\mathbf{x},\mathbf{x}^{\prime}\big{)}\big{)} \tag{1}\]  
where \(\mathbf{x}\in R^{k}\) is a sample with \(k\) input variables; \(\mathbf{f}(\mathbf{x})\) can be constant or polynomial function to approximate a global trend over the entire parameter space; \(\mathbf{Z}(\mathbf{x})\) is a stochastic component to represent a local spatial deviation and assumed to be a stationary Gaussian process with a zero-mean and covariance (\(\mathbf{C}(\cdot,\cdot)\)) as  
\[C\big{(}\mathbf{Z}(\mathbf{x}_{i}),\mathbf{Z}(\mathbf{x}_{j})\big{)}=\sigma^{2}\psi\big{(}\bm {x}_{i},\mathbf{x}_{j}\big{)} \tag{2}\]  
where \(\sigma^{2}\) is the process variance; the subscript \(i\) and \(j\) indicate the \(i\)th and \(j\)th training sample, respectively; \(\psi\big{(}\mathbf{x}_{i},\mathbf{x}_{j}\big{)}\) is a spatial correlation function to control the smoothness of the kriging model. A typical choice of \(\psi\big{(}\mathbf{x}_{i},\mathbf{x}_{j}\big{)}\) is the \(k\)-dimensional Gaussian correlation function given by the following:  
\[\psi\big{(}\mathbf{x}_{i},\mathbf{x}_{j}\big{)}=\exp\left(-\sum_{p=1}^{k}\theta_{p} \|x_{i}^{(p)}-x_{j}^{(p)}\|^{2}\right) \tag{3}\]  
where \(\|x_{i}^{(p)}-x_{j}^{(p)}\|^{2}\) is the Euclidean distance measure between two samples. \(\mathbf{\theta}=\big{[}\theta_{1},\cdots,\theta_{p}\big{]}\) is an unknown correlation parameter vector to scale the correlation length in each input. The correlation length also reflects the relative significance of each input [2]. Different types of the kriging model can be modeled as follows: (1) the ordinary kriging model assumes a constant mean only (\(\mathbf{f}(\mathbf{x})=\mu(\mathbf{x})\)) over the neighborhood of \(\mathbf{x}\), and (2) the universal kriging model assumes a general polynomial trend model (\(\mathbf{f}(\mathbf{x})=\sum_{i=0}^{k}\beta_{i}x_{i}\)). Since the ordinary kriging model is used in this study, the ordinary kriging model is briefly described. Hereafter, the kriging model denotes the ordinary kriging model for simplicity. Suppose that we have \(n\) observed inputs \(\mathbf{X}_{\mathbf{n}}=[\mathbf{x}_{1},\cdots,\mathbf{x}_{n}]^{T}\) and the corresponding outputs \(\mathbf{Y}_{\mathbf{n}}=[y(\mathbf{x}_{1}),\cdots,y(\mathbf{x}_{n})]^{T}\). To construct an ordinary kriging model using \(\mathbf{X}_{\mathbf{n}}\) and \(\mathbf{Y}_{\mathbf{n}}\), the unknown parameters of the kriging model (\(\mu\), \(\sigma^{2}\) and \(\mathbf{\theta}\)) are estimated by maximum likelihood estimation. The log-likelihood function (Equation (4)) is expressed with these unknown parameters:  
\[\ln\Bigl{(}L\Bigl{(}\mu,\sigma^{2},\mathbf{\theta}\big{|}\mathbf{X}_{\mathbf{n}},\mathbf{Y}_{ \mathbf{n}}\Bigr{)}\Bigr{)}\approx-\frac{n}{2}\ln\Bigl{(}\sigma^{2}\Bigr{)}-\frac {1}{2}\ln\lvert\mathbf{\Psi}\rvert-\frac{(\mathbf{Y}_{\mathbf{t}}-\mathbf{1}\mu)^{T}\mathbf{\Psi} ^{-1}(\mathbf{Y}_{\mathbf{t}}-\mathbf{1}\mu)}{2\sigma^{2}} \tag{4}\]  
Figure 1: An illustration of surrogate modeling. where \(\boldsymbol{\
abla}\) is an \(n~{}\times~{}n\) correlation matrix with the element \(\psi\big{(}\boldsymbol{x}_{i},\boldsymbol{x}_{j}\big{)}\) and \(\boldsymbol{\varOmega}\) is an n-by-1 unit vector.