Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2879656

Formats

Article sections

- Abstract
- 1. INTRODUCTION
- 2. CURRENT METHODS
- 3. A HIERARCHICAL BAYESIAN MODEL FOR TUNING AND CALIBRATION
- 4. METHODOLOGY FOR SIMULTANEOUS TUNING AND CALIBRATION
- 5. EXAMPLES AND COMPARISON
- 6. SUMMARY AND DISCUSSION
- Supplementary Material
- REFERENCES

Authors

Related links

Technometrics. Author manuscript; available in PMC 2010 June 2.

Published in final edited form as:

Technometrics. 2009 November 1; 51(4): 464–474.

doi: 10.1198/TECH.2009.08126PMCID: PMC2879656

NIHMSID: NIHMS160446

Gang Han, H. Lee Moffitt Cancer Center & Research Institute, MRC/BIOSTAT, 12902 Magnolia Drive, Tampa, FL 33612, (Email: gro.ttiffom@nah.gnag)

Tuning and calibration are processes for improving the representativeness of a computer simulation code to a physical phenomenon. This article introduces a statistical methodology for simultaneously determining tuning and calibration parameters in settings where data are available from a computer code and the associated physical experiment. Tuning parameters are set by minimizing a discrepancy measure while the distribution of the calibration parameters are determined based on a hierarchical Bayesian model. The proposed Bayesian model views the output as a realization of a Gaussian stochastic process with hyperpriors. Draws from the resulting posterior distribution are obtained by the Markov chain Monte Carlo simulation. Our methodology is compared with an alternative approach in examples and is illustrated in a biomechanical engineering application. Supplemental materials, including the software and a user manual, are available online and can be requested from the first author.

Complex computer simulation codes (computer experiments) are used extensively in scientific and engineering studies (see Santner, Williams, and Notz 2003; Fang, Li, and Sudjianto 2005; and Forrester, Sobester, and Keane 2008). Many computer codes were developed for applications where the results of physical experiments are also available to describe the true input/output relationship of interest. In this article we are interested in settings where both types of data are available, but the number of computer code runs is limited because they have lengthy run times. In such a case, statistical models are needed to extend the set of inputs at which output is available by predicting the response from the computer and physical experiments based on a limited number of training runs.

Such statistical models must incorporate different types of input variables. This article focuses on computer experiments having three types of inputs: control variables, tuning parameters, and calibration parameters. Control variables are (design) inputs to both the simulation code and the corresponding physical experiment (Santner, Williams, and Notz 2003, p. 15). A primary goal of the physical experiment is to study the relationship between the response and control variables.

Tuning parameters are inputs that are present only in the simulation code. Tuning parameters are typically numerical quantities that control the solution of a numerical algorithm implemented by a computer experiment. Thus, tuning parameters have *no* meaning in the physical experiment. For example, Rawlinson et al. (2006) described a finite element analysis (FEA) computer code simulating the forces and movements of a knee prosthesis under a given loading regimen. The two tuning parameters for their model were the mesh density of the tibial tray and the amount of discretization in the curve used to describe the knee loading. Tuning is the process of determining the values of the tuning parameters so that a computer simulation best represents a physical experiment.

Calibration parameters are also present only in the simulation code. Unlike tuning parameters, calibration parameters have meaning in the physical experiment, but their values are either unknown or unmeasured during the running of the physical experiment. For example, Kennedy and O’Hagan (2001) described a computer code for simulating the deposition of ruthenium 106 in the Tomsk-7 chemical plant, which experienced an accident in 1993. In their code, one calibration parameter was the deposition velocity. Calibration is the process of determining plausible values of the calibration parameters to match the results of the physical experiment.

Some computer experiments have both tuning and calibration parameters. For example, in the FEA code of Rawlinson et al. (2006) two calibration parameters were position of the femur relative to the tibia in the initial gait cycle and friction between the bone and the prosthesis.

Although considerable research has been dedicated to setting either tuning or calibration parameters, there has been little study on the problem of setting both simultaneously. In this article, we introduce a methodology for the simultaneous determination of tuning and calibration parameters. Section 2 reviews current methods for tuning and calibration and discusses their limitations. Section 3 introduces a hierarchical Bayesian model describing the response from a physical experiment together with the output from the computer simulation code having both tuning and calibration parameters. Section 4 proposes a methodology for simultaneously setting these parameters. Section 5 compares the proposed methodology with an approach that treats all parameters as calibration parameters.

Previous proposals for calibration have been primarily Bayesian. Their goal was to simulate the posterior distributions of the calibration parameters and to use this posterior distribution to better predict unknown responses and assess the uncertainty in the prediction. Craig et al. (1996, 2001) developed *a Bayesian linear forecasting method*, which can be seen as a pragmatic approximation to a fully Bayesian analysis. Kennedy and O’Hagan (2001) described a Bayesian calibration framework based on a model having a bias function and a (modular) Bayesian analysis. Higdon et al. (2004) proposed a fully Bayesian implementation for the model in Kennedy and O’Hagan (2001). Higdon et al. (2008) developed a methodology for calibration when the outputs are multivariate. Gattiker (2005) implemented the models in Higdon et al. (2004, 2008) in MATLAB. Loeppky, Bingham, and Welch (2006) proved that the previous procedures would lead to correct estimation of the calibration parameters asymptotically, provided the true values of the calibration parameters could reduce the bias of the computer simulation to zero. We discuss a lemma in Loeppky, Bingham, and Welch (2006) in Section 5.

Classical methods used to set tuning parameters, such as Park (1991) and Cox, Park, and Singer (1996), select these inputs to make the computer code output fit the physical observations as closely as possible in an integrated prediction error sense. These proposals for selecting tuning parameters do not provide estimated uncertainties associated with the selected values. However, several recent proposals for “validating” computer experiments are similar to the methodology for Bayesian calibration. For example, Bayarri et al. (2007b) presented a framework with a Bayesian model for validation and Bayarri et al. (2007a) used wavelet decomposition based on the same Bayesian model to validate computer experiments having multivariate outputs.

Because the methods above focus on either tuning or calibration, they may suffer limitations in applications where both tuning and calibration parameters are present. Specifically, if one applies classical tuning methodology to a problem involving both tuning and calibration parameters, one cannot quantify the uncertainties of the predicted calibration parameters. If one applies calibration methodology to such a problem and treats all parameters as calibration parameters, one can obtain biased parameter predictions and/or undesirable estimated uncertainties of the tuning and calibration parameters as we will show below. These misleading predictions and estimated uncertainties can lead to large prediction errors.

To demonstrate the limitations alluded to above, we used the Bayesian calibration program (gpmsa) described in Gattiker (2005) to set tuning and calibration parameters for an application in Rawlinson et al. (2006). Originally, the goal of this project was to compare the damage in two knee implants [the Install-Burstein (IB) manufactured by Zimmer, Inc. and the Optetrak produced by Exactech, Inc.] under certain loading regimens.

Here, we consider only the IB knee implant and a physical experiment in which the measurements consist of kinematics and kinetics of knees tested in an Instron-Stanmore KC1 testing device, a “knee simulator.” The loading (magnitude, angle, and rate), joint stiffness, and knee design were modeled in a finite element analysis (FEA) computer code whose output included the kinematics, kinetics, and stresses experienced by the knee component. The output we considered was the anterior–posterior displacement (APD) of the femoral component with respect to the tibial tray, which is an important kinematic output of fore-aft motion during knee function. The APD was measured in both the knee simulator and the computer simulation; it is roughly proportional to the anterior–posterior force that, coupled with the vertical joint load, contributed to damage of the prosthesis. Thus this computer experiment had one control variable (the percentile in gait cycle), two tuning parameters (finite element mesh density and load discretization), and two calibration parameters (friction and initial position).

The design of the computer and the physical experiments were roughly maximin Latin hypercube designs (McKay, Beckman, and Conover 1979). A pilot study of the APD for this design acquired 439 observations from the computer code and 36 observations from the physical experiment for the purpose of tuning and calibration. Figure 1 depicts the training data. A previous sensitivity analysis (Saltelli, Chan, and Scott 2000) found that APD was neither sensitive to friction nor mesh density. Therefore, we study load discretization (a tuning parameter) and initial position (a calibration parameter) regarding the percentage of the gait cycle as a design input to both physical and computer experiments.

Scatterplot of the measured APD against percentile in the gait cycle, from the knee simulator (triangles) and the FEA computer code (dots).

Using gpmsa (Gattiker 2005), we treated both load discretization and initial position as calibration parameters whose prior distributions were near uniform distribution. We ran gpmsa with 8000 burn-in iterations and 2000 production runs, and based our inferences on 100 equally spaced values of each parameter taken from the production runs. Figure 2 shows the simulated posterior distributions of these two parameters. The posterior distribution of load discretization is bimodal and that of initial position has a lower value in its range as mode, but has large variation. Thus, plausible values of the tuning and calibration parameters remain unclear. We demonstrate in Section 5.2 that our methodology can significantly improve the selection of the tuning and calibration parameters in this application.

We let **x** denote (the vector of) control inputs, **t** denote the vector of tuning parameters, and **t**^{} denote the “best” values of the tuning parameters, where **t**^{} will be defined in Section 4. We let **c** denote the vector of calibration parameters and **c**^{} denote the *unknown* true values of the calibration parameters. We assume **x** [0, 1]^{px}, **t** [0, 1]^{pt}, and **c** [0, 1]^{pc} or they can be so scaled. Information about the value of the unknown **c**^{} will be contained in its (marginal) prior, which will be denoted by [**c**^{}].

We let $\{({\mathbf{x}}_{i}^{s},{\mathbf{c}}_{i},{\mathbf{t}}_{i}),{y}^{s}({\mathbf{x}}_{i}^{s},{\mathbf{c}}_{i},{\mathbf{t}}_{i});i=1,2,\mathrm{\dots},{n}_{s}\}$ denote the training data from the computer experiment and $\{{\mathbf{x}}_{j}^{p},{y}^{p}({\mathbf{x}}_{j}^{p});j=1,2,\mathrm{\dots},{n}_{p}\}$ denote the training data from the physical experiment. Here *n _{s}* and

The model proposed here is similar in spirit to the ones in Kennedy and O’Hagan (2001) and Higdon et al. (2005) except that it explicitly incorporates both the tuning and calibration parameters. In other words, the response from the physical experiment is the sum of the “true response” and a random noise; the output from the computer experiment is the difference between the true response and the code bias.

Formally, we assume that the code output at (**x**, **c**, **t**), *y ^{s}*(

$${Y}^{S}(\mathbf{x},\mathbf{c},\mathbf{t})={\mathbf{f}}_{\mathbf{Z}}^{\top}(\mathbf{x},\mathbf{c},\mathbf{t}){\mathit{\beta}}_{Z}+Z(\mathbf{x},\mathbf{c},\mathbf{t}),$$

(1)

where ${\mathbf{f}}_{Z}^{\top}(\mathbf{x},\mathbf{c},\mathbf{t})$ is a vector of *known* regression functions and **β**_{Z} is a vector of *unknown* regression parameters. If *Y ^{s}*(•) is assumed to have a constant mean, then ${\mathbf{f}}_{Z}^{\top}(\mathbf{x},\mathbf{c},\mathbf{t}){\mathit{\beta}}_{Z}={\beta}_{Z,0}$. In Equation (1),

$$\begin{array}{l}{\text{Cor}(Z(\mathbf{x}}_{1},{\mathbf{c}}_{1},{\mathbf{t}}_{1}),Z{(\mathbf{x}}_{2},{\mathbf{c}}_{2},{\mathbf{t}}_{2})|{\rho}_{Z})\hfill \\ \begin{array}{cc}& ={\displaystyle \prod _{i=1}^{{p}_{x}}{\rho}_{Z,x,i}^{4{({x}_{1,i}-{x}_{2,i})}^{2}}}\times {\displaystyle \prod _{j=1}^{{p}_{c}}{\rho}_{Z,c,j}^{4{({c}_{1,j}-{c}_{2,j})}^{2}}}\times {\displaystyle \prod _{k=1}^{{p}_{t}}{\rho}_{Z,t,k}^{4{({t}_{1,k}-{t}_{2,k})}^{2}}},\end{array}\hfill \end{array}$$

(2)

with correlation parameters ρ_{Z} = (ρ*Z,x*,1, …, ρ*Z,x,p _{x}*, ρ

We regard the output of the physical experiment at **x**, *y ^{p}*(

$${Y}^{p}(\mathbf{x})=\eta (\mathbf{x},{\mathbf{c}}^{\u2605})+\epsilon ,$$

(3)

where η(**x, c**^{}) = *E*(*Y ^{p}*(

We define

$$\delta (\mathbf{x},\mathbf{c},\mathbf{t})=\eta (\mathbf{x},{\mathbf{c}}^{\u2605})-{y}^{s}(\mathbf{x},\mathbf{c},\mathbf{t})$$

(4)

to be the code bias at (**x, c, t**). Our model assumes δ(•) is a draw from the Gaussian stochastic process

$$\mathrm{\Delta}(\mathbf{x},\mathbf{c},\mathbf{t})={\mathbf{f}}_{D}^{\top}(\mathbf{x},\mathbf{c},\mathbf{t}){\mathit{\beta}}_{D}+D(\mathbf{x},\mathbf{c},\mathbf{t}),$$

(5)

where **f**_{D}(**x, c, t**) is a vector of known regression functions and β_{D} is a vector of unknown regression parameters; if Δ(•) is assumed to have a constant mean, then ${\mathbf{f}}_{D}^{\top}(\mathbf{x},\mathbf{c},\mathbf{t}){\mathit{\beta}}_{D}={\beta}_{D,0}$. Similar to *Z*(•), *D*(•) is defined to be a stationary Gaussian stochastic process with mean 0, variance ${\sigma}_{D}^{2}$, and correlation function

$$\begin{array}{l}\text{Cor}(D{(\mathbf{x}}_{1},{\mathbf{c}}_{1},{\mathbf{t}}_{1}),D{(\mathbf{x}}_{2},{\mathbf{c}}_{2},{\mathbf{t}}_{2})|{\mathit{\rho}}_{D})\hfill \\ \begin{array}{cc}& ={\displaystyle \prod _{i=1}^{{p}_{x}}{\rho}_{D,x,i}^{4{({x}_{1,i}-{x}_{2,i})}^{2}}}\times {\displaystyle \prod _{j=1}^{{p}_{c}}{\rho}_{D,c,j}^{4{({c}_{1,j}-{c}_{2,j})}^{2}}}\times {\displaystyle \prod _{k=1}^{{p}_{t}}{\rho}_{D,t,k}^{4{({t}_{1,k}-{t}_{2,k})}^{2}}},\end{array}\hfill \end{array}$$

(6)

where ρ_{D} = (ρ*D,x*,1, …, ρ*D,x,p _{x}*, ρ

In conclusion, conditional on the parameters $({\mathit{\beta}}_{Z},{\sigma}_{Z}^{2},{\mathit{\rho}}_{Z})$, our model assumes that for *given* inputs (**x, c, t**), the simulator output *y ^{s}*(

$${Y}^{s}(\mathbf{x},\mathbf{c},\mathbf{t})={\mathbf{f}}_{Z}^{\top}(\mathbf{x},\mathbf{c},\mathbf{t}){\mathit{\beta}}_{Z}+Z(\mathbf{x},\mathbf{c},\mathbf{t}).$$

(7)

Similarly conditional on unknown parameters, the physical experiment output at control input **x**, *y ^{p}*(

$$\begin{array}{ll}{Y}^{p}(\mathbf{x})\hfill & ={Y}^{s}(\mathbf{x},{\mathbf{c}}^{\u2605},\mathbf{t})=\mathrm{\Delta}(\mathbf{x},{\mathbf{c}}^{\u2605},\mathbf{t})+\mathit{\epsilon}\hfill \\ \hfill & =({\mathbf{f}}_{Z}^{\top}(\mathbf{x},{\mathbf{c}}^{\u2605},\mathbf{t}){\mathit{\beta}}_{Z}+Z(\mathbf{x},{\mathbf{c}}^{\u2605},\mathbf{t}))+({\mathbf{f}}_{D}^{\top}(\mathbf{x},{\mathbf{c}}^{\u2605},\mathbf{t}){\mathit{\beta}}_{D}+D({\mathbf{x},\mathbf{c}}^{\u2605},\mathbf{t})+\mathit{\epsilon},\hfill \end{array}$$

(8)

for any given vector **t** of tuning parameters. The unknown parameters being conditioned in Equation (8) are those of *Y ^{s}*(

Next we introduce our prior model for the unknown model parameters and derive the posteriors given the tuning parameters **t**. The proposed priors are similar to the ones in Higdon et al. (2004, 2005), and Gattiker (2005). For a generic parameter (or a vector of parameters) ** v**, we let [

We assume that the priors for the calibration parameters, **c**^{}, the process means, $\mathit{\beta}={({\mathit{\beta}}_{Z}^{\top},{\mathit{\beta}}_{D}^{\top})}^{\top}$, the process variances, $\mathit{\sigma}={({\sigma}_{Z}^{2},{\sigma}_{D}^{2},{\sigma}_{\epsilon}^{2})}^{\top}$, and process correlations, $\mathit{\rho}={({\mathit{\rho}}_{Z}^{\top},{\mathit{\rho}}_{D}^{\top})}^{\top}$, are independent so that

$$[{\mathbf{c}}^{\u2605},\mathit{\beta},\mathit{\sigma},\mathit{\rho}]=[{\mathbf{c}}^{\u2605}]\phantom{\rule{thinmathspace}{0ex}}\times \phantom{\rule{thinmathspace}{0ex}}[\mathit{\beta}]\phantom{\rule{thinmathspace}{0ex}}\times \phantom{\rule{thinmathspace}{0ex}}[\mathit{\sigma}]\phantom{\rule{thinmathspace}{0ex}}\times \phantom{\rule{thinmathspace}{0ex}}[\mathit{\rho}].$$

(9)

If expert knowledge about the priors is available it should be used. For example, if the researcher knows (roughly) the magnitude of the noise term, it should be used to construct the prior for ${\sigma}_{\epsilon}^{2}$. Lacking such knowledge, our methodology is based on the following priors.

The examples and software take ${\mathbf{f}}_{D}^{\top}(\u2022){\mathit{\beta}}_{Z}={\beta}_{Z,0}$ and ${\mathbf{f}}_{D}^{\top}(\u2022){\mathit{\beta}}_{D}={\beta}_{D,0}$. Our methodology for selecting (β_{Z,0}, β_{D,0}) is modular in the spirit of Bayarri et al. (2007a) in that these mean parameters are determined by non-Bayesian arguments in Section 4.2. All other parameters are determined from draws of the posterior.

Because the calibration parameters are scaled to [0, 1]^{pc}, all the elements of **c**^{} are taken to have independent and identical truncated normal distributions with mean 0.5, standard deviation 2, and support [0, 1]. This prior is close to a uniform distribution on [0, 1]^{pc}.

The variance parameters are taken to have independent inverse Gamma prior distributions. Letting IG(α, γ) denote the inverse Gamma distribution with mean 1/[γ(α – 1)] and variance 1/[γ^{2}(α – 1)^{2}(α – 2)], we take ${\sigma}_{Z}^{2}\sim \text{IG}({\alpha}_{Z},{\gamma}_{Z}),{\sigma}_{\epsilon}^{2}\sim \text{IG}({\alpha}_{\epsilon},{\gamma}_{\epsilon})$, and ${\sigma}_{D}^{2}\sim \text{IG}({\alpha}_{D},{\gamma}_{D})$ where the parameters are set data-adaptively, as follows. We first compute the sample variances of **y**^{s} and **y**^{p}, which we denote by ${\stackrel{\uff3e}{\sigma}}_{s}^{2}$ and ${\stackrel{\uff3e}{\sigma}}_{p}^{2}$ respectively. We take ${\sigma}_{Z}^{2}$ to have a (mildly) informative prior distribution by setting $({\alpha}_{Z},{\gamma}_{Z})=(10,0.1/{\stackrel{\uff3e}{\sigma}}_{s}^{2})$. We select the ${\sigma}_{\epsilon}^{2}$ prior to have a small magnitude and variation by taking $({\alpha}_{\epsilon},{\gamma}_{\epsilon})=(1,1000/{\stackrel{\uff3e}{\sigma}}_{p}^{2})$ as the default. The supplemental software package allows the user to specify a multiplier for ${\sigma}_{\epsilon}^{2}$ so that the program can suit physical experiments with random noise of different magnitudes. We use the value of $({\stackrel{\uff3e}{\sigma}}_{p}^{2}-{\stackrel{\uff3e}{\sigma}}_{s}^{2})$ to determine the prior for ${\sigma}_{D}^{2}$; if ${\stackrel{\uff3e}{\sigma}}_{p}^{2}>{\stackrel{\uff3e}{\sigma}}_{s}^{2},{\sigma}_{D}^{2}$ is set to have a mildly informative prior by taking $({\alpha}_{D},{\gamma}_{D})=(10,0.1/({\stackrel{\uff3e}{\sigma}}_{p}^{2}-{\stackrel{\uff3e}{\sigma}}_{s}^{2}))$ while ${\stackrel{\uff3e}{\sigma}}_{p}^{2}\le {\stackrel{\uff3e}{\sigma}}_{s}^{2},{\sigma}_{D}^{2}$ is set to have a prior with small magnitude and variation by selecting $({\alpha}_{D},{\gamma}_{D})=(1,100/{\stackrel{\uff3e}{\sigma}}_{s}^{2})$.

Finally, we take the priors of the correlation parameters to be independently and identically distributed Beta(1, 0.5) distributions. This Beta distribution has mean 2/3, mode close to 1, and a rather diffuse support over (0, 1); it encourages dependence of the output on each input, but the broad support allows the data a strong voice in determining the posteriors (Linkletter et al. 2006 and Higdon et al. 2008).

Let **ϕ** = (**β, σ, ρ**) [or **ϕ** = (**σ, ρ**) when β is fixed as below] and [**c**^{}, **ϕ**] denote the joint prior of **c**^{} and **ϕ**. For a fixed **t** and given the data (**y**^{s}, **y**^{p}), the joint posterior distribution can be written as

$$[{\mathbf{c}}^{\u2605},\mathit{\varphi}|{\mathbf{y}}^{s},{\mathbf{y}}^{p},\mathbf{t}]\phantom{\rule{thinmathspace}{0ex}}\propto [{\mathbf{y}}^{p},{\mathbf{y}}^{s}|{\mathbf{c}}^{\u2605},\mathit{\varphi},\mathbf{t}]\phantom{\rule{thinmathspace}{0ex}}\times \phantom{\rule{thinmathspace}{0ex}}[{\mathbf{c}}^{\u2605},\mathit{\varphi},\mathbf{t}]$$

(10)

$$\propto [{\mathbf{y}}^{p},{\mathbf{y}}^{s}|{\mathbf{c}}^{\u2605},\mathit{\varphi},\mathbf{t}]\phantom{\rule{thinmathspace}{0ex}}\times \phantom{\rule{thinmathspace}{0ex}}[{\mathbf{c}}^{\u2605},\mathit{\varphi}].$$

(11)

(See Appendix A for details of the posterior distribution.) We simulate draws from [**c**^{}, ϕ|**y**^{s}, **y**^{p}, **t**] by implementing a Metropolis–Hastings (MH) algorithm. Appendix A sketches this algorithm (see Higdon et al. 2004, 2005 for details about related MH algorithms). Using this model and the relationship in Equations (10) and (11), draws can be made from the joint posterior distribution [**c**^{}, ϕ|**y**^{s}, **y**^{p}, **t**] which forms the basis for our proposed methodology for simultaneous tuning and calibration in Section 4.

Our definition of the ideal tuning parameters is based on a function that measures the discrepancy between the physical experiment response and the computer code output. The discrepancy is a function of (**c, t**), but can be viewed as a function of **t** for given **c** = **c**^{}. Three of the many possible discrepancy functions are *L*_{2} discrepancy, ∫(η(**x, c**^{}) – *y ^{s}*(

If **c**^{} is known, the *L*_{2} discrepancy is defined by

$$\int {(\eta ({\mathbf{x},\phantom{\rule{thinmathspace}{0ex}}\mathbf{c}}^{\u2605})-{y}^{s}({\mathbf{x},\phantom{\rule{thinmathspace}{0ex}}\mathbf{c}}^{\u2605},\mathbf{t}))}^{2}}d\mathbf{x}={\displaystyle \int {\delta}^{2}}({\mathbf{x},\phantom{\rule{thinmathspace}{0ex}}\mathbf{c}}^{\u2605},\mathbf{t})\phantom{\rule{thinmathspace}{0ex}}d\mathbf{x},$$

(12)

regarded as a function of **t**. If **c**^{} is unknown, but has a posterior distribution [**c**^{}|**y**^{s}, **y**^{p}, **t**], the *L*_{2} discrepancy becomes

$$\int \int {\delta}^{2}}({\mathbf{x},\phantom{\rule{thinmathspace}{0ex}}\mathbf{c}}^{\u2605},\mathbf{t})[{\mathbf{c}}^{\u2605}|{\mathbf{y}}^{s},{\mathbf{y}}^{p},\mathbf{t}]d{\mathbf{c}}^{\u2605}d\mathbf{x}.$$

Our objective in tuning is to estimate

$${\mathbf{t}}^{\u2605}=\underset{\mathbf{t}}{\text{argmin}}{\displaystyle \int \int {\delta}^{2}}({\mathbf{x},\mathbf{c}}^{\u2605},\mathbf{t})[{\mathbf{c}}^{\u2605}|{\mathbf{y}}^{s},{\mathbf{y}}^{p},\mathbf{t}]d{\mathbf{c}}^{\u2605}d\mathbf{x}.$$

(13)

Our objective in calibration is to compute the posterior distribution [**c**^{}|**y**^{s}, **y**^{p}, **t**^{}]. We achieve the two objectives by minimizing an estimated squared discrepancy (defined in Section 4.2) and simulating the posterior distribution of **c**^{} given the minimizer.

To review, the main idea is to choose the most plausible values of the tuning parameters so that the computer code can best represent the physical experiment and use these tuning parameters in the inferences about calibration parameters. Next, we describe our method for estimating δ^{2}(**x, c**^{},**t**), **t**^{}, and [**c**^{}|**y**^{s}, **y**^{p}, **t**^{}] which we call the STaC procedure, which stands for simultaneous tuning and calibration. Then we summarize the methodology in a three step procedure.

Our first of the three steps is to estimate δ^{2}(**x**, **c**^{},**t**) by the posterior mean *E*(Δ^{2}(**x**, **c**^{}, **t**)|**y**^{s}, **y**^{p}, **t**), which is the Bayes predictor of δ^{2}(**x, c**^{},**t**) under *squared error loss* (Santner, Williams, and Notz 2003). Notice that *E*(Δ^{2}(**x**, **c**^{}, **t**)| **c**^{}, ϕ, **y**^{s}, **y**^{p}, **t**) is equal to

$${[E(\mathrm{\Delta}(\mathbf{x},{\mathbf{c}}^{\u2605},\mathbf{t})|{\mathbf{c}}^{\u2605},\varphi ,{\mathbf{y}}^{s},{\mathbf{y}}^{p},\mathbf{t})]}^{2}+\text{Var}(\mathrm{\Delta}(\mathbf{x},{\mathbf{c}}^{\u2605},\mathbf{t})|{\mathbf{c}}^{\u2605},\varphi ,{\mathbf{y}}^{s},{\mathbf{y}}^{p},\mathbf{t}).$$

(14)

Closed form expresssions for both *E*(Δ(**x, c**^{}, **t**)|**c**^{}, ϕ, **y**^{s}, **y**^{p}, **t**) and Var(Δ(**x, c**^{}, **t**)|**c**^{}, ϕ, **y**^{s}, **y**^{p}, **t**) can be derived analytically based on the probability density of the multivariate normal distribution (see Appendix B).

The two terms in Equation (14) can be interpreted intuitively, as follows. Large values of [*E*(Δ(**x, c**^{}, **t**)|**c**^{}, ϕ, **y**^{s}, **y**^{p}, **t**)]^{2} imply that there is a large discrepancy between the computer and the physical experiments; thus [*E*(Δ(**x, c**^{}, **t**)|**c**^{}, ϕ, **y**^{s}, **y**^{p}, **t**)]^{2} can be interpreted as *squared prediction bias*. The term Var(Δ(**x, c**^{}, **t**)|**c**^{}, ϕ, **y**^{s}, **y**^{p}, **t**) is a measure of the uncertainty of the predicted bias. Thus Equation (14) favors a **t** that minimizes the sum of the squared prediction bias and its predicted uncertainty.

The right-hand side of

$$\begin{array}{c}E({\mathrm{\Delta}}^{2}(\mathbf{x},{\mathbf{c}}^{\u2605},\mathbf{t})|{\mathbf{y}}^{s},{\mathbf{y}}^{p},\mathbf{t})\hfill \\ \hfill ={E}_{[{\mathbf{c}}^{\u2605},\mathit{\varphi}|{\mathbf{y}}^{s},{\mathbf{y}}^{p},\mathbf{t}]}[E({\mathrm{\Delta}}^{2}(\mathbf{x},{\mathbf{c}}^{\u2605},\mathbf{t})|{\mathbf{c}}^{\u2605},\mathit{\varphi},{\mathbf{y}}^{s},{\mathbf{y}}^{p,}\mathbf{t})],\end{array}$$

can be approximated using Equation (14), and based on draws
$(\stackrel{\uff3e}{{\mathbf{c}}_{\ell}^{\u2605}},\stackrel{\uff3e}{{\mathit{\varphi}}_{\ell}})$ from [**c**^{}, ϕ|**y**^{s}, **y**^{p}, **t**], which can be obtained using the MH algorithm. The law of large numbers guarantees that

$$\begin{array}{c}E({\mathrm{\Delta}}^{2}(\mathbf{x},{\mathbf{c}}^{\u2605},\mathbf{t})|{\mathbf{y}}^{s},{\mathbf{y}}^{p},\mathbf{t})\hfill \\ \hfill \approx \frac{1}{{N}_{MC}}{\displaystyle \sum _{\ell =1}^{{N}_{MC}}E({\mathrm{\Delta}}^{2}(\mathbf{x},\underset{\ell}{\overset{\uff3e}{{\mathbf{c}}^{\u2605}}},\mathbf{t})|}\underset{\ell}{\overset{\uff3e}{{\mathbf{c}}^{\u2605}}},\stackrel{\uff3e}{{\mathit{\varphi}}_{\ell}},{\mathbf{y}}^{s},{\mathbf{y}}^{p},\mathbf{t}),\end{array}$$

(15)

for a sufficiently large *N _{MC}*. Thus, we can approximate

The estimation of δ^{2}(**x, c**^{},**t**) using Equation (15) requires knowledge of (β_{Z,0}, β_{D,0}); these parameters are set as follows. The parameter β_{D,0} can be interpreted as the difference between the average of the physical experiment and the average of the computer experiment outputs when the calibration and tuning parameters are equal to (**c**^{},**t**). The value of β_{D,0} can impact both calibration and tuning and thus must be carefully chosen. Assuming expert knowledge about the difference β_{D,0} is not available, we suggest setting β_{D,0} based on the assumption that, for any given **t**, the code inadequacy is decreased as **c** gets closer to **c**^{} so that β_{D,0} becomes smaller. Thus β_{D,0} is minimized when **c** = **c**^{}. Hence we take β_{D,0} to be the minimized difference, over **c** values, between the sample mean of the *n _{p}* physical experiment responses and the sample mean of

In the second step, we estimate **t**^{} by

$$\begin{array}{cc}\stackrel{\uff3e}{{\mathbf{t}}^{\u2605}}\hfill & =\underset{\mathbf{t}}{\text{argmin}}{\displaystyle \int \int E}[{\mathrm{\Delta}}^{2}({\mathbf{x},\mathbf{c}}^{\u2605},\mathbf{t})|{\mathbf{y}}^{s},{\mathbf{y}}^{p},\mathbf{t}]\phantom{\rule{thinmathspace}{0ex}}\times {[\mathbf{c}}^{\u2605},\mathit{\varphi}|{\mathbf{y}}^{s},{\mathbf{y}}^{p},\mathbf{t}]d{\mathbf{c}}^{\u2605}d\mathbf{x}\hfill \\ \hfill & \approx {\displaystyle \underset{\mathbf{t}}{\text{argmin}}\int \frac{1}{{N}_{MC}}}\times {\displaystyle \sum _{\ell =1}^{{N}_{MC}}E}[{\mathrm{\Delta}}^{2}(\mathbf{x},\stackrel{\uff3e}{{\mathbf{c}}_{\ell}^{\u2605}},\mathbf{t})|\stackrel{\uff3e}{{\mathbf{c}}_{\ell}^{\u2605}},{\stackrel{\uff3e}{\mathit{\varphi}}}_{\ell},{\mathbf{y}}^{s},{\mathbf{y}}^{p},\mathbf{t}]d\mathbf{x}.\hfill \end{array}$$

(16)

The integral over **x** in Equation (16) can be approximated by a further Monte Carlo integration that averages Equation (15) over a grid of **x** inputs. In our examples, we took

$$\stackrel{\uff3e}{{\mathbf{t}}^{\u2605}}\approx \underset{\mathbf{t}}{\text{argmin}}\frac{1}{{N}_{x}}\frac{1}{{N}_{MC}}\times {\displaystyle \sum _{i=1}^{{N}_{x}}{\displaystyle \sum _{\ell =1}^{{N}_{MC}}E[{\mathrm{\Delta}}^{2}({\mathbf{x}}_{i},\stackrel{\uff3e}{{\mathbf{c}}_{\ell}^{\u2605}},\mathbf{t})|\stackrel{\uff3e}{{\mathbf{c}}_{\ell}^{\u2605}},\stackrel{\uff3e}{{\mathit{\varphi}}_{\ell}},{\mathbf{y}}^{s},{\mathbf{y}}^{p},\mathbf{t}],}}$$

(17)

with *N _{x}*= 101 and an equally spaced grid of values of

In the third step, we estimate the posterior distribution [**c**^{}|**y**^{s}, **y**^{p}, **t**^{}]. Because
$\stackrel{\uff3e}{{t}^{\u2605}}$ converges to **t**^{} almost surely, [**c**^{}|**y**^{s}, **y**^{p}, $\stackrel{\uff3e}{{t}^{\u2605}}$] converges to [**c**^{}|**y**^{s}, **y**^{p}, **t**^{}] in distribution. We therefore estimate the distribution [**c**^{}|**y**^{s}, **y**^{p}, **t**^{}] using the draws of **c**^{} from [**c**^{}, ϕ|**y**^{s}, **y**^{p},
$\stackrel{\uff3e}{{t}^{\u2605}}$].

In summary, we conduct simultaneous tuning and calibration as follows:

- Step 1. For each possible
**t**in a grid of the tuning parameter vectors, set the values of (β_{Z,0}, β_{D,0}), make draws $\{(\stackrel{\uff3e}{{\mathbf{c}}_{l}^{\u2605}},\stackrel{\uff3e}{{\mathit{\varphi}}_{l}});l=1,\mathrm{\dots},{N}_{MC}\}$ from [**c**^{},ϕ |**y**^{s},**y**^{p},**t**], and compute the estimated squared discrepancy at a set of inputs {**x**_{i}:*i*=1, …,*N*}._{x} - Step 2. Compute $\stackrel{\uff3e}{{\mathbf{t}}^{\u2605}}=\text{argmin}{}_{\mathbf{t}}\frac{1}{{N}_{x}}\frac{1}{{N}_{MC}}{\displaystyle {\sum}_{i=1}^{{N}_{x}}{\displaystyle {\sum}_{\ell =1}^{{N}_{MC}}E({\mathrm{\Delta}}^{2}({\mathbf{x}}_{i},}}\stackrel{\uff3e}{{\mathbf{c}}_{\ell}^{\u2605}},\mathbf{t})|\stackrel{\uff3e}{{\mathbf{c}}_{\ell}^{\u2605}},\stackrel{\uff3e}{{\mathit{\varphi}}_{\ell}},{\mathbf{y}}^{s},{\mathbf{y}}^{p},\mathbf{t})$
- Step 3. Estimate $[{\mathbf{c}}^{\u2605}|{\mathbf{y}}^{s},{\mathbf{y}}^{p},\stackrel{\uff3e}{{\mathbf{t}}^{\u2605}}]$.

In addition to tuning and calibration, we can use $\stackrel{\uff3e}{{t}^{\u2605}}$ and the draws ${\{(\stackrel{\uff3e}{{\mathbf{c}}_{\ell}^{\u2605}},\stackrel{\uff3e}{{\mathit{\varphi}}_{\ell}})\}}_{\ell}$ from $[({\mathbf{c}}^{\u2605},\mathit{\varphi})|{\mathbf{y}}^{s},{\mathbf{y}}^{p},\stackrel{\uff3e}{{\mathbf{t}}^{\u2605}}]$ to predict *y ^{s}*(•) and η(•) at arbitrary inputs and to construct predictive intervals for these quantities. Given (

$$\stackrel{\uff3e}{{y}^{s}}(\mathbf{x},\mathbf{c},\mathbf{t})=\frac{1}{{N}_{MC}}{\displaystyle \sum _{\ell =1}^{{N}_{MC}}E[{Y}^{s}}(\mathbf{x},\mathbf{c},\mathbf{t})|{\mathbf{y}}^{s},{\mathbf{y}}^{p},\stackrel{\uff3e}{{\mathbf{t}}^{\u2605}},\stackrel{\uff3e}{{\mathbf{c}}_{\ell}^{\u2605}},\u2605{\stackrel{\uff3e}{\mathit{\varphi}}}_{\ell})$$

(18)

is an approximation to the BLUP of *y ^{s}*(

$$\stackrel{\uff3e}{\eta}({\mathbf{x},\mathbf{c}}^{\u2605})=\frac{1}{{N}_{MC}}{\displaystyle \sum _{\ell =1}^{{N}_{MC}}E({Y}^{s}}(\mathbf{x},\stackrel{\uff3e}{{\mathbf{c}}_{\ell}^{\u2605}},\stackrel{\uff3e}{{\mathbf{t}}^{\u2605}})+\mathrm{\Delta}(\mathbf{x},\phantom{\rule{thinmathspace}{0ex}}\stackrel{\uff3e}{{\mathbf{c}}_{\ell}^{\u2605}},\stackrel{\uff3e}{{\mathbf{t}}^{\u2605}})|{\mathbf{y}}^{s},{\mathbf{y}}^{p},\stackrel{\uff3e}{{\mathbf{t}}^{\u2605}},\stackrel{\uff3e}{{\mathbf{c}}_{\ell}^{\u2605}},\stackrel{\uff3e}{{\mathit{\varphi}}_{\ell}})$$

(19)

is an approximation to the BLUP of η(**x, c**^{}) by Equations (3) through (5).

We quantify the uncertainty in the $\stackrel{\uff3e}{{y}^{s}}(\mathbf{x},\mathbf{c},\mathbf{t})$ predictor using

$$\begin{array}{c}\text{Var(}{Y}^{s}(\mathbf{x},\mathbf{c},\mathbf{t})|{\mathbf{y}}^{s},{\mathbf{y}}^{p},\stackrel{\uff3e}{{\mathbf{t}}^{\u2605}})\hfill \\ \hspace{1em}={E}_{[{\mathbf{c}}^{\u2605},\mathit{\varphi}|{\mathbf{y}}^{s},{\mathbf{y}}^{p},\stackrel{\uff3e}{{\mathbf{t}}^{\u2605}}]}[\text{Var}({Y}^{s}(\mathbf{x},\mathbf{c},\mathbf{t})|{\mathbf{y}}^{s},{\mathbf{y}}^{p},\stackrel{\uff3e}{{\mathbf{t}}^{\u2605}},{\mathbf{c}}^{\u2605},\mathit{\varphi})]+{\text{Var}}_{[{\mathbf{c}}^{\u2605},\mathit{\varphi}|{\mathbf{y}}^{s},{\mathbf{y}}^{p},\stackrel{\uff3e}{{\mathbf{t}}^{\u2605}}]}[E({Y}^{s}(\mathbf{x},\mathbf{c},\mathbf{t})|{\mathbf{y}}^{s},{\mathbf{y}}^{p},\stackrel{\uff3e}{{\mathbf{t}}^{\u2605}},{\mathbf{c}}^{\u2605},\mathit{\varphi})],\hfill \end{array}$$

(20)

which can be approximated using draws from the MH algorithm. We provide details in Appendix C. Similarly, we quantify the uncertainty in the (**x, c**^{}) predictor using $\text{Var}({Y}^{s}(\mathbf{x},{\mathbf{c}}^{\u2605},\stackrel{\uff3e}{{\mathbf{t}}^{\u2605}})+\mathrm{\Delta}(\mathbf{x},{\mathbf{c}}^{\u2605},\stackrel{\uff3e}{{\mathbf{t}}^{\u2605}})|{\mathbf{y}}^{s},{\mathbf{y}}^{p},\stackrel{\uff3e}{{\mathbf{t}}^{\u2605}})$. Details of this computation are also given in Appendix C.

Based on Equations (18), (19), (20), and the results in Appendix C, we can compute point-wise 100(1 – α)% prediction bands for *y ^{s}*(

$$\stackrel{\uff3e}{{y}^{s}}(\mathbf{x},\phantom{\rule{thinmathspace}{0ex}}\mathbf{c},\phantom{\rule{thinmathspace}{0ex}}\mathbf{t})\pm {z}^{\alpha /2}\sqrt{\text{Var}({Y}^{s}(\mathbf{x},\phantom{\rule{thinmathspace}{0ex}}\mathbf{c},\phantom{\rule{thinmathspace}{0ex}}\mathbf{t})|{\mathbf{y}}^{s},{\mathbf{y}}^{p},\stackrel{\uff3e}{{\mathbf{t}}^{\u2605}})}$$

and

$$\stackrel{\uff3e}{\eta}(\mathbf{x},{\mathbf{c}}^{\u2605})\pm {z}^{\alpha /2}\sqrt{\text{Var}({Y}^{s}(\mathbf{x},{\mathbf{c}}^{\u2605},\phantom{\rule{thinmathspace}{0ex}}\stackrel{\uff3e}{{\mathbf{t}}^{\u2605}}+\mathrm{\Delta}({\mathbf{x},\phantom{\rule{thinmathspace}{0ex}}\mathbf{c}}^{\u2605}\stackrel{\uff3e}{{\mathbf{t}}^{\u2605}})|{\mathbf{y}}^{s},{\mathbf{y}}^{p},\stackrel{\uff3e}{{\mathbf{t}}^{\u2605}})},$$

where *z*^{α/2} is the upper α/2 critical point of the standard normal distribution.

In this section, we compare Bayesian calibration as implemented by the gpmsa software with STaC. Recall that the former performs a fully Bayesian calibration for the tuning as well as the calibration parameters, whereas STaC sets the tuning parameters to minimize a predicted *L*_{2} discrepancy between the computer and the physical experiment outputs. Lemma 1 in Loeppky, Bingham, and Welch (2006) shows that, roughly, for their Gaussian stochastic process model, under the conditions (a) the discrepancy between the computer and the physical experiments can be minimized to be zero and (b) infinitely many runs from the computer simulation and the physical experiment are available, the maximum likelihood estimator (MLE) of the discrepancy is zero and the MLE of the calibration parameters minimizes the discrepancy asymptotically. Based on this lemma, Bayesian calibration and STaC can give comparable estimates of **t**^{} in our model if (A) the discrepancy δ(•) can be made sufficiently close to zero, (B) the number of observations from the computer experiment is reasonably large, and (C) the posterior mode of every calibration parameter is close to the MLE of the parameter.

However, (A), (B), and (C) can be hard or impossible to verify/implement in many applications for at least three reasons. It may not be possible to drive the the discrepancy δ(•) to zero because of the simplified physics or biology used to construct the simulation code. The number of computer and physical experiment runs are typically limited. The posterior mode of a calibration parameter can differ from the MLE if either the prior for **c**^{} or Δ(•) is informative. Therefore, the Bayesian calibration approach may fail to minimize the difference between the computer and the physical experiments. Next we compare STaC with gpmsa in a numerical example where the bias can be minimized to zero, yet the training data are scarce.

This example has three real-valued inputs to a computer code: one control variable *x*, one tuning parameter *t*, and one calibration parameter *c*. All three inputs have the support [0, 1]. The output from the computer code was generated as

$${y}^{s}(x,c,t={x}^{2}+c\times (x+1)+0.5\times t,$$

and the response from the physical experiment as

$${y}^{p}(x)={x}^{2}+0.1\times x+0.5\times \epsilon ,$$

where ε denotes noise generated from normal distribution having mean zero and standard deviation 0.01. In this example the bias term can be driven to zero by setting *c* = 0.1 and *t* = 0.8. Thus any approach that suggests *t*^{} = 0.8 and *c*^{} = 0.1 is favorable.

We generated 30 inputs for the computer experiment and 5 for the physical experiment with both designs being maximin Latin hypercube designs. The Bayesian calibration program gpmsa was run treating *c* and *t* as calibration parameters. Following 8000 burn-in iterations and 2000 production ones, we regarded 100 equally spaced samples drawn from the 2000 production runs as draws from the joint posterior distribution of (*t*^{}, *c*^{}). The estimated posterior distributions of *t*^{} and *c*^{} are shown in Figure 3. The estimated [*t*^{}|**y**^{p}, **y**^{s}] is not able to pinpoint that *t*^{} = 0.8 and the estimated [*c*^{}|**y**^{p}, **y**^{s}] fails to suggest values close to 0.1. We next used gpmsa to predict η(*x*, *c*^{}) over the grid, *x* = 0, 0.02, …, 1, of equally spaced inputs. Figure 4 depicts the training data, the true response curve, and the predictions. We see that the predictions are biased when *x* is near 0. The root mean squared prediction error (RMSPE) was computed as a measure of the predictive accuracy; the RMSPE of the generic predictor (•) with respect to the *x* grid above is defined to be $\sqrt{\frac{1}{51}{\displaystyle {\sum}_{i=1}^{51}{(\stackrel{\uff3e}{\eta}({x}_{i},{c}^{\u2605})-\eta ({x}_{i},{c}^{\u2605}))}^{2}}}$. The RMSPE of the predictor obtained by gpmsa is 0.0693.

Simulated posterior distributions of *t*^{} (the left panel) and *c*^{} (the right panel) using the Bayesian calibration program.

The training data (solid circles), the true response curve (the solid line), and the predictions (pluses) obtained by gpmsa.

We applied STaC to this example with *N _{MC}* = 100 and

Results of tuning and calibration from the STaC program: The left panel plots of the estimated *L*_{2} discrepancy against *t* where *t* = 0, 0.1, …, 0.9, 1; the right panel is a histogram of the draws from the posterior distribution of *c*^{}.

Figure 6 plots the predictions and the 99% two-sided prediction bands for true input-output function η(*x*, *c*^{}) using STaC. We see that the predictions are close to the true responses and the 99% prediction band contains the true curve for all *x* [0, 1]. The RMSPE of the predictor obtained by the STaC program is 0.0382. The relative improvement of STaC compared with the gpmsa-based predictor is 44.8% [= (0.0693 – 0.0382)/0.0693 × 100%].

The training data (solid circles), true response curve (the solid line), predictions (pluses), and 99% prediction bands (dashes) using the STaC program.

Thus, for this example, the inferences about the tuning and calibration parameters and the predictions with STaC are more informative and more accurate than using a fully Bayesian procedure.

We apply STaC to the biomechanics application described in Section 2. The number of burn-in and production runs are identical to those used in our first example. The left panel of Figure 7 plots the estimated squared discrepancies for the tuning parameter equal to 15, 16, …, 30 (with *N _{MC}* = 100 and

Scatterplot of the estimated squared discrepancy versus the value of the tuning parameter (left panel); histogram of the simulated posterior distribution of the calibration parameter when $\stackrel{\uff3e}{{t}^{\u2605}}=19$ (right panel).

Figure 8 depicts the predictions and 99% prediction bands for the anterior–posterior displacement in the physical experiment using the current data. Once the initial position and load discretization are set using STaC, the predictions clearly match those provided by the knee simulator.

This article introduces a modular Bayesian methodology for simultaneous tuning and calibration (STaC). It demonstrates, with examples, that STaC performs well for tuning, calibration, and prediction of the true input–output relationship. The examples also show that treating tuning and calibration parameters differentially is beneficial whether or not the code bias can be reduced to zero.

There are several areas of continuing research interest concerning STaC. Currently, determining the tuning parameters in the first step of STaC can be time consuming; additional work on speeding up the computation is needed. Second, a measure of the uncertainty in the estimated tuning parameter $\stackrel{\uff3e}{{t}^{\u2605}}$, perhaps based on the curvature of the integral ${\int}_{0}^{1}{\displaystyle {\int}_{0}^{1}{\delta}^{2}}(\mathbf{x},{\mathbf{c}}^{\u2605},\mathbf{t}})[{\mathbf{c}}^{\u2605}|{\mathbf{y}}^{s},{\mathbf{y}}^{p},\mathbf{t}]d{\mathbf{c}}^{\u2605}d\mathbf{x$, is needed. Third, it is of interest to extend STaC to other settings, for example, to cases where outputs are multivariate and to applications where both the computer code and the physical experiment have mixed quantitative and qualitative inputs. Fourth, there are other estimates of δ that can produce superior results for tuning and thus deserve further investigation. For example, another estimate of δ(**x, c**^{}, **t**) is *E*(*Y ^{p}*(

Lastly, we sketch two research areas related to calibration and tuning. One is to determine better experimental designs for the inputs to the computer and physical experiments when STaC is the objective, i.e., inputs $\{({\mathbf{x}}_{i}^{s},{\mathbf{c}}_{i},{\mathbf{t}}_{i});i=1,2,\mathrm{...},{n}_{s}\}$. We use two separate Latin hypercube designs for the physical and computer experiments in our examples. We believe nested Latin hypercube designs, for example, orthogonal array based nested Latin hypercube designs proposed by Qian, Ai, and Wu (2009), Qian (2009), and Qian and Wu (2009), may better determine tuning and calibration parameters. The second research area is to solve the potential issue of nonidentifiability. Note that tuning parameters, calibration parameters, and the code bias δ may be correlated. For example, if *y ^{p}*(

This research was sponsored, in part, by the National Science Foundation under Agreements DMS-0112069 (SAMSI/Duke University) and DMS-0406026 (The Ohio State University). Research support for the development of the biomechanics code and knee simulator used to produce the data for the example in Sections 2 and 5.3 was provided by NIH grant AR049793, Exactech, Inc., and the Clark and Kirby Foundations. The authors thank the editor, associate editor, and two reviewers for their constructive suggestions and valuable comments.

The posterior density is proportional to the product of the joint prior density [**c**^{}, **ϕ**|**t**] and the likelihood [**y**^{s}, **y**^{p}|(**c**^{}, **ϕ**, **t**)]. We can determine [**c**^{}, **ϕ**|**t**] by computing the prior densities of [**c**^{}], [**β**], [**σ**], and [**ρ**] and applying Equation (9). The density [**y**^{s}, **y**^{p}|(**c**^{}, **ϕ**, **t**)] is analytically derivable because **y**^{s} and **y**^{p} are viewed as observations from Gaussian stochastic processes. Specifically, our model is ${Y}^{s}({\mathbf{x}}_{i}^{s},{\mathbf{c}}_{i},{\mathbf{t}}_{i})={\mathbf{f}}_{Z}^{\top}({\mathbf{x}}_{i}^{s},{\mathbf{c}}_{i},{\mathbf{t}}_{i}){\mathit{\beta}}_{Z}+Z({\mathbf{x}}_{i}^{s},{\mathbf{c}}_{i},{\mathbf{t}}_{i})$ and ${Y}^{p}({\mathbf{x}}_{j}^{p})={Y}^{s}({\mathbf{x}}_{j}^{p},{\mathbf{c}}^{\u2605},\mathbf{t})+\mathrm{\Delta}({\mathbf{x}}_{j}^{p},{\mathbf{c}}^{\u2605},\mathbf{t})+\epsilon $ for all *i* = 1, …, *n _{s}* and

$$\begin{array}{c}\text{Cov}[{Y}^{s}({\mathbf{x}}_{i}^{s},{\mathbf{c}}_{i},{\mathbf{t}}_{i}),{Y}^{p}({\mathbf{x}}_{j}^{p})|({\mathbf{c}}^{\u2605},\mathit{\varphi},\mathbf{t})]\hfill \\ \hfill =\text{Cov}[Z({\mathbf{x}}_{i}^{s},{\mathbf{c}}_{i},{\mathbf{t}}_{i}),Z({\mathbf{x}}_{j}^{p},{\mathbf{c}}^{\u2605},\mathbf{t})|({\mathbf{c}}^{\u2605},\mathit{\varphi},\mathbf{t})],\end{array}$$

which can be computed using Equation (2).

With the joint prior density and the likelihood, the Metropolis–Hastings algorithm updates each of the parameters using the specifications listed in Table 1. For all the parameters, the proposal distributions are uniform distributions with centers set to be the previous draws and widths specified by the last column of Table 1. However, for all correlation parameters ρ, the proposal draws are made on ξ = −4 log(ρ).

Using Equations (5) and (6) we compute
${[E(\mathrm{\Delta}(\mathbf{x},{\mathbf{c}}^{\u2605},\mathbf{t})|{\mathbf{c}}^{\u2605},\mathit{\varphi},{\mathbf{y}}^{s},{\mathbf{y}}^{p},\mathbf{t})]}^{2}={[{\mathbf{f}}_{Z}^{\top}(\mathbf{x},{\mathbf{c}}^{\u2605},\mathbf{t}){\mathit{\beta}}_{D}+{\displaystyle {\sum}_{\mathit{\text{DZ}}}^{\top}}{\displaystyle {\sum}_{\mathit{\text{ZZ}}}^{-1}}{\mathbf{D}}_{\mathbf{t}}]}^{2}$ and $\text{Var}(\mathrm{\Delta}(\mathbf{x},{\mathbf{c}}^{\u2605},\mathbf{t})|{\mathbf{c}}^{\u2605},\mathit{\varphi},{\mathbf{y}}^{s},{\mathbf{y}}^{p},\mathbf{t})={\sigma}_{D}^{2}-{\displaystyle {\sum}_{\mathit{\text{DZ}}}^{\top}{\displaystyle {\sum}_{\mathit{\text{ZZ}}}^{-1}{\mathrm{\Sigma}}_{\mathit{\text{DZ}}}}}$ where **∑**_{DZ} is the (*n _{p}* +

- (B.1) For
*i*= 1,2, …,*n*, the_{p}*i*th element of**∑**_{DZ}is $\text{Cov}(\mathrm{\Delta}(\mathbf{x},{\mathbf{c}}^{\u2605},\mathbf{t}),{Y}^{p}({\mathbf{x}}_{i}^{p})|({\mathbf{c}}^{\u2605},\mathit{\varphi},\mathbf{t}))=\text{Cov}(\mathrm{\Delta}(\mathbf{x},\mathbf{t},{\mathbf{c}}^{\u2605}),\mathrm{\Delta}({\mathbf{x}}_{i}^{p},{\mathbf{c}}^{\u2605},\mathbf{t})|(\mathbf{t},{\mathbf{c}}^{\u2605},\mathit{\varphi}))$. For*i*=*n*+ 1,_{p}*n*+ 2, …,_{p}*n*+_{p}*n*, the_{s}*i*th element of**∑**_{DZ}is $\text{Cov}(\mathrm{\Delta}(\mathbf{x},{\mathbf{c}}^{\u2605},\mathbf{t}),{Y}^{s}({\mathbf{x}}_{i-{n}_{p}}^{s},{\mathbf{c}}_{i-{n}_{p}},{\mathbf{t}}_{i-{n}_{p}})|({\mathbf{c}}^{\u2605},\mathit{\varphi},\mathbf{t}))=0$. - (B.2) where$${\mathbf{\Sigma}}_{\mathit{\text{ZZ}}}=\left(\begin{array}{cc}{\mathbf{\Sigma}}_{\mathit{\text{PP}}}\hfill & {\mathbf{\Sigma}}_{\mathit{\text{PS}}}\hfill \\ {\mathbf{\Sigma}}_{\mathit{\text{PS}}}^{\top}\hfill & {\mathbf{\Sigma}}_{\mathit{\text{SS}}}\hfill \end{array}\right),$$(B.1)
- –
**Σ**_{PP}is an*n*×_{p}*n*matrix whose (_{p}*i*_{1},*i*_{2})th entry is $\text{Cov}({Y}^{p}({\mathbf{x}}_{{i}_{1}}^{p}),{Y}^{p}({\mathbf{x}}_{{i}_{2}}^{p})|({\mathbf{c}}^{\u2605},\mathit{\varphi},\mathbf{t}))$, for*i*,*i*_{1},*i*_{2}= 1, 2, …,*n*;_{p} - –
**Σ**_{PS}is an*n*×_{p}*n*matrix whose (_{s}*i,j*)th entry is $\text{Cov}({Y}^{p}({\mathbf{x}}_{i}^{p}),{Y}^{s}({\mathbf{x}}_{j}^{s},{\mathbf{c}}_{j},{\mathbf{t}}_{j})|({\mathbf{c}}^{\u2605},\mathit{\varphi},\mathbf{t}))$, for*i*= 1, 2, …,*n*and_{p}*j*= 1, 2, …,*n*;_{s} - –
**Σ**_{SS}is an*n*×_{s}*n*matrix whose (_{s}*j*_{1},*j*_{2})th entry is $\text{Cov}({Y}^{s}({\mathbf{x}}_{j1}^{s},{\mathbf{c}}_{j1},{\mathbf{t}}_{j1}),{Y}^{s}({\mathbf{x}}_{j2}^{s},{\mathbf{c}}_{j2},{\mathbf{t}}_{j2})|({\mathbf{c}}^{\u2605},\mathit{\varphi},\mathbf{t}))$, for*j*_{1},*j*_{2}= 1, 2, …,*n*._{s}

- (B.3) For
*i*= 1, 2, …,*n*, the_{p}*i*th element of**D**is ${y}^{p}({\mathbf{x}}_{i}^{p})-{\mathbf{f}}_{Z}^{\top}({\mathbf{x}}_{i}^{p},{\mathbf{c}}^{\u2605},\mathbf{t}){\mathit{\beta}}_{Z}-{\mathbf{f}}_{D}^{\top}({\mathbf{x}}_{i}^{p},{\mathbf{c}}^{\u2605},\mathbf{t}){\mathit{\beta}}_{D}$. For_{t}*i*=*n*+ 1,_{p}*n*+ 2, …,_{p}*n*+_{p}*n*, the_{s}*i*th element of**D**is ${y}^{s}({\mathbf{x}}_{i-{n}_{p}}^{s},{\mathbf{c}}_{i-{n}_{p}},{\mathbf{t}}_{i-{n}_{p}})-{\mathbf{f}}_{Z}^{\top}({\mathbf{x}}_{i-{n}_{p}}^{s},{\mathbf{c}}_{i-{n}_{p}},{\mathbf{t}}_{i-{n}_{p}}){\mathit{\beta}}_{Z}$_{t}

By conditioning on $(\stackrel{\uff3e}{{\mathbf{c}}_{\ell}^{\u2605}},\stackrel{\uff3e}{{\mathit{\varphi}}_{\ell}})$ and the law of large numbers,

$$\begin{array}{c}\text{Var}({Y}^{s}(\mathbf{x},\mathbf{c},\mathbf{t})|{\mathbf{y}}^{s},{\mathbf{y}}^{p},\stackrel{\uff3e}{{\mathbf{t}}^{\u2605}})\hfill \\ \hspace{1em}\approx \frac{1}{{N}_{MC}}{\displaystyle \sum _{\ell =1}^{{N}_{MC}}\text{Var}({Y}^{s}}(\mathbf{x},\mathbf{c},\mathbf{x})|{\mathbf{y}}^{s},{\mathbf{y}}^{p},\stackrel{\uff3e}{{\mathbf{t}}^{\u2605}},\stackrel{\uff3e}{{\mathbf{c}}_{\ell}^{\u2605}},\stackrel{\uff3e}{{\mathit{\varphi}}_{\ell}})+\frac{1}{{N}_{MC}}{\displaystyle \sum _{\ell =1}^{{N}_{MC}}}{[E({Y}^{s}(\mathbf{x},\mathbf{c},\mathbf{t})|{\mathbf{y}}^{s},{\mathbf{y}}^{p},\stackrel{\uff3e}{{\mathbf{t}}^{\u2605}},\stackrel{\uff3e}{{\mathbf{c}}_{\ell}^{\u2605}},\stackrel{\uff3e}{{\mathit{\varphi}}_{\ell}})]}^{2}-\phantom{\rule{thinmathspace}{0ex}}{\left[\frac{1}{{N}_{MC}}{\displaystyle \sum _{\ell =1}^{{N}_{MC}}E({Y}^{s}(\mathbf{x},\mathbf{c},\mathbf{t})|{\mathbf{y}}^{s},{\mathbf{y}}^{p},\stackrel{\uff3e}{{\mathbf{t}}^{\u2605}},\stackrel{\uff3e}{{c}_{\ell}^{\u2605}},\stackrel{\uff3e}{\mathit{\varphi}\ell})}\right]}^{2}\hfill \end{array}$$

(C.1)

and

$$\begin{array}{c}\text{Var}({Y}^{s}({\mathbf{x},\mathbf{c}}^{\u2605},\stackrel{\uff3e}{{\mathbf{t}}^{\u2605}}+\mathrm{\Delta}({\mathbf{x},\mathbf{c}}^{\u2605},\stackrel{\uff3e}{{\mathbf{t}}^{\u2605}}|{\mathbf{y}}^{s},{\mathbf{y}}^{p},\stackrel{\uff3e}{{\mathbf{t}}^{\u2605}})\hfill \\ \approx \frac{1}{{N}_{MC}}{\displaystyle \sum _{\ell =1}^{{N}_{MC}}\text{Var}({Y}^{s}}(\mathbf{x},\stackrel{\uff3e}{{\mathbf{c}}_{\ell}^{\u2605}},\stackrel{\uff3e}{{\mathbf{t}}^{\u2605}})+\mathrm{\Delta}(\mathbf{x},\stackrel{\uff3e}{{\mathbf{c}}_{\ell}^{\u2605}},\stackrel{\uff3e}{{\mathbf{t}}^{\u2605}})|{\mathbf{y}}^{s},{\mathbf{y}}^{p},\stackrel{\uff3e}{{\mathbf{t}}^{\u2605}},\stackrel{\uff3e}{{\mathbf{c}}_{\ell}^{\u2605}},\stackrel{\uff3e}{{\mathit{\varphi}}_{\ell}})+\frac{1}{{N}_{MC}}{\displaystyle \sum _{l=1}^{{N}_{MC}}}{\left[E({Y}^{s}(\mathbf{x},\stackrel{\uff3e}{{\mathbf{c}}_{l}^{\u2605}},\stackrel{\uff3e}{{\mathbf{t}}^{\u2605}})+\mathrm{\Delta}(\mathbf{x},\stackrel{\uff3e}{{\mathbf{c}}_{l}^{\u2605}},\stackrel{\uff3e}{{\mathbf{t}}^{\u2605}})|{\mathbf{y}}^{s},{\mathbf{y}}^{p},\stackrel{\uff3e}{{\mathbf{t}}^{\u2605}},\stackrel{\uff3e}{{\mathbf{c}}_{l}^{\u2605}},\stackrel{\uff3e}{{\mathit{\varphi}}_{l}})\right]}^{2}-{\left[\frac{1}{{N}_{MC}}{\displaystyle \sum _{l=1}^{{N}_{MC}}}E({Y}^{s}(\mathbf{x},\stackrel{\uff3e}{{\mathbf{c}}_{l}^{\u2605}},\stackrel{\uff3e}{{\mathbf{t}}^{\u2605}})+\mathrm{\Delta}(\mathbf{x},\stackrel{\uff3e}{{\mathbf{c}}_{l}^{\u2605}},\stackrel{\uff3e}{{\mathbf{t}}^{\u2605}})|{\mathbf{y}}^{s},{\mathbf{y}}^{p},\stackrel{\uff3e}{{\mathbf{t}}^{\u2605}},\stackrel{\uff3e}{{\mathbf{c}}_{l}^{\u2605}},\stackrel{\uff3e}{{\mathit{\varphi}}_{l}})\right]}^{2}\hfill \end{array}$$

(C.2)

for large *N _{MC}*. Using Equations (1), (2), (4), (5), (6), and the following formulas for $E({Y}^{s}(\mathbf{x},\mathbf{c},\mathbf{t})|{\mathbf{y}}^{s},{\mathbf{y}}^{p},\stackrel{\uff3e}{\mathbf{t}},\stackrel{\uff3e}{{\mathbf{c}}^{\u2605}},\stackrel{\uff3e}{\mathit{\varphi}}),E({Y}^{s}(\mathbf{x},{\mathbf{c}}^{\u2605},\mathbf{t})+\mathrm{\Delta}(\mathbf{x},{\mathbf{c}}^{\u2605},\mathbf{t})|{\mathbf{y}}^{s},{\mathbf{y}}^{p},\stackrel{\uff3e}{\mathbf{t}},\stackrel{\uff3e}{{\mathbf{c}}^{\u2605}},\stackrel{\uff3e}{\mathit{\varphi}})$, $\text{Var}({Y}^{s}(\mathbf{x},\mathbf{c},\mathbf{t})|{\mathbf{y}}^{s},{\mathbf{y}}^{p},\stackrel{\uff3e}{\mathbf{t}},\stackrel{\uff3e}{{\mathbf{c}}^{\u2605}},\stackrel{\uff3e}{\mathit{\varphi}})$, and
$\text{Var}({Y}^{s}(\mathbf{x},{\mathbf{c}}^{\u2605},\mathbf{t})+\mathrm{\Delta}(\mathbf{x},{\mathbf{c}}^{\u2605},\mathbf{t})|{\mathbf{y}}^{s},{\mathbf{y}}^{p},\stackrel{\uff3e}{\mathbf{t}},\stackrel{\uff3e}{{\mathbf{c}}^{\u2605}},\stackrel{\uff3e}{\mathit{\varphi}})$ one can compute Equations (18), ((19)), ((C.1)), and ((C.2)) by taking $(\stackrel{\uff3e}{\mathbf{t}},\stackrel{\uff3e}{{\mathbf{c}}^{\u2605}},\stackrel{\uff3e}{\mathit{\varphi}})$ to be $(\stackrel{\uff3e}{{\mathbf{t}}^{\u2605}},\stackrel{\uff3e}{{\mathbf{c}}_{\ell}^{\u2605}},\stackrel{\uff3e}{{\mathit{\varphi}}_{\ell}})$.

- (C.1) The conditional expectation $E({Y}^{s}(\mathbf{x},\mathbf{c},\mathbf{t})|{\mathbf{y}}^{s},{\mathbf{y}}^{p},\stackrel{\uff3e}{\mathbf{t}},\stackrel{\uff3e}{{\mathbf{c}}^{\u2605}},\stackrel{\uff3e}{\mathit{\varphi}}),\text{is}\phantom{\rule{thinmathspace}{0ex}}{\mathbf{f}}_{Z}^{\top}(\mathbf{x},\mathbf{c},\mathbf{t}){\mathit{\beta}}_{Z}+{\stackrel{\uff3e}{\mathrm{\Sigma}}}_{0Z}^{\top}{\stackrel{\uff3e}{\mathrm{\Sigma}}}_{ZZ}^{-1}{\stackrel{\uff3e}{\mathbf{D}}}_{t}$ where
_{0Z}is the (*n*+_{p}*n*) × 1 vector of covariances between_{s}*Z*(**x, c, t**) and the entries of (**Y**^{p┬},**Y**^{s┬})^{┬}given $(\stackrel{\uff3e}{\mathbf{t}},\stackrel{\uff3e}{{\mathbf{c}}^{\u2605}},\stackrel{\uff3e}{\mathit{\varphi}})$. For*i*= 1, …,*n*, the_{p}*i*th element of**Σ**_{0Z}is $\text{Cov}Z(\mathbf{x},\mathbf{c},\mathbf{t}),\phantom{\rule{thinmathspace}{0ex}}Z({\mathbf{x}}_{i}^{p},{\mathbf{c}}^{\u2605},\mathbf{t})|\stackrel{\uff3e}{\mathbf{t}},\stackrel{\uff3e}{{\mathbf{c}}^{\u2605}},\stackrel{\uff3e}{\mathit{\varphi}})$ while for*i*=*n*+ 1, …,_{p}*n*+_{p}*n*the_{s}*i*th element is $\text{Cov}(Z(\mathbf{x},\mathbf{c},\mathbf{t}),\phantom{\rule{thinmathspace}{0ex}}Z({\mathbf{x}}_{i-{n}_{p}}^{s},{\mathbf{c}}_{i-{n}_{p}},{\mathbf{t}}_{i-{n}_{p}})|\stackrel{\uff3e}{\mathbf{t}},\stackrel{\uff3e}{{\mathbf{c}}^{\u2605}},\stackrel{\uff3e}{\mathit{\varphi}})$. Similarly,_{ZZ}is Equation (B.1) and_{t}is given by (B.3), both conditional on $(\stackrel{\uff3e}{\mathbf{t}},\stackrel{\uff3e}{{\mathbf{c}}^{\u2605}},\stackrel{\uff3e}{\mathit{\varphi}})$. - (C.2) The conditional variance $\text{Var}({Y}^{s}(\mathbf{x},\mathbf{c},\mathbf{t})|{\mathbf{y}}^{s},{\mathbf{y}}^{p},\stackrel{\uff3e}{\mathbf{t}},\stackrel{\uff3e}{{\mathbf{c}}^{\u2605}},\stackrel{\uff3e}{\mathit{\varphi}})={\stackrel{\uff3e}{\sigma}}_{Z}^{2}-{\stackrel{\uff3e}{\mathbf{\Sigma}}}_{0Z}^{\top}{\stackrel{\uff3e}{\mathbf{\Sigma}}}_{ZZ}^{-1}{\stackrel{\uff3e}{\mathbf{\Sigma}}}_{0Z}$ where
_{0Z},_{ZZ}are defined in (C.1) while ${\stackrel{\uff3e}{\sigma}}_{Z}^{2}$ is taken from . - (C.3) The conditional mean of η(
**x, c**^{}) is $E({Y}^{s}(\mathbf{x},\mathbf{c},\mathbf{t})|{\mathbf{y}}^{s},{\mathbf{y}}^{p},\stackrel{\uff3e}{\mathbf{t}},\stackrel{\uff3e}{{\mathbf{c}}^{\u2605},}\stackrel{\uff3e}{\mathit{\varphi}})+E(\mathrm{\Delta}(\mathbf{x},\mathbf{c},\mathbf{t})|{\mathbf{y}}^{s},{\mathbf{y}}^{p},\stackrel{\uff3e}{\mathbf{t}},\stackrel{\uff3e}{{\mathbf{c}}^{\u2605},}\stackrel{\uff3e}{\mathit{\varphi}})$, where $E({Y}^{s}(\mathbf{x},\mathbf{c},\mathbf{t})|{\mathbf{y}}^{s},{\mathbf{y}}^{p},\stackrel{\uff3e}{\mathbf{t}},\stackrel{\uff3e}{{\mathbf{c}}^{\u2605},}\stackrel{\uff3e}{\mathit{\varphi}})$ and $E(\mathrm{\Delta}(\mathbf{x},\mathbf{c},\mathbf{t})|{\mathbf{y}}^{s},{\mathbf{y}}^{p},\stackrel{\uff3e}{\mathbf{t}},\stackrel{\uff3e}{{\mathbf{c}}^{\u2605},}\stackrel{\uff3e}{\mathit{\varphi}})$ is derived in (C.1) and Appendix B, respectively. - (C.4) The conditional variance $\text{Var}({Y}^{s}(\mathbf{x},\mathbf{c},\mathbf{t})+\mathrm{\Delta}(\mathbf{x},\mathbf{c},\mathbf{t})|{\mathbf{y}}^{s},{\mathbf{y}}^{p},\stackrel{\uff3e}{\mathbf{t}},\stackrel{\uff3e}{{\mathbf{c}}^{\u2605}},\stackrel{\uff3e}{\mathit{\varphi}})={\stackrel{\uff3e}{\sigma}}_{Z}^{2}+{\stackrel{\uff3e}{\sigma}}_{D}^{2}-({\stackrel{\uff3e}{\mathrm{\Sigma}}}_{0Z}^{\top}+{\stackrel{\uff3e}{\mathrm{\Sigma}}}_{DZ}^{\top}){\stackrel{\uff3e}{\mathrm{\Sigma}}}_{ZZ}^{-1}({\stackrel{\uff3e}{\mathrm{\Sigma}}}_{0Z}+{\stackrel{\uff3e}{\mathrm{\Sigma}}}_{DZ}^{\top})$, where ${\stackrel{\uff3e}{\sigma}}_{Z}^{2}$ and ${\stackrel{\uff3e}{\sigma}}_{D}^{2}$ are taken from ,
_{0Z}and_{ZZ}given in (C.1), and_{DZ}is given by (B.1) where conditioning is given $(\stackrel{\uff3e}{\mathbf{t}},\stackrel{\uff3e}{{\mathbf{c}}^{\u2605}},\stackrel{\uff3e}{\mathit{\varphi}})$.

SUPPLEMENTAL MATERIALS

**A MATLAB Routine:** The supplemental materials include a user guide to this MATLAB package (STaCguide.pdf), the software package (directory STaCprogram), and an example subdirectory having the driver file for running the example in Section 5.3. (STaC_software_version1.0.zip, zip archive)

Gang Han, H. Lee Moffitt Cancer Center & Research Institute, MRC/BIOSTAT, 12902 Magnolia Drive, Tampa, FL 33612, (Email: gro.ttiffom@nah.gnag)

Thomas J. Santner, Department of Statistics, The Ohio State University, 415 Cockins Hall, 1958 Neil Avenue, Columbus, OH 43210, (Email: ude.uso.tats@sjt)

Jeremy J. Rawlinson, C2-527 Clinical Programs Center, Department of Clinical Sciences, College of Veterinary Medicine, Cornell University, Ithaca, NY 14853, (Email: ude.llenroc@4rjj)

- Bayarri MJ, Berger JO, Cafeo J, Garcia-Donato G, Liu F, Palomo J, Parthasarathy RJ, Paulo R, Sacks J, Walsh D. Computer Model Validation With Functional Output. The Annals of Statistics. 2007a;35(5):1874–1906.
- Bayarri MJ, Berger JO, Paulo R, Sacks J, Cafeo JA, Cavendish J, Lin C-H, Tu J. A Framework for Validation of Computer Models. Technometrics. 2007b;49(2):138–154.
- Cox DD, Park JS, Singer CE. A Statistical Method for Tuning a Computer Code to a Database. Rice University, Dept. of Statistics; 1996. Technical Report 96-3.
- Craig PC, Goldstein M, Rougier JC, Seheult AH. Bayesian Forecasting for Complex Systems Using Computer Simulators. Journal of the American Statistical Association. 2001;96:717–729.
- Craig PS, Goldstein M, Seheult AH, Smith JA. Bayes Linear Strategies for History Matching of Hydrocarbon Reservoirs. In: Bernardo JM, Berger JO, Dawid AP, Smith AFM, editors. Bayesian Statistics. Vol. 5. Oxford: Oxford University Press; 1996. pp. 69–95.
- Fang KT, Li R, Sudjianto A. Design and Modeling for Computer Experiments. London: Chapman & Hall; 2005.
- Forrester A, Sobester A, Keane A. Engineering Design via Surrogate Modelling: A Practical Guide. Chichester, U.K: Wiley; 2008.
- Gattiker J. Using the Gaussian Process Model for Simulation Analysis (GPM/SA) Code. Los Alamos National Laboratory; 2005. technical report.
- Higdon D, Gattiker J, Williams B, Rightley M. Computer Model Calibration Using High Dimensional Output. Journal of the American Statistical Association. 2008;103:570–583.
- Higdon D, Kennedy M, Cavendish J, Cafeo J, Ryne R. Combining Field Data and Computer Simulations for Calibration and Prediction. SIAM Journal of Scientific Computing. 2004;26:448–466.
- Higdon D, Williams B, Moore L, McKay M, Keller-McNulty S. Uncertainty Quantification for Combining Experimental Data and Computer Simulations. Los Alamos National Laboratory; 2005. Technical Report LA-UR 05-4951.
- Kennedy MC, O’Hagan A. Bayesian Calibration of Computer Models. Journal of the Royal Statistical Society, Ser. B. 2001;63:425–464. (with discussion)
- Linkletter C, Bingham D, Hengartner N, Higdon D, Ye KQ. Variable Selection for Gaussian Process Models in Computer Experiments. Technometrics. 2006;48:478–490.
- Loeppky J, Bingham D, Welch W. Computer Model Calibration or Tuning in Practice. University of British Columbia; 2006. technical report.
- McKay MD, Beckman RJ, Conover WJ. A Comparison of Three Methods for Selecting Values of Input Variables in the Analysis of Output From a Computer Code. Technometrics. 1979;21:239–245.
- Park JS. Ph.D. thesis. Champaign/Urbana, IL: University of Illinois; 1991. Tuning Complex Computer Codes to Data and Optimal Designs.
- Qian PZG. Nested Latin Hypercube Designs. Biometrika. 2009 to appear, doi:
*10.1093/biomet/asp045*. - Qian PZG, Wu CFJ. Sliced Space-Filling Designs. Biometrika. 2009 to appear, doi:
*10.1093/biomet/asp044*. - Qian PZG, Ai M, Wu CFJ. Construction of Nested Space-Filling Designs. The Annals of Statistics. 2009;37:3616–3643.
- Rawlinson JJ, Furman BD, Li S, Wright TM, Bartel DL. Retrieval, Experimental, and Computational Assessment of the Performance of Total Knee Replacements. Journal of Orthopaedic Research. 2006;24(7):1384–1394. [PubMed]
- Saltelli A, Chan K, Scott E. Sensitivity Analysis. Chichester: Wiley; 2000.
- Santner TJ, Williams BJ, Notz WI. The Design and Analysis of Computer Experiments. New York: Springer Verlag; 2003.

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |