arXiv:1703.01326v1 [] 3 Mar 2017

Prediction based on the Kennedy-O’Hagan calibration
model: asymptotic consistency and other properties
Rui Tuo
Academy of Mathematics and Systems Sciences
Chinese Academy of Sciences
C. F. Jeff Wu
School of Industrial and Systems Engineering
Georgia Institute of Technology
March 7, 2017
Abstract
Kennedy and O’Hagan (2001) propose a model for calibrating some unknown parameters in a computer model and estimating the discrepancy between the computer
output and physical response. This model is known to have certain identifiability
issues. Tuo and Wu (2016) show that there are examples for which the KennedyO’Hagan method renders unreasonable results in calibration. In spite of its unstable
performance in calibration, the Kennedy-O’Hagan approach has a more robust behavior in predicting the physical response. In this work, we present some theoretical
analysis to show the consistency of predictor based on their calibration model in the
context of radial basis functions.
Key words and phrases: computer experiments, kriging, Bayesian inference

1

Introduction

With the development of mathematical models and computational technique, simulation programs or softwares are shown to be increasingly powerful for the prediction,
validation and control of many physical processes. A computer simulation run, based
on a virtual platform, requires only computational resources that are rather inexpensive in today’s computing environment. In contrast, a physical experiment usually
requires more facilities, materials, and human labor. As a consequence, a typical
computer simulation run is much cheaper than its corresponding physical experiment

1

trial. The economic benefits of computer simulations make them particularly useful
and attractive in scientific and engineering research. As a branch of statistics, design
of experiments mainly studies the methodologies on the planning, analysis and optimization of physical experiments (Wu and Hamada, 2011). Given the rapid spread of
computer simulations, it is beneficial to develop theory and methods for the design and
analysis of computer simulation experiments. This emerging field is commonly referred
to as computer experiments. We refer to Santner et al. (2003) for more details.
The input variables of a computer experiment normally consist of factors which
can be controlled in the physical process, referred to as the control variables, as well as
some model parameters. These model parameters represent certain intrinsic properties
of the physical system. For example, to simulate a heat transfer process, we need to
solve a heat equation. The formulation of the equation requires the environmental settings and the initial conditions of the system which can be controlled physically, as well
as the thermal conductivity which is uncontrollable and cannot be measured directly
in general. For most computer simulations, the prediction accuracy of the computer
model is closely related to the choice of the model parameters. A standard method
for determining the unknown model parameters is to estimate them by comparing the
computer outputs and the physical responses. Such a procedure is known as the calibration for computer models, and the model parameters to be identified are called
the calibration parameters. Kennedy and O’Hagan (2001) first study the calibration
problem using ideas and methods in computer experiments. They propose a Bayesian
hierarchical model to estimate the calibration parameters by computing their posterior distributions. Tuo and Wu (2016) show that the Kennedy-O’Hagan method may
render unreasonable estimates for the calibration parameters. Given the widespread
use of the Kennedy-O’Hagan method, it will be desirable to make a comprehensive
assessment about this method. For brevity, we sometimes refer to Kennedy-O’Hagan
as KO.
This paper endeavors to study the prediction performance of the Kennedy-O’Hagan
approach. First we adopt the framework of Tuo and Wu (2016) which assumes the
physical observations to be non-random. Interpolation theory in the native spaces
becomes the key mathematical tool in this part. Then, we study the more realistic
situation where the physical data are noisy. We employ the asymptotic theory of the
smoothing splines in the Sobolev spaces to obtain the rate of convergence of the KO
predictor in this case.
This article is organized as follows. In Section 2 we review the Bayesian method
proposed by Kennedy and O’Hagan (2001) for calibrating the model parameters and
predicting for new physical responses. In Section 3 we present our main results on
the asymptotic theory on the prediction performance of the KO method. Concluding
remarks and further discussions are made in Section 4. Some technical proofs are given
in Appendix A.

2

2

Review on the Kennedy-O’Hagan Method

In this section we review the Bayesian method proposed by Kennedy and O’Hagan
(2001). The formulation of this approach can be generalized to some extend. See, for
example, Higdon et al. (2004).
Denote the experimental region for the control variables as Ω. We suppose that Ω
is a convex and compact subset of Rd . Let {x1 , . . . , xn } ⊂ Ω be the set of design points
for the physical experiment. Denote the responses of the n physical experimental runs
by y1p , . . . , ynp respectively, with p standing for “physical”. Let Θ be the domain of the
calibration parameter. In this article, we suppose the computer model is deterministic,
i.e., the computer output is a deterministic function of the control variables and the
calibration parameters, denoted by y s (x, θ) for x ∈ Ω, θ ∈ Θ with s standing for
“simulation”.
We consider two types of computer models. The first is called “cheap computer
simulations”. In these problems each run of the computer code takes only a short time
so that we can call the computer simulation code inside our statistical analysis program
which is usually based on an iterative algorithm like the Markov Chain Monte Carlo
(MCMC). The second is called “expensive computer simulations”. In these problems
each run of the computer code takes a long time so that it is unrealistic to embed
the computer simulation code into an iterative algorithm. A standard approach in
computer experiments is to run the computer code over a set of selected points, and
build a surrogate model based on the obtained computer outputs to approximate the
underlying true function. The surrogate model can be evaluated much faster. In the
statistical analysis, the response values from the surrogate model are used instead of
those from the original computer model.

2.1

The Case of Cheap Computer Simulations

We model the physical response y p in the following nonparametric manner
yip = ζ(xi ) + ei ,

(1)

where ζ(·) is an underlying function, referred to as the true process, and ei ’s are the
observation error. We assume ei ’s are independent and identically distributed normal
random variables with mean zero and unknown variance σ 2 . The computer output
function and the physical true process are linked by
ζ(·) = y s (·, θ0 ) + δ(·),

(2)

where θ0 denotes the “true” calibration parameter (from a physical point of view),
and δ denotes an underlying discrepancy function between the physical process and
the computer model under the true calibration parameters. It is reasonable to believe
that in most computer experiment problems, the discrepancy function δ should be
nonzero and possibly highly nonlinear because the computer codes are usually built
under assumptions or simplifications that do not hold true in reality.

3

To estimate θ0 and δ, we follow a standard Bayesian procedure by imposing certain
prior distributions on the unknown parameters θ0 and σ 2 and the unknown function
δ(·). In the computer experiment literature, a prominent method is to use a Gaussian
process as the prior for an unknown function (Santner et al., 2003). There are two
major reasons for choosing Gaussian processes. First, the sample paths of a Gaussian
process are smooth if a smooth covariance function is chosen, which can be beneficial when the target function is smooth as well. Second, the computational burden
of the statistical inference and prediction for a Gaussian process model is relatively
low. Specifically, we use a Gaussian process with mean zero and covariance function
τ 2 Cγ (·, ·) as the prior of δ(·), where Cγ is a stationary kernel with hyper-parameter γ.
In view of the finite-dimensional distribution of a Gaussian process, given τ 2 and
γ, δ(x) = (δ(x1 ), . . . , δ(xn ))T follows the multivariate normal distribution N (0, τ 2 Σγ ),
where Σγ = (Cγ (xi , xj ))ij . In order to discuss the prediction problem later, we apply the data augmentation algorithm of Tanner and Wong (1987) and consider the
posterior distribution of (θ0 , δ(x), σ 2 , γ) given by
π(θ0 , δ(x), τ 2 , σ 2 , γ|yp )
∝ π(yp |θ0 , δ(x), τ 2 , σ 2 , γ)π(δ(x)|θ0 , τ 2 , σ 2 , γ)π(θ0 , τ 2 , σ 2 , γ)


1
p
s
2
−n/2
∝ σ
exp − 2 ky − y (x, θ0 ) − δ(x)k
2σ
)
(
δ(x)T Σ−1
γ δ(x)
−1/2
−n/2
π(θ0 , τ 2 , σ 2 , γ),
×τ
(det Σγ )
exp −
2τ 2

(3)

where yp = (y1p , . . . , ynp )T , y s (x, θ0 ) = (y s (x1 , θ0 ), . . . , y s (xn , θ0 ))T . It is not timeconsuming to evaluate the posterior density function π(·, ·, ·, ·, ·|yp ) because the computer code is cheap to run. A standard MCMC procedure can then be employed to
draw samples from the posterior distribution. We refer to Higdon et al. (2004) for
further details.
In this work, we pay special attention to the prediction for a new physical reponse
at an untried point xnew , denoted as y p (xnew ). Samples from the posterior predictive
distribution of y p (xnew ) can be drawn along with the MCMC sampling. To see this,
we note that in view of the Gaussian process assumption, given δ(x) and γ, δ(xnew )
follows the normal distribution
2
T −1
−1
N (ΣT
1 Σγ δ(x), τ (Cγ (xnew , xnew ) − Σ1 Σγ Σ1 )),

where Σ1 = (Cγ (x1 , xnew ), . . . , Cγ (xn , xnew ))T . Because in each iteration of the MCMC
procedure a sample of (δ(x), θ0 , γ, σ 2 ) is drawn, we can draw a sample of y p (xnew ) from
its posterior distribution π(y p (xnew )|yp , δ(x), θ0 , τ 2 , γ, σ 2 ), which is the multivariate
normal distribution
−1
2
T −1
2
N (y s (xnew , θ0 ) + ΣT
1 Σγ δ(x), τ (Cγ (xnew , xnew ) − Σ1 Σγ Σ1 ) + σ ).

4

(4)

2.2

The Case of Expensive Computer Simulations

When the computer code is expensive to run, it is intractable to run MCMC based
on (3) directly. Instead, we need a surrogate model to approximate the computer
output function y s (·, ·). In this setting Kennedy and O’Hagan (2001) use the Gaussian
process modeling again. Suppose we first run the computer simulation over a set of
design points {(xs1 , θ1s ), . . . , (xsl , θls )} ⊂ Ω × Θ. We choose a Gaussian process with
mean mβ (·) and covariance function τ ′2 Cγs′ (·, ·) as the prior for y s , where β, τ ′ and
γ ′ are hyper-parameters. Besides, the prior processes of y s and δ are assumed to be
independent.
The Bayesian analysis for the present model is similar to that in Section 2.1 but
with more cumbersome derivations. We write ys := (y s (xs1 , θ1s ), . . . , y s (xsl , θls ))T and
define (n + l)-dimensional vectors
T
E
:= (x1 , . . . , xn , xs1 , . . . , xsl )T ,
xE = (xE
1 , . . . , xn+l )
E
)T := (θ0 , . . . , θ0 , θ1s , . . . , θls )T .
θ E = (θ1E , . . . , θn+l

By (1) and (2), the joint distribution of yp and ys conditional on θ0 , σ 2 , γ, β and τ is



Σ11 + σ 2 In 0
p
s
2
2 ′2 ′
E
(y , y )|σ , γ, β, τ , τ , γ ∼ N mβ (x ), ΣE +
,
0
0
T
E
where mβ (xE ) = (mβ (xE
1 ), . . . , mβ (xn+l )) and

ΣE =
Σ11 =

τ ′2 Cγs′



E
E E
xE
,
i , θ i , xj , θ j
ij

τ 2 Cγ (xi , xj ) ij .

Then the posterior distribution of the parameters is given by

π(θ0 , σ 2 , γ, β, τ 2 , τ ′2 , γ ′ |yp , ys ) ∝ π(yp , ys |θ0 , σ 2 , γ, β, τ 2 , τ ′2 , γ ′ )π(θ0 , σ 2 , γ, β, τ 2 , τ ′2 , γ ′ ).
The parameter estimation proceeds in a similar manner to the MCMC scheme discussed
in Section 2.1. As before, the prediction for the true process can be done along with
the MCMC iterations. Noting the fact that (y p (xnew ), yp , ys ) follows a multivariate
normal distribution given the model parameters, the posterior predictive distribution
of y p (xnew ) can be obtained using the Bayes’ theorem.
It can be seen that the modeling and analysis for the KO method with expensive
computer code is much more complicated than that with cheap computer code. For
the ease of mathematical analysis, our theoretical studies in the next section considers
only the cases with cheap code. Hence, we omit the detailed formulae of the posterior
density of the model parameters and the posterior predictive distribution of y p (xnew )
in this section.

5

3

Theoretical Studies

In this section we conduct some theoretical study on the power of prediction of the
KO method. For the ease of the mathematical treatment, we only consider the case of
cheap computer code, because the formulae for the case of expensive computer code
are much more complicated and cumbersome as shown in Section 2.2. We believe that
this simplification does not affect our general conclusion.
The mathematical treatment to develop the asymptotic theory for the KO method
also depend on the choice of the correlation family Cγ . In the present work, we restrict
ourselves with the Matérn family of kernel functions (Stein, 1999), defined as
Cυ,γ (s, t) =

υ

√
√
1
2 υγks − tk Kυ 2 υγks − tk ,
υ−1
Γ(υ)2

(5)

where Kυ is the modified Bessel function of the second kind. In Matérn family, the
model parameter υ dominates the smoothness of the process and γ is a scale parameter.
Because the smoothness parameter υ has an effect on the rate of convergence of the
prediction, for simplicity we suppose υ is fixed in the entire data analysis.
All proofs in this section are postponed to Appendix A.

3.1

A Function Approximation Perspective

In this section we follow the theoretical framework of Tuo and Wu (2016) to study
the prediction performance of the KO method. Under this framework, the physical
responses are assumed to have no random error, i.e., ǫi ’s in (1) are zero. This is an
unrealistic assumption in practice. But this assumption simplifies the model structure,
so that we are able to find some mathematical tools which help us to understand certain
intrinsic properties of the KO method.
From (1), we have yip = ζ(xi ). We remind that ζ is indeed a deterministic function
(as the expectation of the physical response). Therefore, we will regard the Gaussian
process modeling technique used in the KO method as a way of reconstructing the
function ζ based on samples ζ(xi ).
An immediate consequence of the deterministic assumption is δ(x) = yp − y s (x, θ0 ),
i.e., δ(x) is determined by θ0 given the observations. Thus (3) is not applicable. Instead,
we have
π(θ0 , τ 2 , γ|yp ) ∝ π(yp |θ0 , τ 2 , γ)π(θ0 , τ 2 , γ)


1 p
s
T −1 p
s
−1/2
∝ (det Σγ )
exp − (y − y (x, θ0 )) Σγ (y − y (x, θ0 ) π(θ0 , τ 2 , γ).
2
To differentiate between the true process ζ and its estimate based on the observations, we denote a draw from the predictive distribution π(ζ(xnew )) by ζ rep (xnew ).
Then the posterior predictive distribution π(ζ rep (xnew )|θ0 , γ, yp ) is

−1 p
s
2
T −1
N y s (xnew , θ0 ) + ΣT
(6)
1 Σγ (y − y (x, θ0 )), τ (Cυ,γ (xnew , xnew ) − Σ1 Σγ Σ1 ) .

6

We now suppose the prior distribution π(θ0 , τ 2 , γ) is separable, i.e., π(θ0 , τ 2 , γ) =
π(θ0 )π(τ 2 )π(γ). Let Sθ , Sτ 2 and Sγ denote the supports of the distributions π(θ0 ), π(τ 2 )
and π(γ) respectively. For the ease of mathematical treatment, we further suppose that
Sθ is a compact subset of R, and Sτ 2 ⊂ [0, τ02 ], Sγ ⊂ [γ1 , γ2 ] for some 0 < τ02 < +∞, 0 <
γ1 < γ2 < +∞. The independence assumption of the prior distributions can be replaced
with a more general assumption, which would not affect the validity of our theoretical
analysis. However, the compact support assumption is technically unavoidable in the
current treatment. Because here we only focus on the posterior mode, the use of the
compact support assumption does not affect the practical applicability of the results.
The aim of this section is to study the asymptotic behavior of
−1 p
s
µ̂θ,γ = y s (xnew , θ) + ΣT
1 Σγ (y − y (x, θ)),

−1
ςˆτ22 ,γ = τ 2 (Cυ,γ (xnew , xnew ) − ΣT
1 Σγ Σ1 ),

as the design points become dense in Ω, for (θ, τ 2 , γ) ∈ Sθ , Sτ 2 , Sγ . Clearly, the true
posterior mean of ζ rep (xnew ) given by (6) is
E[ζ rep (xnew )|yp ] = E[µ̂θ̂,γ̂ |yp ],
where (θ̂, γ̂) follows the posterior distribution π(θ0 , γ|yp ). Note that
|E[ζ rep (xnew )|yp ] − ζ(xnew )|
o
n
= E E[ζ rep (xnew ) − ζ(xnew )|yp , θ̂, γ̂] yp
≤

=

sup

θ∈Sθ ,γ∈Sγ

sup
θ∈Sθ ,γ∈Sγ

|E[ζ rep (xnew ) − ζ(xnew )|yp , θ, γ]|

|µ̂θ,γ − ζ(xnew )|,

i.e., the bias of the posterior predictive mean can be bounded by the supremum of
|µ̂θ,γ − ζ(xnew )|. Similarly, we find
Var(ζ rep (xnew )|yp ) ≤

sup
τ 2 ∈S

τ 2 ,γ∈Sγ

ςˆτ22 ,γ .

In this section we will bound supθ∈Sθ ,γ∈Sγ |µ̂θ,γ − ζ(xnew )| and supτ 2 ∈Sτ 2 γ∈Sγ ςˆτ22 ,γ .
To this end, we resort to the theory of native spaces. We refer to Wendland (2005)
for detailed discussions. For a symmetric and positive definite function Φ over Ω × Ω,
consider the linear space
)
(m
X
αi Φ(si , ·) : m ∈ N+ , αi ∈ R ,
FΦ (Ω) :=
i=1

equipped with the inner product
+
*m
l
m X
l
X
X
X
αi βj Φ(si , tj ).
βj Φ(tj , ·) =
αi Φ(si , ·),
i=1

i=1 j=1

j=1

7

(7)

The completion of FΦ (Ω) with respect to its inner product is called the native space
generated by Φ, denoted by NΦ (Ω). Denote the inner product and the norm of NΦ (Ω)
by h·, ·iNΦ (Ω) and k · kNΦ (Ω) respectively.
Now we state the interpolation scheme in the native space. Let f ∈ NΦ (Ω) and
x = {x1 , . . . , xn } a set of distinct points in Ω. Let y = (f (x1 ), . . . , f (xn ))T be the
observerd data. Define
sf,x (x) =

n
X

ui Φ(xi , x),

(8)

i=1

where u = (u1 , . . . , un )T is given by the linear equation
y = Φ(x, x)u
for (Φ(x, x))ij = Φ(xi , xj ).
Clearly, sf,x ∈ FΦ and thus sf,x ∈ NΦ (Ω). The next lemma can be found in
Wendland (2005). For the completeness of the present article, we provide its proof in
Appendix A.
Lemma 1. For f ∈ NΦ (Ω) and a set of design points x ⊂ Ω,
hsf,x , f − sf,x iNΦ (Ω) = 0.
From Lemma 1 we can deduce the Pythagorean identity
ksf,x k2NΦ (Ω) + kf − sf,x k2NΦ (Ω) = kf k2NΦ (Ω) .

(9)

Now we consider an arbitrary function h ∈ NΦ (Ω) which interpolates f over x, denoted
as f |x = h|x . Then we have sf,x = sh,x and thus (9) also holds true if we replace f
with h. This suggests ksf,x kNΦ (Ω) ≤ khkNΦ (Ω) , which yields the following optimality
condition
sf,x = argmin khkNΦ (Ω) ,

(10)

h∈NΦ (Ω)
h|x =f |x

i.e., sf,x has the minimum native norm among all functions in NΦ (Ω) that interpolate
f over x.
It can be shown that the native space generated by the Matérn kernel Cυ,γ for υ ≥ 1
coincides with the (fractional) Sobolev space H υ+d/2 (Ω) (Adams and Fournier, 2003),
and the norms are equivalent. See Tuo and Wu (2016) for details. Moreover, we can
also prove that the norms of the native spaces generated by Cυ,γ for a set of γ values
bounded away from 0 and +∞ are equivalent.
Lemma 2. Suppose υ ≥ 1. There exist constants c1 , c2 > 1, so that
c1 kf kH υ+d/2 (Ω) ≤ kf kNCυ,γ (Ω) ≤ c2 kf kH υ+d/2 (Ω)
holds for all f ∈ H υ+d/2 (Ω) and all γ ∈ [γ1 , γ2 ].

8

(11)

Next, we turn to the error estimate of the interpolant sf,x . Wendland (2005) shows
that for u ∈ H µ (Ω) with u|x = 0 and ⌊µ⌋ > d/2,
µ−d/2

kukL∞ (Ω) ≤ Chx,Ω

kukH µ (Ω) ,

provided that x is “sufficiently dense”, where C is independent of x and u; hx,Ω is the
fill distance of the design x defined as
hx,Ω = sup min kx − xj k.
x∈Ω xj ∈x

Here “x is sufficiently dense” means that its fill distance hx,Ω is less than a constant
h0 depending only on Ω and µ. Noting the fact that (f − sf,x )|x = 0 and f − sf,x ∈
H υ+d/2 (Ω), we obtain that for υ ≥ 1,
kf − sf,x kL∞ (Ω) ≤ Chυx,Ω kf − sf,x kH υ+d/2 (Ω) ,
which, together with (9), yields
kf − sf,x kL∞ (Ω) ≤ Chυx,Ω kf kH υ+d/2 (Ω) .

(12)

Then we apply Lemma 2 to prove Lemma 3.
Lemma 3. Suppose υ ≥ 1. For f ∈ H υ+d/2 (Ω), let sf,x be the interpolant of f over x
with the kernel Cγ,υ , γ ∈ [γ1 , γ2 ]. Then for sufficiently dense x
kf − sf,x kL∞ (Ω) ≤ Chυx,Ω kf kNCυ,γ (Ω) ,
where C is independent of the choices of f , x and γ.
Following the notation of Tuo and Wu (2016), we define ǫ(x, θ) = ζ(x) − y s (x, θ).
It is commented by Tuo and Wu (2016) that in general θ0 is not estimable due to the
identifiability problem, and thus neither is δ(·) = ǫ(·, θ0 ). However, as will be shown
later, the function ǫ(·, ·) can be consistently estimated using KO calibration. Suppose
ǫ(·, θ) ∈ H υ+d/2 (Ω) for each θ ∈ Sθ . Let ǫ(x, θ) = (ǫ(x1 , θ), . . . , ǫ(xn , θ))T . Clearly,
yp − y s (x, θ) = ǫ(x, θ) and thus
−1 p
s
sǫ(·,θ),x (xnew ) = ΣT
1 Σγ (y − y (x, θ)).

By (12) we obtain
|µ̂θ,γ − ζ(xnew )| = |ǫ(xnew , θ) − sǫ(·,θ),x (xnew )|
≤ Chυx,Ω kǫ(·, θ)kH υ+d/2 (Ω)

≤ Chυx,Ω sup kǫ(·, θ)kH υ+d/2 (Ω) .

(13)

θ∈Sθ

The error bound for the variance term can be obtained similarly. Elementary calculations show that
−1
ΣT
1 Σγ Σ1 = sCυ,γ (·,xnew ),x (xnew ).

9

Hence we apply Lemma 3 to find
−1
2
|τ 2 (Cυ,γ (xnew , xnew ) − ΣT
1 Σγ Σ1 )| = τ |Cυ,γ (xnew , xnew ) − sCυ,γ (·,xnew ),x (xnew )|

≤ τ02 Chυx,Ω kCυ,γ (·, xnew )kNCυ,γ (Ω)
= τ02 Chυx,Ω ,

(14)

where the last equality follows from the fact that kCυ,γ (·, xnew )kNCυ,γ (Ω) = 1. We
summarize our findings in (13) and (14) as Theorem 1.
Theorem 1. Suppose υ ≥ 1, γ ∈ [γ1 , γ2 ], τ ≤ τ0 . Then for a sufficiently dense design
x, we have the upper bound for the predictive mean as
sup
θ∈Sθ ,γ∈Sγ

|µ̂θ,γ − ζ(xnew )| ≤ Chυx,Ω sup kǫ(·, θ)kH υ+d/2 (Ω) ,
θ∈Sθ

and the upper bound for the predictive variance as
sup
τ 2 ∈Sτ2 ,γ∈Sγ

ςˆτ22 ,γ ≤ τ02 Chυx,Ω ,

with a constant C depending only on Ω, υ, γ1 , γ2 .
From Theorem 1, the rate of convergence is O(hυx,Ω ), which is known to be optimal in the current setting (Wendland, 2005). It is worth noting that the predictive
behavior of the KO calibration is more robust than in the case of estimation as shown
by Tuo and Wu (2016), Theorem 4.2. Specifically, they show the KO calibration estimator tends to the minimizer of a norm involving the prior assumption, i.e., the KO
calibration can reply heavily on the prior specification. By comparison, the predictive
performance as shown in Theorem 1 above does not depend on the choice of the prior
asymptotically.

3.2

A Nonparametric Regression Perspective

Now we turn to a more realistic case, where the physical observations have random
measurement errors. As before, we treat the true process ζ(·) as a deterministic function. For the ease of mathematical treatment, in this section we fix the value of γ. Our
analysis later will show in Theorems 3 and 4 that the resulting rate of convergence is
not influenced by the choice of γ. Other parameters are either estimated or chosen to
vary along with the sample size n.
To study the predictive behavior of the KO method asymptotically, the key is
to understand the posterior mode of δ(x) in (3). We first introduce the representer
theorem (Schölkopf et al., 2001; Wahba, 1990). We also give a proof of the representer
theorem using Lemma 1 in Appendix A.

10

Lemma 4 (Representer Theorem). Let x1 , . . . , xn be a set of distinct points in Ω and
L : Rn → R be an arbitrary function. Denote the minimizer of the optimization
problem
min L(f (x1 ), f (x2 ), . . . , f (xn )) + kf k2NΦ (Ω)

f ∈NΦ (Ω)

by fˆ. Then fˆ possesses the representation
fˆ =

n
X
i=1

αi Φ(xi , ·),

with coefficients αi ∈ R, i = 1, . . . , n.
Similar to Section 3.1, we first fix the values of τ 2 , σ 2 , γ in their domain. Then
we consider the profile posterior density function of δ(x), which, according to (3), is
proportional to
(
)
T Σ−1 δ(x)
δ(x)
1
γ
πτ 2 ,σ2 ,γ (θ, δ(x)) = exp − 2 ky2 − y s (x, θ) − δ(x)k2 −
.
(15)
2σ
2τ 2
The profile posterior mode (θ̂KO , δ̂(x)) maximizes πτ 2 ,σ2 ,γ (·, ·). Using the representer theorem, we show an equality between δ̂(x) and the solution to a penalized least
squares problem.
ˆ be the solution to
Theorem 2. Let (θ̂, ∆)
argmin

n
X
σ2
(yip − y s (xi , θ) − f (xi ))2 + 2 kf k2NC (Ω) .
υ,γ
τ

(16)

θ∈Θ
f ∈NCυ,γ (Ω) i=1

ˆ 1 ), . . . , ∆(x
ˆ n ))T =: ∆(x)
ˆ
Then θ̂ = θ̂KO and (∆(x
= δ̂(x).
Now we are ready to state the main asymptotic theory. We will first investigate the
asymptotic properties of the predictive mean. Next we will consider the consistency of
the predictive variance.
From (4), the predictive mean of the KO model is
−1
ζ rep (xnew ) = y s (xnew , θ̂KO ) + ΣT
1 (xnew )Σγ δ̂(x),

(17)

where Σ1 (xnew ) = (Cυ,γ (xnew , x1 ), . . . , Cυ,γ (xnew , xn ))T . Invoking Theorem 2, we have
ˆ
ˆ defined in (16). Using the formula of the kernel interpolant
δ̂(x) = ∆(x)
with ∆
given by (8), it can be seen that ζ̂(·) − y s (·, θ̂KO ) is the kernel interpolant of the data
ˆ
(x, ∆(x)).
Hence, from Lemma 4 and Theorem 2 we have
ˆ
ζ̂(·) − y s (·, θ̂KO ) = ∆(·).

11

We note that the ratio of the variances σ 2 /τ 2 plays an important role in (16). In
the nonparametric regression literature, such a quantity is commonly referred to as
the smoothing parameter. The smoothing parameter is a tuning parameter to balance
the bias and variance of the estimator. It can be seen that as σ 2 /τ 2 → ∞, ǫ̂ tends
to 0, which has the smallest variance but a large bias; as σ 2 /τ 2 ↓ 0, ǫ̂ will eventually
interpolate (xi , yip − y s (xi , θ̂KO )), which typically leads to an over-fitting problem. We
denote σ 2 /τ 2 by rn when the sample size is n. According to Theorem 3, the optimal
4υ+5d

rate for rn is rn ∼ n 4υ+4d . There is a theory by van der Vaart and van Zanten (2008)
which says that the optimal tuning rate can be automatically achieved by following a
standard Bayesian analysis procedure. To save space, we do not pursue this approach
here.
Some asymptotic theory for the penalized least squares problem (16) is available
in the literature (van der Geer, 2000). In order to employ such a theory, we need to
choose the smoothing parameter rn to diverge at an appropriate rate as n goes to
infinity. For convenience, we suppose that the design points are randomly chosen.
We consider the rate of convergence of the penalized least squares estimator under
the L2 metric. We assume that y s is Lipschitz continuous. Then the metric entropy
of {y s (·, θ) : θ ∈ Θ} is dominated by that of the unit ball of the nonparametric class
NCυ,γ (Ω). See van der Vaart and Wellner (1996) for details. Theorem 3 then is a direct
consequence of Theorem 10.2 of van der Geer (2000), where the required upper bound
for the metric entropy is obtained from (3.6) of Tuo and Wu (2015).
Theorem 3. Suppose the design points {xi } are independent samples from the uniform
distribution over Ω. We assume υ ≥ 1 and y s is Lipschitz continuous. Choose rn
4υ+5d
appropriately so that rn ∼ n 4υ+4d . Under model (1) with σ 2 > 0, the KO predictor ζ̂
defined in (17) has the approximation properties
n

2υ+d
1X
(ζ̂(xi ) − ζ(xi ))2 = Op (n− 2υ+2d ),
n

(18)

i=1

υ+d/2

kζ̂ − ζkL2 (Ω) = Op (n− 2υ+2d ),

(19)

and
kζ̂ − ζkH υ+d/2 (Ω) = Op (1).

(20)

Since the native space NCυ,γ (Ω) is equivalent to the Sobolev space H υ+d/2 (Ω), the
rate of convergence in (19) is optimal according to the theory by Stone (1982). The interpolation inequality of the Sobolev spaces (see (12) in Chapter 5 of Adams and Fournier
(2003)) claims that
d

1−

d

2(υ+d/2)
2(υ+d/2)
kζ̂ − ζkL∞ (Ω) ≤ Kkζ̂ − ζkH
,
υ+d/2 (Ω) kζ̂ − ζkL2 (Ω)

with constant K depending only on Ω and υ. In view of (19) and (20), we have
υ

kζ̂ − ζkL∞ (Ω) = Op (n− 2υ+2d ),

12

(21)

which gives the rate of convergence of the predictive mean under the uniform metric.
Now we turn to the consistency of the predictive variance. To avoid ambiguity, we
denote the true value of σ 2 by σ02 . First we consider the profile posterior mode of σ 2 in
(3). For simplicity, we only consider the non-informative prior for σ 2 with π(σ 2 ) ∝ 1.
However, we note that the limiting value of the posterior mode of σ 2 is not affected by
the choice of π(σ 2 ), provided that σ02 is contained in the support of π(σ 2 ). It is easily
seen from (3) that the posterior mode of σ 2 is
kyp − y s (x, θ̂KO ) − δ̂(x)k2
n
n
X
1
{ei + (ζ̂(xi ) − ζ(xi ))}2
n

σ̂ 2 =
=

i=1

n
n
n
1X 2 2X
1X
ei +
ei (ζ̂(xi ) − ζ(xi )) +
(ζ̂(xi ) − ζ(xi ))2 ,
n
n
n

=

i=1

i=1

i=1

which yields the following asymptotic property
n

|σ̂ 2 − σ02 | =

σ̂ 2 −

1X 2
ei + Op (n−1/2 )
n
i=1

n
n
2X
1X
=
ei (ζ̂(xi ) − ζ(xi )) +
(ζ̂(xi ) − ζ(xi ))2 + Op (n−1/2 )
n
n
i=1
i=1
!1/2
!1/2
n
n
1X 2
1X
≤ 2
ei
(ζ̂(xi ) − ζ(xi ))2
n
n
i=1

i=1

n
1X
+
(ζ̂(xi ) − ζ(xi ))2 + Op (n−1/2 )
n
i=1

υ+d/2

= Op (n− 2υ+2d ),

(22)

where the inequality follows from Cauchy-Schwarz inequality and the last equality
follows from (18).
From (4), the predictive variance of the KO model is
−1
2
ςˆ2 (xnew ) = τ 2 (Cγ (xnew , xnew ) − ΣT
1 Σγ Σ1 ) + σ̂ .

(23)

−1
As discussed in Section 3.1, Cγ (xnew , xnew ) − ΣT
1 Σγ Σ1 is the approximation error
of the kernel interpolation for the function Cγ (·, xnew ). Clearly, the error from the
interpolation problem discussed in Section 3.1 should be no more than that from the
smoothing problem discussed in the current section, because of the presence of the

13

random error in the latter situation. Thus we have
−1
sup |τ 2 (Cγ (xnew , xnew ) − ΣT
1 Σγ Σ1 )|

xnew ∈Ω

−1
sup |rn σ̂ 2 (Cγ (xnew , xnew ) − ΣT
1 Σγ Σ1 )|

=

xnew ∈Ω

d

υ

υ+d/2

= Op (n− 4υ+4d n− 2υ+2d ) = Op (n− 2υ+2d ),

(24)
d

where the second equality follows from the assumption rn ∼ n 4υ+4d in Theorem 3, (23)
and (21). Combining (23) and (24) we obtain Theorem 4.
Theorem 4. Under the conditions of Theorem 3, we have the error bound for the
predictive variance under the uniform metric
υ+d/2

kˆ
ς 2 (·) − σ02 kL∞ (Ω) = Op (n− 2υ+2d ).
Noting that σ02 is the variance of the random noise, which is present in prediction
for a new physical response. In other word, no predictor has a mean square error less
than σ02 . Theorems 3 and 4 reveal that the predictive distribution given by the KO
method can capture the true uncertainty of the physical data in the asymptotic sense.

4

Discussions

In this work, we prove some error bounds for the predictive error given by the KennedyO’Hagan method in two cases: 1) the physical observations have no random error
and 2) the physical observations are noisy. For the ease of mathematical analysis,
we only consider the Matérn correlation family. If a different covariance structure
is used, we believe that the consistency for the predictive mean and the predictive
variance still holds. However, additional study is required to obtain the appropriate
rate of convergence. In our entire analysis, we ignore the estimation for some model
parameters like γ and τ 2 . One may consider the error estimate in a fully Bayesian
procedure. But the analysis will then become rather complicated and it is unclear
whether a new theory can be developed along the same lines.
Throughout this work, we assume that the smoothness parameter υ is given. From
Theorems 1 and 3, a better rate of convergence can be obtained by using a greater υ
provided that the target function still lies in NCυ,γ (Ω). Thus ideally, one should choose
υ to be close to, but no more than the true degree of smoothness of the target function.
There are different ways of choosing data-dependent υ, but the mathematical analysis
will become much more involved. We refer to Loh et al. (2015) and the references
therein for some related discussions.
In Section 3.2, we assume that the design points xi ’s are random samples over Ω.
In practice, one may also wish to choose design points using a systematic (deterministic) scheme. In general, if a sequence of fixed designs is used, the same (optimal)
rate of convergence is retained, provided that these designs satisfy certain space-filling

14

conditions. We refer to Utreras (1988) for the results and necessary mathematical
tools.
Finally, we discuss how the calibration procedure can have an effect on the prediction for the true process. In this article, we allow the number of physical measurements
to grow to infinity and obtain the rate of convergence. By comparing the results presented here and the standard ones using radial basis functions or smoothing spline
approximation, we find that the rate of convergence is not elevated by doing calibration. But we can use the following heuristics to show that by doing KO calibration the
predictive error can be improved by a constant factor. To see this, we review the proof
of Theorem 1, from which it can be seen that if we fix Φ and γ, the predictive error is
bounded by
|µ̂θ,γ − ζ(xnew )| ≤ Chυx,Ω kǫ(·, θ)kH υ+d/2 (Ω) ,

(25)

for an arbitrarily chosen θ ∈ Θ. So the rate of convergence is given by O(hυx,Ω ), and
kǫ(·, θ)kH υ+d/2 (Ω) acts as a constant factor. Tuo and Wu (2016) show that under certain
conditions, the KO estimator for the calibration parameter converges to
θ ′ = argmin kǫ(·, θ)kNΦ (Ω) ,
θ∈Θ

as the design points become dense over Ω. Since k·kNΦ (Ω) is equivalent to k·kH υ+d/2 (Ω) ,
estimating the calibration parameter via the KO method is apparently beneficial for
prediction in the sense that the upper error bound is reduced because kǫ(·, θ ′ )kNΦ (Ω) ≤
kǫ(·, θ)kNΦ (Ω) for all θ ∈ Θ. There is a similar phenomenon for the stochastic case, by
using the arguments in the proof of Theorem 10.2 of van der Geer (2000).

Appendix
A

Technical Proofs

In this section, we present the proofs for Lemma 1, Lemma 2, Lemma 4 and Theorem
2.
Proof of Lemma 1. We first assume f ∈ FΦ . If f = sf,x , there is nothing to prove. If
f 6= sf,x , without loss of generality, we write
f (x) =

n+m
X

αi Φ(x, xi ),

i=1

for an extra set of distinct points {xn+1 , . . . , xn+m } ⊂ Ω. Now partition (Ai,j ) =
Φ(xi , xj ), 1 ≤ i, j ≤ n + m into


(A1 )n×n (A2 )n×m
,
A=
(A3 )m×n (A4 )m×m

15

where A3 = AT
2 because Φ is symmetric.
Let y = (f (x1 ), . . . , f (xn ))T , a1 = (α1 , . . . , αn )T , a2 = (αn+1 , . . . , αn+m )T . Clearly,
y = A1 a1 + A2 a2 . By the definition of sf,x , we have
sf,x (x) =

n
X

ui Φ(x, xi ),

i=1

with u = (u1 , . . . , un )T satisfying y = A1 u. Then from (7) we obtain
hsf,x , f − sf,x iNΦ (Ω)
+
* n
n+m
n
X
X
X
αi Φ(x, xi )
(αi − ui )Φ(x, xi ) +
ui Φ(x, xi ),
=


 A1 A2
0
A3 A4

uT

=

i=n+1

i=1

i=1



NΦ (Ω)


a1 − u
a2

= uT (A1 a1 + A2 a2 − A1 u)

= uT (y − y) = 0.

(26)

For a general f ∈ NΦ (Ω), we can find a sequence fn ∈ FΦ with fn → f in NΦ (Ω)
as n → ∞. The desired result then follows from a limiting form of (26).
Proof of Lemma 2. For any g ∈ L2 (Rd ) ∩ C(Rd ), its native norm admits the representation
Z
|g̃(ω)|2
2
−d/2
dω,
(27)
kgkNΦ (Rd ) = (2π)
Rd Φ̃(ω)
where g̃ and Φ̃ denote the Fourier transforms of g and Φ respectively. See Theorem
10.12 of Wendland (2005). The (fractional) Sobolev norms have a similar representation
Z
kgk2H s (Rd ) = (2π)−d/2
|g̃(ω)|2 (1 + kωk2 )s dω.
(28)
Rd

See Adams and Fournier (2003) for details. Tuo and Wu (2016) show that
C̃υ,γ (ω) = 2d/2 (4υγ 2 )υ

Γ(υ + d/2)
(4υγ 2 + kωk2 )−(υ+d/2) .
Γ(ν)

Using the inequality
(1 + b) min(1, a) ≤ a + b ≤ (1 + b) max(1, a),
for a, b ≥ 0, we obtain
C̃υ,γ (ω)

≤
≤

o
n
Γ(υ + d/2)
max 1, (4υγ 2 )−(υ+d/2) (1 + kωk2 )−(υ+d/2)
Γ(υ)
n
o
Γ(υ + d/2)
max (4υγ22 )υ , (4υγ12 )−d/2 (1 + kωk2 )−(υ+d/2)
2d/2
Γ(υ)

2d/2 (4υγ 2 )υ

=: C1 (1 + kωk2 )−(υ+d/2) ,

(29)

16

and
C̃υ,γ (ω)

≥
≥

o
n
Γ(υ + d/2)
min 1, (4υγ 2 )−(υ+d/2) (1 + kωk2 )−(υ+d/2)
Γ(υ)
o
n
Γ(υ + d/2)
2d/2
min (4υγ12 )υ , (4υγ22 )−d/2 (1 + kωk2 )−(υ+d/2)
Γ(υ)

2d/2 (4υγ 2 )υ

=: C2 (1 + kωk2 )−(υ+d/2) ,

(30)

hold for all ω ∈ Rd .
Now we apply the extension theorem of the native spaces (Theorem 10.46 of Wendland,
2005) to obtain a function f E ∈ NCυ,γ (Rd ) such that f E |Ω = f and kf kNCυ,γ (Ω) =
kf E kNCυ,γ (Rd ) for each γ ∈ [γ1 , γ2 ]. We use (27)-(29) to obtain
kf k2NC

υ,γ (Ω)

= kf E k2NC

υ,γ (R

= (2π)−d/2

d)

≥ C1−1 (2π)−d/2

Z

Rd

Z

Rd

|f˜E (ω)|2
dω
C̃υ,γ (ω)

|f˜E (ω)|2 (1 + kωk2 )υ+d/2 dω

= C1−1 kf E k2H υ+d/2 (Rd ) ≥ C1−1 kf k2H υ+d/2 (Ω) ,

(31)

where the last inequality follows from the fact that f E |Ω = f . On the other hand,
because Ω is convex, f has an extension fE ∈ H υ+d/2 (Rd ) satisfying kfE kH k (Rd ) ≤
ckf kH k (Ω) for some constant c independent of f . Then we use (27), (28) and (30) to
obtain
kfE k2H k (Ω) ≥ c−2 kf k2H k (Ω)
−2

= c

−2

≥ c

−2

= c

−d/2

(2π)

Z

Rd

−d/2

C2 (2π)

Z

|f˜E (ω)|2 (1 + kωk2 )υ+d/2 dω
|f˜E (ω)|2
dω
C̃υ,γ (ω)

Rd
2
C2 kfE kNC (Rd )
υ,γ

≥ c−2 C2 kf k2NCυ,γ (Ω) ,

where the last inequality follows from the restriction theorem of the native space, which
states that the restriction f = fE |Ω is contained in NCυ,γ (Ω) with a norm that is less
than or equal to the norm kfE kNCυ,γ (Rd ) . See Theorem 10.47 of Wendland (2005). The
desired result is proved by combining (31) and (32).
Proof of Lemma 4. For f ∈ NΦ (Ω), define
M (f ) = L(f (x1 ), . . . , f (xn )) + kf k2NΦ (Ω) .
Now consider sfˆ,X , i.e., the interpolant of fˆ over X = {x1 , . . . , xn } using the kernel
function Φ. Because fˆ(xi ) = s ˆ (xi ) for i = 1, . . . , n, we have
f ,X

L(f (x1 ), . . . , f (xn )) = L(sfˆ,X (x1 ), . . . , sfˆ,X (xn )).

17

(32)

In addition, it is easily seen from Lemma 1, (9) and (10) that
ksfˆ,X k2NΦ (Ω) ≤ kfˆk2NΦ (Ω) ,

(33)

and the equality holds if and only if sfˆ,X = fˆ. By combining (32) and (33) we obtain
M (sfˆ,X ) ≤ M (fˆ).

(34)

= fˆ,
Because fˆ minimizes M (f ), the reverse of (34) also holds. Hence we deduce sf,X
ˆ
which proves the theorem according to the definition of the interpolant.
Proof of Theorem 2. We first rewrite the minimization problem (16) as the following
iterated form
n
X
σ2
(yip − y s (xi , θ) − f (xi ))2 + 2 kf k2NC (Ω)
min
υ,γ
θ∈Θ
τ
f ∈NCυ,γ (Ω) i=1

= min

min

θ∈Θ f ∈NCυ,γ (Ω)

n
X
i=1

(yip − y s (xi , θ) − f (xi ))2 +

σ2
kf k2NC (Ω)
υ,γ
τ2

(35)

Now we apply Lemma 4 to the inner minimization problem in (35) and obtain the
ˆ
following representation for ∆:
ˆ =
∆

n
X
i=1

αi Cυ,γ (xi , ·),

with an undetermined vector of coefficients α = (α1 , . . . , αn )T . Using the definition
Σγ = (Cυ,γ (xi , xj ))ij , clearly we have the matrix representation
ˆ
∆(x)
= Σγ α.

(36)

Now using (7) we have
+
* n
n
X
X
2
ˆ
αi Cυ,γ (xi , ·)
αi Cυ,γ (xi , ·),
=
k∆k

= αT Σγ α.

NCυ,γ (Ω)

i=1

i=1

NCυ,γ (Ω)

The minimization problem (16) then reduces to
argmin kyp − y s (x, θ) − αΣγ k2 +
θ∈Θ
α∈Rn

σ2 T
α Σγ α.
τ2

Applying a change-of-variable argument using (36) we obtain the following optimization
formula
argmin kyp − y s (x, θ) − ∆(x)k2 +

θ∈Θ
∆(x)∈Rn

σ2
∆(x)T Σ−1
γ ∆(x).
τ2

Elementary calculations show its equivalence to the definition of (θ̂KO , δ̂(x)).

18

Acknowledgements
Tuo’s work is supported by the National Center for Mathematics and Interdisciplinary
Sciences in CAS and NSFC grant 11501551, 11271355 and 11671386. Wu’s work is
supported by NSF grant DMS 1564438. The authors are grateful to the Associate
Editor and referees for very helpful comments.

References
Adams, R. A. and J. J. Fournier (2003). Sobolev Spaces, Volume 140. Access Online
via Elsevier.
Higdon, D., M. Kennedy, J. Cavendish, J. Cafeo, and R. Ryne (2004). Combining
field data and computer simulations for calibration and prediction. SIAM Journal
of Scientific Computing 26, 448–466.
Kennedy, M. and A. O’Hagan (2001). Bayesian calibration of computer models. Journal
of the Royal Statistical Society: Series B 63 (3), 425–464.
Loh, W.-L. et al. (2015). Estimating the smoothness of a gaussian random field from
irregularly spaced data via higher-order quadratic variations. The Annals of Statistics 43 (6), 2766–2794.
Santner, T., B. Williams, and W. Notz (2003). The Design and Analysis of Computer
Experiments. Springer Verlag.
Schölkopf, B., R. Herbrich, and A. J. Smola (2001). A generalized representer theorem.
In Computational Learning Theory, pp. 416–426. Springer.
Stein, M. (1999). Interpolation of Spatial Data: Some Theory for Kriging. Springer
Verlag.
Stone, C. J. (1982). Optimal global rates of convergence for nonparametric regression.
The Annals of Statistics 10 (4), 1040–1053.
Tanner, M. A. and W. H. Wong (1987). The calculation of posterior distributions
by data augmentation. Journal of the American Statistical Association 82 (398),
528–540.
Tuo, R. and C. F. J. Wu (2015). Efficient calibration for imperfect computer models.
The Annals of Statistics 43, 2331–2352.
Tuo, R. and C. F. J. Wu (2016). A theoretical framework for calibration in computer models: parametrization, estimation and convergence properties. SIAM/ASA
Journal on Uncertainty Quantification 4 (1), 767–795.

19

Utreras, F. I. (1988). Convergence rates for multivariate smoothing spline functions.
Journal of Approximation Theory 52 (1), 1–27.
van der Geer, S. A. (2000). Empirical Processes in M-estimation, Volume 6. Cambridge
university press.
van der Vaart, A. W. and J. H. van Zanten (2008). Rates of contraction of posterior
distributions based on gaussian process priors. The Annals of Statistics 36 (3), 1435–
1463.
van der Vaart, A. W. and J. A. Wellner (1996). Weak Convergence and Empirical
Processes: With Applications to Statistics. Springer.
Wahba, G. (1990). Spline Models for Observational Data, Volume 59. Society for
Industrial Mathematics.
Wendland, H. (2005). Scattered Data Approximation. Cambridge University Press.
Wu, C. F. J. and M. S. Hamada (2011). Experiments: Planning, Analysis, and Optimization, Volume 552. John Wiley & Sons.

20

