Best Linear Unbiased Prediction for Multifidelity Computer Experiments

Hindawi Mathematical Problems in Engineering Volume 2018, Article ID 8525736, 7 pages https://doi.org/10.1155/2018/8525736 Research Article Best Linear Unbiased Prediction for Multifidelity Computer Experiments Weiyan Mu 1 ,1 Qiuyue Wei,1 Dongli Cui,1 and Shifeng Xiong2 School of Science, Beijing University of Civil Engineering and Architecture, Beijing 100044, China NCMIS, Academy of Mathematics and Systems Science, Chinese Academy of Sciences, Beijing 100190, China 2 Correspondence should be addressed to Weiyan Mu; Received 26 September 2017; Revised 26 April 2018; Accepted 8 May 2018; Published 7 June 2018 Academic Editor: Elisa Francomano Copyright © 2018 Weiyan Mu et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. Recently it becomes a growing trend to study complex systems which contain multiple computer codes with different levels of accuracy, and a number of hierarchical Gaussian process models are proposed to handle such multiple-fidelity codes. This paper derives the best linear unbiased prediction for three popular classes of multiple-level Gaussian process models. The predictors all have explicit expressions at each untried point. Empirical best linear unbiased predictors are also provided by plug-in methods with generalized maximum likelihood estimators of unknown parameters. 1. Introduction With the rapid development of computer technology, computer experiments have been widely used in engineering and science [1]. Consider a computer code with inputs x ∈ [0, 1]𝑑 . The common statistical modeling method is to model the response of the code as 𝑝 𝑌 (x) = ∑𝑓𝑗 (x) 𝛽𝑗 + 𝑍 (x) = 𝑓󸀠 (x) 𝛽 + 𝑍 (x) , (1) 𝑗=1 where 𝑓1 (⋅), . . . , 𝑓𝑝 (⋅) are known regression functions, 𝛽 = (𝛽1 , . . . , 𝛽𝑝 )󸀠 is a vector of unknown regression coefficients, and 𝑍(⋅) is a stationary Gaussian process on [0, 1]𝑑 having zero mean, variance 𝜎𝑧2 , and correlation function 𝑅(⋅). Let Y = (𝑦1 , . . . , 𝑦𝑛 )󸀠 be the responses corresponding to the design points x1 , . . . , x𝑛 ∈ [0, 1]𝑑 . With a known correlation function 𝑅(⋅), the best linear unbiased predictor (BLUP) is given by ̂ ̂ + r󸀠 𝑅−1 (Y𝑛 − F𝛽) ̂ (x0 ) = 𝑌 ̂0 ≡ f0 󸀠 𝛽 𝑌 0 (2) for x0 ∈ [0, 1]𝑑 , where r0 = (𝑅(x0 − x1 ), . . . , 𝑅(x0 − x𝑛 ))󸀠 , ̂ = R = (𝑅(x𝑖 − x𝑗 ))𝑖,𝑗=1,...,𝑛 , F = (f(x1 ), . . . , f(x𝑛 ))󸀠 , and 𝛽 (F󸀠 R−1 F)−1 F󸀠 R−1 Y is the generalized least squares estimator of 𝛽 [1]. It is a growing trend to study the complex system which contains multifidelity computer codes with different levels of accuracy. For example, a Bayesian approach is described to predict and analyze complex computer codes which can be run at different levels of sophistication [2]. A novel approach is taken to integrate data from approximate and detailed simulations to build a surrogate model that describes the relationship between output and input parameters [3]. The Bayesian hierarchical Gaussian process (BHGP) models are introduced to integrate low-accuracy and high-accuracy [4]. A class of nonstationary Gaussian process models are proposed to link the computer outputs of different mesh densities [5]. However, there are few papers that provide BLUPs for multifidelity computer experiments. The purpose of this article is to find BLUPs for multifidelity computer experiments. The structure of the article is as follows. In Section 2, BLUPs for two levels of accuracy [4] are discussed. In Section 3, BLUPs for general 𝑘-level cases in autoregressive model described by Kennedy and O’Hagan [2] are illustrated. In Section 4, BLUPs for continuous level in nonstationary Gaussian process model described by Tuo, Wu, and Yu [5] are demonstrated. We present a real application in Section 5. Concluding remarks are given in Section 6. 2 Mathematical Problems in Engineering 2. BLUPs for Two-Level Cases + 𝜎𝑍2 a󸀠1 R𝜃11 a1 + 𝜎𝑍2 𝜎𝜌2 a󸀠2 V1 a2 + 𝜌02 𝜎𝑍2 a󸀠2 R𝜃21 a2 The two experiments considered in this section are named as the low-accuracy experiment (LE) and high-accuracy experiment (HE). Let GP(𝜇, 𝜎2 , 𝜃) denote the Gaussian process with mean 𝜇, variance 𝜎2 , and correlation parameters 𝜃. Let 𝐷𝑙 = (x1 , . . . , x𝑛 ) and 𝐷ℎ = (s1 , . . . , s𝑚 ) denote the design set for the LE and HE, respectively. Following Qian and Wu [4], for any x𝑖 ∈ 𝐷𝑙 , the LE is described by + (𝜎𝜌2 + 𝜌02 ) 𝜎𝑍2 + 2𝜌0 𝜎𝑍2 a󸀠1 R𝜃31 a2 − 2𝜌0 𝜎𝑍2 a󸀠1 r1𝜃1 𝑌𝑙 (x𝑖 ) = f 󸀠 (x𝑖 ) 𝛽 + 𝑍 (x𝑖 ) , (3) where 𝑍(⋅) ∼ GP(0, 𝜎𝑍2 , 𝜃1 ) and f(x) = [𝑓1 (x), . . . , 𝑓𝑞 (x)]󸀠 is a set of prespecified regression functions. For any s𝑖 ∈ 𝐷ℎ , the HE can be described by 𝑌ℎ (s𝑖 ) = 𝜌 (s𝑖 ) 𝑦𝑙 (s𝑖 ) + 𝛿 (s𝑖 ) + 𝜀 (s𝑖 ) , (4) where the scale changes from LE to HE 𝜌(⋅) ∼ GP(𝜌0 , 𝜎𝜌2 , 𝜃3 ), the location adjustment 𝛿(⋅) ∼ GP(𝛿0 , 𝜎𝛿2 , 𝜃2 ), the measurement error 𝜀 ∼ 𝑁(0, 𝜎𝜀2 ) and 𝑍(⋅), 𝛿(⋅), and 𝜀 are jointly independent. Let Y𝑙 = (𝑌𝑙 (x1 ), . . . , 𝑌𝑙 (x𝑛 ))󸀠 and Yℎ = (𝑌𝑙 (s1 ), . . . , 𝑌𝑙 (s𝑚 ))󸀠 . 𝑑 − 2𝜎𝑍2 𝜎𝜌2 a󸀠2 k1 − 2𝜌02 𝜎𝑍2 a󸀠2 r2𝜃1 + 𝜎𝛿2 (a󸀠2 R𝜃2 a2 − 2a󸀠2 r𝜃2 + 1) + 𝜎𝜖2 (a󸀠2 a2 + 1) , (8) where V = (𝑅𝜃3 (s𝑖 − s𝑗 )f(s𝑖 )f 󸀠 (s𝑗 )), V1 = (𝑅𝜃1 (s𝑖 − s𝑗 )𝑅𝜃3 (s𝑖 − s𝑗 )), R𝜃21 = (𝑅𝜃1 (s𝑖 − s𝑗 )), R𝜃2 = (𝑅𝜃2 (s𝑖 − s𝑗 )), 1 ≤ 𝑖, 𝑗 ≤ 𝑚; k = (𝑅𝜃3 (x0 − s𝑖 )f 󸀠 (s𝑖 )), k1 = (𝑅𝜃1 (x0 − s𝑖 )𝑅𝜃3 (x0 − s𝑖 )), r2𝜃1 = (𝑅𝜃1 (x0 − s𝑖 ))󸀠 , r𝜃2 = (𝑅𝜃2 (x0 − s𝑖 ))󸀠 , 1 ≤ 𝑖 ≤ 𝑚; R𝜃11 = (𝑅𝜃1 (x𝑖 − x𝑗 )),1 ≤ 𝑖, 𝑗 ≤ 𝑛; R𝜃31 = (𝑅𝜃1 (x𝑖 − s𝑗 )),1 ≤ 𝑖 ≤ 𝑚, 1 ≤ 𝑗 ≤ 𝑛; r1𝜃1 = (𝑅𝜃1 (x0 − x𝑖 ))󸀠 , 1 ≤ 𝑖 ≤ 𝑛. Thus the Lagrange multipliers can be used to solve BLUP corresponding to a1 and a2 that minimize (5) subject to a󸀠1 F1 + 𝜌0 a󸀠2 F2 = 𝜌0 f 󸀠 (x0 ) , a󸀠2 1𝑚 = 1. The Lagrange function is 󸀠 Theorem 1. For x0 ∈ [0, 1] , the BLUP of 𝑌ℎ (x0 ) is a Y, where 𝐿 = 𝜎𝜌2 (a󸀠2 𝛽󸀠 V𝛽a2 + f 󸀠 (x0 ) 𝛽𝛽󸀠 f (x0 ) a −1 −1 = B−1 A (A󸀠 B−1 A) k1 + [B−1 − B−1 A (A󸀠 B−1 A) A󸀠 B−1 ] k2 , − 2a󸀠2 k𝛽𝛽󸀠 f (x0 )) + 𝜎𝑍2 a󸀠1 R𝜃11 a1 𝜎𝑍2 𝜎𝜌2 a󸀠2 V1 a2 (5) + 𝜌02 𝜎𝑍2 a󸀠2 R𝜃21 a2 + (𝜎𝜌2 + 𝜌02 ) 𝜎𝑍2 + 2𝜌0 𝜎𝑍2 a󸀠1 R𝜃31 a2 F1 0 A=[ ], 𝜌0 F2 1𝑚 k1 = [ 𝜌0 f (x0 ) 1 (9) − 2𝜌0 𝜎𝑍2 a󸀠1 r1𝜃1 − 2𝜎𝑍2 𝜎𝜌2 a󸀠2 k1 − 2𝜌02 𝜎𝑍2 a󸀠2 r2𝜃1 (10) + 𝜎𝛿2 (a󸀠2 R𝜃2 a2 − 2a󸀠2 r𝜃2 + 1) + 𝜎𝜖2 (a󸀠2 a2 + 1) ], 𝜌0 𝜎𝑍2 r1𝜃1 k2 = [ 2 ], 𝜎𝜌 k𝛽𝛽󸀠 f (x0 ) + 𝜎𝜌2 𝜎𝑍2 k1 + 𝜌02 𝜎𝑍2 r2𝜃1 + 𝜎𝛿2 r𝜃2 + 2𝜆󸀠1 (F󸀠1 a1 + 𝜌0 F󸀠2 a2 − 𝜌0 f (x0 )) + 2𝜆 2 (1󸀠𝑚 a2 (6) Let the gradient with respect to a1 , a2 , 𝜆1 , 𝜆 2 be zero, and we have 𝜆 k1 0 A󸀠 [ (11) ] [ ] = [ ], k2 a A B and B 𝜎𝑍2 R𝜃11 𝜌0 𝜎𝑍2 R𝜃31 ]. =[ 2 3 󸀠 2 󸀠 𝜌𝜎 R 𝜎𝜌 𝛽 V𝛽 + 𝜎𝜌2 𝜎𝑍2 V1 + 𝜌02 𝜎𝑍2 R𝜃21 + 𝜎𝛿2 R𝜃2 + 𝜎𝜖2 I𝑚 [ 0 𝑍 𝜃1 ] ̂ℎ (x0 ) = a󸀠 Y𝑙 + a󸀠 Yℎ + 𝐶 based Proof. The linear predictor 𝑌 1 2 on training data Y𝐿 and Y𝐻 at an untried point x0 is unbiased for 𝑌ℎ (x0 ) provided a󸀠1 F1 + 𝜌0 a󸀠2 F2 = 𝜌0 f 󸀠 (x0 ) , a󸀠2 1𝑚 = 1, (7) 𝐶 = 0. For any linear unbiased predictor (LUP) of 𝑌0 = 𝑌ℎ (x0 ), ̂ℎ (x0 ) = a󸀠 Y𝑙 + a󸀠 Yℎ , the mean squared prediction error 𝑌 1 2 (MSPE) of a󸀠1 Y𝑙 + a󸀠2 Yℎ is 𝐸{𝑀2 } = 𝐸{(a󸀠1 Y𝑙 + a󸀠2 Yℎ − 𝑌0 )2 }. We have 𝐸 {𝑀2 } = 𝜎𝜌2 (a󸀠2 𝛽󸀠 (...truncated)