Evan Baker
email: e.baker@exeter.ac.uk

Pierre Barbillon

Arindam Fadikar

Robert B Gramacy

Radu Herbei

David Higdon

Jiangeng Huang

Leah R Johnson

Pulong Ma

Anirban Mondal

Bianica Pires

Jerome Sacks

Vadim Sokolov

Analyzing Stochastic Computer Models: A Review with Opportunities

Keywords: Computer Model, Gaussian Process, Uncertainty Quantification, Emulator, Computer Experiment, Agent Based Model, Surrogates, Calibration

In modern science, computer models are often used to understand complex phenomena, and a thriving statistical community has grown around analyzing them. This review aims to bring a spotlight to the growing prevalence of stochastic computer models -providing a catalogue of statistical methods for practitioners, an introductory view for statisticians (whether familiar with deterministic computer models or not), and an emphasis on open questions of relevance to practitioners and statisticians. Gaussian process surrogate models take center stage in this review, and these, along with several extensions needed for stochastic settings, are explained. The basic issues of designing a stochastic computer experiment and calibrating a stochastic computer model are prominent in the discussion. Instructive examples, with data and code, are used to describe the implementation of, and results from, various methods.

Introduction

Computer models, also known as simulators, are in use everywhere. These are programs which describe and approximate a process of interest. The code typically takes a set of inputs and produces some output. Stochastic simulators, unlike deterministic ones, can produce different output with the same inputs due to the presence of random elements.1 Such computer models are in wide use. For example, agent-based models (ABMs) deal with large populations of individuals, where specific actions taken at any one time-step have complexities and uncertainties that do not allow deterministic modeling. ABMs are prevalent [START_REF] Johnson | Implications of dispersal and life history strategies for the persistence of linyphiid spider populations[END_REF][START_REF] Johnson | Parameter inference for an individual based model of chytridiomycosis in frogs[END_REF][START_REF] Ramsey | Management of bovine tuberculosis in brushtail possums in new zealand: predictions from a spatially explicit, individual-based model[END_REF][START_REF] Smieszek | Reconstructing the 2003/2004 h3n2 influenza epidemic in switzerland with a spatially explicit, individual-based model[END_REF][START_REF] Grimm | A standard protocol for describing individual-based and agent-based models[END_REF] and used to explore complex phenomena in sociology, transportation, ecology, epidemiology, and other phenomena.

The following is a basic model of a stochastic simulator experiment. If the code is run at a (vector) input x producing a (scalar) output y(x), this could be represented as:

y(x) = M (x) + v, v ∼ N (0, σ 2 v (x)), (1.1)
where M (x) is the expected value, E[y(x)], of the output. The variability v accounts for the randomness of the stochastic simulator, ultimately caused by pseudo-random number generation within the code. Its variance, σ 2 v , often depends on x, with constant variance subsumed as a special case. For deterministic simulators, σ 2 v = 0. Randomness in stochastic simulators invariably requires many simulations thereby limiting the complexity (including the size of the input dimension) that can be effectively treated. The prospect of replicate runs in stochastic simulators introduces a trade-off between replication and exploration, a challenging design issue. The noise, v, makes additional demands on the analysis when its variance is non-constant. This article examines these basic issues, identifies accessible and effective methods, and points to unresolved questions that should be addressed.

Equation 1.1 is often used to model physical experiments, where an observation y(x) is truth, M (x), plus measurement error (and, possibly, intrinsic variability as well) or, for an observational study, where M (x) is fit to the observations with residual variance. Because they are structurally the same, physical experiments can be analyzed with methods used for stochastic simulators [START_REF] Gao | Predicting urban ozone levels and trends with semiparametric modeling[END_REF]. However, the contexts and goals are often different, leading to different problem formulations and different interpretations of results.

The choice of method, with its assumptions and limitations, is crucial for any analysis of an experiment. An inclination for simplicity and availability of software would encourage the use of a standard statistical regression model (for example, linear regression) for M with a constant σ 2 v . That this approach is effective under some circumstances, especially when the space, X, of possible inputs is small, begs the question of how reliable it can be as a general prescription. Complex systems modeled by a simulator may neither suggest nor allow much simplification. The methods described in this review allow the simulated data to guide the choice of method under general conditions with little, or no simplification. Statistics [START_REF] Sacks | Design and analysis of computer experiments[END_REF][START_REF] Kennedy | Bayesian calibration of computer models[END_REF] and Applied Mathematics [START_REF] Sullivan | Introduction to Uncertainty Quantification[END_REF] play prominent roles in the design and analysis of deterministic computer experiments. Unsurprisingly, some methods developed for deterministic simulators have modifications that can be used in the stochastic context. Alternatives, driven by the stochasticity, are necessary in many contexts. These structural differences will be noted in the narrative below.

Goals

We have three primary goals; all related to the cross-disciplinary nature of this topic.

One goal is to bring effective statistical methods to the attention of subject scientists and enable a deeper understanding of stochastic simulators in use. The descriptions below of statistical tools used (or cited) try to avoid being bogged down in mathematical intricacies. Some details of individual methods are included to help in understanding the strengths and weaknesses of the methods. Application of a number of methods is exemplified on testbed cases (Section 2), and available software for methods are identified where possible.

A second goal is to familiarize statisticians with an area of major importance that is crucial to the formation of evidence-based policy. Statisticians are sorely needed in the study and application of agent-based models (ABMs) and stochastic simulators in general. Researchers familiar with deterministic simulation techniques will see immediate opportunities, but statistical expertise of all kinds is essential to advance the study of stochastic simulators.

The analysis of stochastic simulators is a developing field with many unsolved problems. Challenges are often driven by the scale of the problems and a range of issues whose resolution requires close cooperation between statisticians, subject scientists, and computer scientists. A third goal of this paper is to spur that process.

The review is structured as follows: Section 3 describes the models that form the basis for the analyses; Section 4 is devoted to the fundamental question of what simulator runs to make. Section 5 addresses a common objective of simulation experiments: calibration. Section 6 discusses other models and objectives that are important, but are more on the "boundaries" of this review and are therefore less detailed. Finally, Section 7 summarizes conclusions and poses unanswered questions. The references here do not cover the entire body of work on stochastic simulators but, together with this overview, should provide adequate coverage of the problems discussed.

Example Simulators

Three stochastic simulators will be discussed throughout this review to aid understanding. Two are deliberately simplified and used to exhibit key features of the methods presented. In some cases simpler strategies could be equally effective because the complexity of the models has been greatly reduced. Since the data/generating mechanisms used are available, others can compare different strategies, but the demonstration purpose is the one that is relevant in the discussion and reported computations. The third is a model which we use to anchor and motivate methods. The specific model in question is an epidemiological model developed in response to the Ebola epidemic of 2014. For the Ebola model, a synthetic population representing the individuals in Liberia (population ∼ 4.5 million) and their activity schedules, inducing a time-varying contact network of individuals and locations, was developed [START_REF] Mortveit | Synthetic populations and interaction networks for guinea. liberia and sierra leone[END_REF], and paired with an agent-based model [START_REF] Bisset | Epifast: a fast algorithm for large scale realistic epidemic simulations on distributed memory systems[END_REF]. Together, this ABM models a contagion spreading from one individual to another in Liberia. Since the parameter for contagion, transmissibility, only controls the probability of infection given an interaction occurs, this model (and many like it) is stochastic. The model is updated daily, with the progress of the disease determined by the activity schedule, contact details, and other epidemiological characteristics. This model is complex, with high dimensional outputs, multiple unknown inputs, and non-normality all present. The analysis performed by [START_REF] Fadikar | Calibrating a stochastic, agent-based model using quantile-based emulation[END_REF] tackles all of these using ideas discussed within this article (see Sections 3.3.1, 3.4, and 5 ).

Fish Capture-Recapture

The first simplified stochastic simulator we consider mimics the movements and schooling behavior of fish in a mark-recapture application. Mark and recapture involves capturing a sample of the population, marking and releasing them, and following up by capturing another sample and counting how many are marked -the recaptured. The number recaptured allows estimation of the population size [START_REF] Begon | Investigating Animal Abundance: Capture-Recapture for Biologists[END_REF]. The process is modeled by initializing a population of fish at random locations in a 2-d, rectangular lake with boundary conditions. The fish begin moving and schooling according to simple, agent-based rules. After an initial period of time, 100 fish are marked as they pass through a "net" in the lake. After a second period of time, 100 fish are captured using the same net and the number of "recaptured" are recorded.

This agent-based model is a modified version of the flocking model developed in NetLogo [START_REF] Wilensky | NetLogo[END_REF]. The collective behavior that emerges in the flocking model is the result of providing each individual agent with the same set of simple rules [START_REF] Reynolds | Flocks, herds and schools: A distributed behavioral model[END_REF]. The flocking model is modified to include the mark-recapture dynamics described above. Given an observed count of recaptured fish, this model can be used to estimate the total size of the fish population (see Section 5.3). The only input considered is the number of fish in the total population and the output is the number of recaptured fish. Other inputs for this model control the individual movement rules of the fish; for simplicity these are ignored here and set to default values. Supplementary code, and compiled Rmarkdown documents, corresponding to our analysis of this simulator can be found at https://github.com/jhuang672/fish. Running the simulator afresh will require the installation of NetLogo from https://ccl.northwestern.edu/netlogo/.

Ocean Circulation

The second simplified example is a stochastic simulator that models the concentration of oxygen in a thin water layer (around 2000m deep) in the South Atlantic ocean [START_REF] Mckeague | Statistical inversion of south atlantic circulation in an abyssal neutral density layer[END_REF][START_REF] Herbei | Estimating ocean circulation: An MCMC approach with approximated likelihoods via the bernoulli factory[END_REF]. The physical model is described via an advection-diffusion equation (equation ( 4) of [START_REF] Mckeague | Statistical inversion of south atlantic circulation in an abyssal neutral density layer[END_REF]), i.e., a non-linear partial differential equation (PDE) describing the dynamics of oxygen concentration in terms of the water velocities and diffusion coefficients. For a given set of inputs, the solution of the advection-diffusion equation is not available in closed form. However, using theoretical results [START_REF] Feynman | Space-time approach to non-relativistic quantum mechanics[END_REF][START_REF] Kac | On distributions of certain Wiener functionals[END_REF], the solution can be closely approximated through an associated random process [START_REF] Herbei | Estimating ocean circulation: An MCMC approach with approximated likelihoods via the bernoulli factory[END_REF]. For a specific location within the domain, random paths of the process are generated, producing noisy outcomes that approximate the solution to the PDE at that location. This example is simplified by taking the oxygen concentration output to only depend on four inputs: two unknown diffusion constants (K x and K y ) and the two location variables (latitude and longitude). All other inputs are held fixed at nominal values. Such stochastic approximations are numerous in physical sciences, either due to computational limitations, a lack of complete understanding of the underlying system, or because the system under study is itself believed to be random. Supplementary code, and compiled Rmarkdown documents, corresponding to our analysis of this simulator can be found at https://github.com/Demiperimetre/Ocean.

Statistical Models

An experiment of running a simulator and producing data whose output is described by equation 1.1 can have a multitude of goals. A principal objective, and the one we focus on here, is using the simulated data to predict values of the simulator, M (x) + v, and the uncertainties of the predictions, at untried xs in a context where getting new runs of the simulator is not cost-free. When M is believed to be "simple" (for example, a polynomial function of the coordinates of x) there are many standard "classical" techniques that can be used to approximate M . For example, linear regression models and generalised linear models have been used by [START_REF] Andrianakis | Efficient history matching of a high dimensional individual-based HIV transmission model[END_REF] and [START_REF] Marrel | Global sensitivity analysis of stochastic computer models with joint metamodels[END_REF]. Complex problems such as those in Section 2 are less easily managed: specifying a functional form for complex M requires sufficient prior knowledge or a huge abundance of data, both of which are often lacking. A prime emphasis of this article is on methods that have been developed to cope with such concerns; adequate references for a variety of standard methods are available for simpler circumstances.

There are a range of factors that need to be taken into account before choosing a statistical model (hereon referred to as a surrogate model, as it acts as a surrogate for the computer model). In addition to methodological assumptions, it is important to consider the "context", that is, the conditions of the particular problem being studied, leading to equation 1.1 and its extensions. Some important contexts include:

• The space of inputs is usually a hyper-rectangle: each coordinate of an input x is constrained by upper and lower bounds. Section 2.2 simplifies issues by taking a rectangular input space even though the Atlantic Ocean is not rectangular.

• The output y in equation 1.1 is scalar, but multiple output, such as time-series, is also common.

• Some inputs may be categorical rather than numerical.

• The probability distribution of v, the variability, is often taken to be normal, but often invalid, as with the Ebola model.

Stretching back to [START_REF] Sacks | Design and analysis of computer experiments[END_REF]; [START_REF] Currin | Bayesian prediction of deterministic functions, with applications to the design and analysis of computer experiments[END_REF], a vast literature, mostly on deterministic simulators, has found that a Gaussian Process (GP) model produces a flexible, effective surrogate for M . This approach, and its modifications needed to address the presence of input dependent σ 2 v (x) in equation 1.1, can be effective for stochastic simulation as has been documented in the literature (e.g. [START_REF] Kleijnen | Kriging metamodeling in simulation: A review[END_REF][START_REF] Kleijnen | Regression and kriging metamodels with their experimental designs in simulation: a review[END_REF] and will be apparent below. A thorough intuitive explanation (for deterministic computer models) can be found in [START_REF] O'hagan | Bayesian analysis of computer code outputs: A tutorial[END_REF]. More technical descriptions of GPs from a statistical perspective can be found in [START_REF] Santner | The Design and Analysis of Computer Experiments[END_REF] and [START_REF] Gramacy | Surrogates: Gaussian Process Modeling, Design and Optimization for the Applied Sciences[END_REF]; for a machine learning perspective, see [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF]. In brief, the use of GPs allows computer model runs to play the key role in selecting a surrogate and assessments of its uncertainty in prediction. Deep learning methods, such as neural networks, and other general-purpose predictors are also in wide use. These modern learning machines have difficulties in producing uncertainties and identifying critical inputs but there is active research directed towards that end [START_REF] Neal | Bayesian Learning for Neural Networks[END_REF][START_REF] Graves | Practical variational inference for neural networks[END_REF][START_REF] Welling | Bayesian learning via stochastic gradient Langevin dynamics[END_REF][START_REF] Papamakarios | Normalizing flows for probabilistic modeling and inference[END_REF][START_REF] Gal | Dropout as a Bayesian approximation: representing model uncertainty in deep learning[END_REF][START_REF] Lakshminarayanan | Simple and scalable predictive uncertainty estimation using deep ensembles[END_REF].

Gaussian Process Surrogates

Suppose that the input space X is a hyper-rectangle in d-dimensions; the output y(x) is univariate (scalar); and that variability is normally distributed. Additionally, assume:

A1 The variability, v, has constant variance σ 2 v A2 The mean M (x) = µ + Z(x) A3 µ is constant A4 Z(•)
is a Gaussian Process on X with mean 0 and covariance function K, deconstructed as a product of a variance σ 2 Z and a correlation function C. The technical definition of a GP (Assumption A4) is: for any finite N and collection of inputs X N = (x 1 , . . . , x N ), Z N = (Z(x 1 ), . . . , Z(x N )) is a multivariate normal random variable with mean 0 and N × N covariance matrix K N , whose entries are K(x i , x j ). It follows that the simulator output, Y N = (y(x 1 ), . . . , y(x N )) , is also multivariate normal but with mean µ1 and covariance matrix K N + σ 2 v I N , where I N is the identity N × N matrix and 1 is the N -vector of 1s.

One interpretation is that these assumptions describe a prior distribution on all possible functions for the mean M . Different choices for the GP allow for different classes of possible M ; the power of a GP is that these classes can be big enough to allow for all reasonable possibilities. After specifying µ, K, and σ 2 v , a Bayesian analysis can then be carried out, resulting in a posterior distribution for all the functions that can still represent M after accounting for the observed simulator runs.

Another interpretation of M and Assumptions A2 and A4 is to think of M as a random function, with µ being a regression function (as in linear regression), and the GP for Z modeling the deviation from µ. Both formulations have the same mathematical structure but with differing interpretations.

The predictive distribution for any new run, y(x new ), given the observed simulator data {X N , Y N } is also normal, and has a known analytical form. The mean µ N (x new ) and variance σ 2 N (x new ) of predictions are:

µ N (x new ) = µ + k N (x new ) (K N + σ 2 v I N ) -1 (Y N -µ1) (3.1) σ 2 N (x new ) = σ 2 v + σ 2 Z -(k N (x new ) (K N + σ 2 v I N ) -1 k N (x new ), (3.2)
with k N (x new ) denoting the N -vector (K(x new , x 1 ), . . . , K(x new , x N )) of covariances between the desired prediction and observed data. Once the correlation function C is specified, parameters (µ, σ 2 Z , and σ 2 v ) can be estimated from the data. For specifying C, the approach taken for deterministic simulators can be adopted here: specify a parameterized family C θ and use the data to estimate θ, thereby tailoring C to observations. One example for C θ is the family of squared-exponential correlation functions (also known as the Gaussian kernel):

C θ (x, w) = exp - d j=1 (x j -w j ) 2 θ j . (3.3)
This correlation function is suited for approximating very smooth, infinitely differentiable, functions over dimension d. Alternative correlation functions exist and are used; one commonly used alternative is the Matérn 5/2 correlation function [START_REF] Stein | Interpolation of Spatial data: Some Theory for Kriging[END_REF], which is appropriate for approximating less-smooth functions (only 2 derivatives). 2

2 The Matérn 5/2 correlation function has the form 1

+ √ 5(xj -wj ) ρ + 5(xj -wj ) 2 3ρ 2 exp - √ 5(xj -wj ) ρ
. Further discussion of the features of different kernels can be found in Chapter 4.2 of [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF], Chapter 2.2 of [START_REF] Santner | The Design and Analysis of Computer Experiments[END_REF], or Chapter 5.3 of [START_REF] Gramacy | Surrogates: Gaussian Process Modeling, Design and Optimization for the Applied Sciences[END_REF].

With a choice of the family C θ and Assumptions A1-A4, the likelihood of the observed output is available and maximum likelihood estimates (MLEs) μ, σ2 v , σ2 Z , and θ can be calculated. Henceforth µ N (x new ) and σ 2 N (x new ) will be used to denote the mean and variance of the predictive distribution even when the parameters in equations 3.1 and 3.2 are estimated. The predictive probability distribution for the computer model output y(x new ) is then:

y(x new ) ∼ N (µ N (x new ), σ 2 N (x new )). (3.4)
A proper assessment of uncertainty is lost by plugging-in estimated parameters without accounting for their uncertainty. Accordingly, the predictive variance, σ 2 N (x new ), obtained this way is called the plugin (or nominal ) predictive variance. The alternative of a full Bayesian analysis to estimate the parameters can be computationally impractical in many circumstances, though not impossible (intermediate schemes and approximations have proven to be useful, e.g., [START_REF] Spiller | Automating emulator construction for geophysical hazard maps[END_REF]).

For the correlation function in equation 3.3, and for others such as the Matérn 5/2, the correlation between Z(x) and Z(w) depends only on x -w, the difference between the two vectors of inputs. That is, Z is assumed to be a stationary GP (and, consequently, so is y). For functions exhibiting markedly different behavior in one region of input space than in another part, stationarity is problematic. This issue is tackled and discussed in [START_REF] Gramacy | Bayesian treed Gaussian process models with an application to computer modeling[END_REF], [START_REF] Ba | Composite Gaussian process models for emulating expensive functions[END_REF], [START_REF] Kersaudy | A new surrogate modeling technique combining kriging and polynomial chaos expansionsapplication to uncertainty analysis in computational dosimetry[END_REF], and [START_REF] Chen | Analysis methods for computer experiments: how to assess and what counts[END_REF], among others, and Section 6.1 discusses one solution.

Despite the fairly complex mathematical expressions above, Gaussian processes are easily accessible thanks to numerous available packages (for example: DiceKriging in R [START_REF] Roustant | DiceKriging: Kriging Methods for Computer Experiments[END_REF], the hetGP R package [START_REF] Binois | hetGP: Heteroskedastic Gaussian Process Modeling and Design under Replication. R package version 1[END_REF] mentioned later, and the GaussianProcessRegressor function from scikit-learn in Python [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF]). In general, a GP is a flexible method for estimating the mean M (x) of the simulator output, despite lack of prior knowledge. This is illustrated in the top panels of Figures 1 and2, but we first introduce a vital modeling twist to cope with a common feature of stochastic computer simulations.

Heteroscedastic GP Surrogates

The constant variance Assumption (A1) simplifies the construction of a statistical model because only one intrinsic variance parameter σ 2 v needs to be estimated. When σ 2 v (x) is believed to vary over the input space more must be done. Boukouvalas et al. (2014a) model σ 2 v (x) as exp(h(x)) for simple functions h (e.g., polynomials), a simple extension to assuming just one variance parameter σ 2 v . (The exponential transform ensures positivity of the variance.) Like analogous approaches to predicting the mean (briefly discussed in Section 3.1), it isn't clear what to use for h, and its simplicity may not meet the complexities found in many applications.

GPs are used for σ 2 v by several authors [START_REF] Goldberg | Regression with inputdependent noise: a Gaussian process treatment[END_REF][START_REF] Kersting | Most likely heteroscedastic Gaussian process regression[END_REF][START_REF] Boukouvalas | Learning heteroscedastic Gaussian processes for complex datasets[END_REF][START_REF] Ankenman | Stochastic kriging for simulation metamodeling[END_REF]Binois et al., 2018a). The difficulty is that doing so directly depends on observing σ 2 v (x) at the inputs X n , but these values are not observed. If there are enough replicated simulation runs, r i , at the inputs x i , then the sample variances (s 2 (x i ) = 1 r i -1 r i j=1 (y(x ij )-ȳi ) 2 , for i = 1 . . . n) at the x i s can be used to estimate the σ 2 v (x) at the inputs X n . Equations 3.1 and 3.2 can then be used to predict σ 2 v (x new ). (Working with the logarithm of the sample variances and then exponentiating the results avoids negative predictions of the variance.) But this approach, called stochastic kriging (SK, [START_REF] Ankenman | Stochastic kriging for simulation metamodeling[END_REF], is limited by the need for adequate numbers of replicates at each input and the possible inefficiency of treating the variance and mean processes separately.

Those limitations can be removed by considering the intrinsic variances at the inputs, (σ 2 v (x 1 ), . . . , σ 2 v (x n )), as unknown parameters (a.k.a., latent variables) to be estimated in the same manner as all the other unknown parameters. [START_REF] Goldberg | Regression with inputdependent noise: a Gaussian process treatment[END_REF] do so in a fully Bayesian, but computationally taxing, way. Efforts to reduce these costs form the essence of approaches by [START_REF] Kersting | Most likely heteroscedastic Gaussian process regression[END_REF] and [START_REF] Boukouvalas | Learning heteroscedastic Gaussian processes for complex datasets[END_REF]. A recent variant, proposed in Binois et al. (2018a) along with accessible software hetGP [START_REF] Binois | hetGP: Heteroskedastic Gaussian Process Modeling and Design under Replication. R package version 1[END_REF], resolves the computational hazards and is the method described and used in this review.

The technical details addressing the computational barriers of a heteroscedastic GP (hetGP) have three elements. One, hetGP models the log variances as the mean output of a GP on latent (hidden) variables. The second uses Woodbury matrix identities [START_REF] Harville | Matrix algebra from a statistician's perspective[END_REF] to reduce computations from treating all N observations to computations involving only the n unique inputs, a reduction of computational complexity from O(N 3 ) to O(n 3 ), especially relevant when there are many replicates. The third element uses MLE to set all parameters.

While full details are provided by Binois et al. (2018a), some specifics of the first element of the description above are worth noting. With λ(x) = σ 2 v (x)/σ 2 Z and Λ n = (λ(x 1 ), . . . , λ(x n )) for the n distinct inputs, log Λ n is taken to be the predictive mean of a GP on latent (hidden) variables, ∆ n = (δ 1 , . . . , δ n ). For ease of exposition assume the GP has 0-mean (a constant mean is actually the default setting in hetGP) and take the covariance function for ∆ n to be σ 2 g (C g + gR -1 ) where g > 0, R = diag(r 1 , . . . , r n ), and C g is a correlation function with parameters θ g . Then log Λ n = C g (C g + gR -1 ) -1 ∆ n . This latent ∆ n approach facilitates smooth estimates of Λ n and provides a fixed functional form for λ(x), but does not incorporate the resulting uncertainty due to the estimates of the intrinsic σ 2 v (x) in predictions. Given Λ n , the Woodbury identities [START_REF] Harville | Matrix algebra from a statistician's perspective[END_REF] reduce the likelihood of Y N , the output at all inputs including replicates, to depend on quantities of size n. Maximum likelihood estimates for the unknown parameters can then be computed at a cost of O(n 3 ). Derivatives are also computable at a cost of O(n 3 ), further facilitating optimization for maximizing likelihood.

As a side note, heteroscedastic measurement error is sometimes present in spatial statistics models (which are often related to surrogate models); however we know of no such models which allow for the full modeling and predictions of the intrinsic variance process in the same way as a hetGP. For example, the model in [START_REF] Nguyen | Multivariate spatial data fusion for very large remote sensing datasets[END_REF] allows for non-constant measurement error at different sites, but it does not estimate these measurement errors jointly with the other model parameters, nor does it allow for the prediction of the measurement errors at new unseen sites. This is mostly because there is little interest in predicting the measurement error process in spatial statistics (the "true" underlying signal is the objective), whereas with stochastic simulators the intrinsic variability can be of direct modeling interest. Fish Example. We apply both an ordinary homoscedastic GP (homGP) and a hetGP surrogate to the fish example from Section 2.1. The simulation budget is constrained to 400 runs and focuses on the relationship between the total number, x, of fish in a population and the number, y(x), of fish recaptured in the second round of capture. The total population is an integer between 150 and 4000. The simulator is run 20 times at each of 20 unique x locations in [150,4000], chosen via a maximin Latin hypercube design (see Section 4). The number of fish counted cannot be less than zero, but the normality assumption would allow negative fish counts, so we square root the simulated output before performing our analysis, squaring the resulting predictions to return to the original scale afterwards. In addition we estimate "truth" by generating another data set; replicating 500 times at each of the same 20 sites.

Applying a homGP surrogate with squared exponential correlation function produces the results in the upper left panel of Figure 1; the upper right panel shows the results of hetGP. The predicted intervals for the fish model are obtained in the transformed (square-root) space, and squared to get back to the original space. 3 The lower panels are plots with the "true" 2.5%, 50%, and 97.5% quantiles superimposed.

The key conclusion is that both homGP and hetGP capture the non-linear trend (though a bit off in the region near 800). The presence of non-constant intrinsic variability is clear from the truth plot, with the region near 800 showing higher variability than elsewhere. The hetGP surrogate does not fully resolve the non-constant predictive variability, which includes both the intrinsic variability and that from the surrogate, but does improve on homGP. Full resolution is largely a matter of simulation budget though alternative designs may further improve hetGP. Our supplementary material includes improved results using the sequential design scheme of Section 4.3. The takeaway message is that the trend is readily treated by both homGP and hetGP; heteroscedasticity encourages use of hetGP perhaps with added simulations or improved designs.

Ocean Example. For the ocean model (Section 2.2), we take each simulation run to be the average of 6 simulation runs. The true simulator is known to be non-normal; this adjustment makes the example more Gaussian. For now, we fix the two diffusion coefficients, K x = 700 and K y = 200, leaving the two spatial coordinates as the only Homoscedastic Heteroscedastic [150,4000], rounded down to nearest integers). The left panels use homGP -the solid red line is the median of the predictive distribution and the dashed red-lines form the 95% uncertainty intervals. The right panels use hetGP. The upper panels include the data used to fit the surrogates; the lower panels omit the data but include the "true" values in black. varying inputs. Using 1000 simulations (50 sites each replicated 20 times), we obtain, for surrogates homGP and hetGP, the predictive mean surface and the predictive standard deviation surface (that is, the standard deviations for prediction of the simulator output, accounting for both the uncertainty around the predictive mean and the intrinsic variance estimate σ 2 v ). These surfaces are plotted in Figure 2, with the left column for homGP and the right column for hetGP.

The mean surfaces for both surrogates are similar. The predictive standard deviation for homGP (bottom-left) is relatively constant across the input area (clearly affected by the constraint that the intrinsic variance σ 2 v is constant). The standard deviation surface for hetGP is markedly different, evidence that intrinsic variance is non-constant. The "truth" is obtained using replicate runs (up to 100,000) of the simulator at 500 sites (chosen via a LHD, Section 4), averaging the replicates at each site to get the true mean and the square deviations from the mean to get the true variance). These are plotted in the appendix (Section A) and the supplementary material; they confirm the presence of non-constant intrinsic variance. Moreover, the standard Homoscedastic Heteroscedastic deviation plot for hetGP exhibits a structure similar to the truth plot, leading to the conclusion that hetGP is the better surrogate for this problem. However, this conclusion comes with a caution: repeating this experiment reveals a great deal of variability in the standard deviation plot, due to variability in the design and the simulations (discussed further in Section 4.3).

Overall, reliable predictions of the mean are achieved, but the uncertainties are less certain. This is similar to the the Fish example, and improving the uncertainties would require more simulation. These results point to the superiority of hetGP to homGP. This is confirmed via a numerical comparison in Section 4.3, where a sequential design is also examined and compared.

Non-Normal Variability

In many applications, assuming the variability v to be normally distributed is inappropriate. For example, count data, as in the Ebola or fish model, is non-normal and cannot be less than 0. Additionally, the distribution of v at a given input x may not be unimodal: in the Ebola example, even with the inputs x fixed, repeated simulations can lead to two distinct groups of possible infection counts, implying bimodality. In some simulators, there may be a greater tendency for extreme values (fatter tails) in the distribution of v. With these possibilities normality can be a strong assumption to be used with caution.

Transformation of the data is a time-honored device that sometimes induces "enough" normality in the data to permit the use of Gaussian-based methodology (as in Section 3.2 for the fish model). For example, [START_REF] Henderson | Bayesian emulation and calibration of a stochastic computer model of mitochondrial DNA deletions in substantia nigra neurons[END_REF] uses the logit transformation (log y/(1 -y)) in analyzing the proportion of deletions in mitochondrial DNA. [START_REF] Plumlee | Building accurate emulators for stochastic simulations via quantile kriging[END_REF] take a different route by focusing on the quantiles of the output distribution -normality is not needed. Both of these approaches have the appeal of leading to relatively simple modifications of the methods in Sections 3.1 and 3.2.

There are also more complex methods that generally lack the same ease of implementation. For example, [START_REF] Moutoussamy | Emulators for stochastic simulation codes[END_REF] attempt to model the underlying probability density function itself, rather than the output y. Xie and Chen (2017) devise a Student t-process that is not much different than the GP process while at the same time allowing heavier tails in the distribution of the data.4

Quantile Kriging

Quantile Kriging (QK) is an increasingly popular tool for the emulation of stochastic computer models [START_REF] Rannou | Kriging the quantile: application to a simple transmission line model[END_REF][START_REF] Plumlee | Building accurate emulators for stochastic simulations via quantile kriging[END_REF][START_REF] Zhang | Asymmetric kriging emulator for stochastic simulation[END_REF][START_REF] Fadikar | Calibrating a stochastic, agent-based model using quantile-based emulation[END_REF]. These approaches are a natural extension of spatial kriging formulations [START_REF] Zhang | Loss function approaches to predict a spatial quantile and its exceedance region[END_REF][START_REF] Zhou | Estimating the health impact of climate change with calibrated climate model output[END_REF][START_REF] Opitz | INLA goes extreme: Bayesian tail regression for the estimation of high spatio-temporal quantiles[END_REF] used in environmental applications, often with modeling further tailored to account for rare events and extreme quantiles.

The QK method directly models specific quantiles of interest, such as the median and the lower/upper 95% quantiles at each input. Minimal assumptions about the distribution of the simulator output are required. Q q (x), the q th quantile of the simulator output at input x, is modeled with a GP. Given values Q q (x i ) at inputs x 1 , . . . , x n , the quantile, Q q (x new ) for x new , can be predicted using equations 3.1 and 3.2. This framework allows the distribution of the variability v to take on almost any shape. Although a true generative process for the output y is lost, we can describe its distribution.

To implement QK, values of the targeted quantiles at the inputs are needed. Just as in Section 3.2, where sample variance estimates at the inputs can be used, sample quantiles can be used here. Said sample quantiles are calculable given enough repli-cates r i at each x i . The GPs used to predict new quantile values, Q q (x new ), should also include a noise term σ 2 q to acknowledge that the sample quantiles are estimates. Assuming the variability of the sample quantiles is normally distributed may also be invalid, but is at a level further removed from the quantity of interest, y, and is often acceptable in practice.

Including the quantile q as an additional input to the GP model can be a useful modification. The quantile Q q (x) can be reformulated as Q(x, q), increasing the dimensionality of the inputs from d to d + 1. This strategy allows for the prediction of Q(x, q) for any desired quantile q, not just those that were empirically estimated, and is used by [START_REF] Fadikar | Calibrating a stochastic, agent-based model using quantile-based emulation[END_REF] for the Ebola model.

Alternative QK-based approaches are also under development. For example; a promising variant of QK called Asymmetric Kriging (AK, [START_REF] Zhang | Asymmetric kriging emulator for stochastic simulation[END_REF] does not require sample quantiles by leveraging quantile regression methods [START_REF] Koenker | Regression quantiles[END_REF].

Fish Example. For the fish simulator, QK is implemented with the same simulated dataset as before since many replicates are available. The sample 5%, 27.5%, 50%, 72.5% and 95% quantiles at each of the 20 population sizes form the observed data, and the modification using the quantile q as an added input dimension is adopted. Figure 1 presents the predicted Q(x, q) mean for 5 different quantiles along with the data (the left plot) and compares the "true" values with predictions at the 5%, 50%, and 95% quantiles (the right plot). 1 but with a QK surrogate. Mean predictions of 5 quantiles (5%, 25%, 50%, 75% and 95%) are provided on the left along with data, and mean predictions of 3 quantiles (5%, 50%, and 95%) are provided on the right along with the "true" values.

The center purple curve in Figure 3 is the predicted median. The outer red lines are the predicted 5% and 95% quantiles; the inner blue curves are the predicted 25% and 75% quantiles. The non-monotone "wavy" lines for the 5% and 95% quantiles reflect the natural variability of extreme quantiles based on only 20 observations. Without an abundance of replicates, accurately capturing extreme quantiles is difficult, a drawback of QK. The other quantiles presented display more regularity.

The results from QK do not differ much from those in Figure 1 where the squareroot transformation was sufficient. With more complex problems, such as the Ebola model, the method is suitable while other approaches may be less so. In any case, QK can be a good robust choice given adequate data for estimating quantiles.

Multiple Outputs

The discussion thus far has assumed that the simulator outputs a single scalar quantity of interest. For multivariate output a more comprehensive model would be ideally used. Sophisticated approaches building multivariate GPs have been employed [START_REF] Conti | Bayesian emulation of complex multi-output and dynamic computer models[END_REF][START_REF] Fricker | Multivariate Gaussian process emulators with nonseparable covariance structures[END_REF][START_REF] Paulo | Calibration of computer models with multivariate output[END_REF]. Bespoke, problem-specific formulations for time-series and other outputs have also been entertained [START_REF] Farah | Bayesian emulation and calibration of a dynamic epidemic model for A/H1N1 influenza[END_REF][START_REF] Sun | Synthesizing simulation and field data of solar irradiance[END_REF]. If there are a small number of outputs, treating each independently, with its own surrogate model, often suffices. This method can be effective, despite ignoring any correlation between the different outputs and thus wasting information. For example, [START_REF] Spiller | Automating emulator construction for geophysical hazard maps[END_REF] deploy independent surrogates at each of a multitude of sites in a region to good effect.

Alternatively, by treating the index, t, of the T outputs as an additional input dimension (changing the dimension of the input space from d to d + 1) a GP surrogate on d + 1 dimensions can be formed [START_REF] Bayarri | Predicting vehicle crashworthiness: Validation of computer models for functional and hierarchical data[END_REF]. This method allows correlation structures between the different outputs to be modeled. This is similar to the QK modification where quantile levels are treated as an added input (Section 3.3.1). A drawback of this technique is that, if T is very large, computational issues will arise because the GP must be trained on N T data points rather than just N . Intrinsic variability prevents simplifications of the sort used in [START_REF] Bernardo | Integrated circuit design optimization using a sequential strategy[END_REF] for deterministic simulators in this setting.

A different approach reduces the effect of the size of T to a smaller K 0 by representing the output through the use of basis functions, ψ(t):

y(x, t) = K 0 k=1 w k (x)ψ k (t) + δ(x, t). (3.5) Coefficients w k (x i ) k = 1, .
. . , K 0 are determined by the data; and δ(x, t) is the residual error between the basis function representation and the data y.

If K 0 = T then δ = 0.
Typically, K 0 is taken to be much less than T but large enough so that the error, δ, is sufficiently small. Each w k (x) can be independently modeled with a surrogate and predictions for y(x, t) are obtained from equation 3.5, ignoring δ.

Different choices for the bases can be appropriate in different settings. For example, Bayarri et al. (2007a) use wavelets for the ψ k s in a deterministic setting where t is time. A common choice of basis functions are principal components: the ψs are the eigenvectors of the matrix Y N Y N , the first K 0 of which are in correspondence with the first K 0 eigenvalues in decreasing order. It is often the case that the first few (five or less) principal components are enough to capture sufficient information about the full (T ) data set. Coefficients w k (x i ) are then equal to T t=1 y(x i , t)ψ k (t). More information about principal components can be found in [START_REF] Jolliffe | Principal Component Analysis[END_REF] and software for obtaining ψ and w k is prevalent.

Further discussion about using principal components to model high-dimensional simulator output can be found in [START_REF] Higdon | Computer model calibration using high-dimensional output[END_REF]. Principal components are also utilized in [START_REF] Fadikar | Calibrating a stochastic, agent-based model using quantile-based emulation[END_REF] to model the time-series output of the stochastic Ebola simulator. While principal components are a common default, there is concern that key features of the data set may be left within the discarded δ(x, t) preventing reliable prediction. [START_REF] Salter | Uncertainty quantification for computer models with spatial output using calibration-optimal bases[END_REF] document these concerns with regards to calibration and suggest an alternative.

For problems with functional outputs, with potentially missing data and/or irregularly spaced data (such as irregularly spaced timesteps or spatial locations), a functional decomposition can also be useful. For example, [START_REF] Ma | Computer model emulation with high-dimensional functional output in large-scale observing system uncertainty experiments[END_REF] use functional principal component analysis to model satellite observation simulations.

Experimental Design

For an experiment, the design (the choice of x values) and analysis (the assessment of the output y(x)) are, in principle, closely connected. Other considerations can also enter. For physical experiments, controlling for external influences or nuisance factors by blocking and randomization is often a vital part of the design. External influences are absent in computer experiments and so controlling for nuisance factors is usually irrelevant. However, many minor parameters are often fixed which could instead be randomized over, with a consequent addition to intrinsic error.

With a specific goal (e.g., predicting simulator output) and a criterion of accuracy (e.g., the average prediction uncertainty: the integrated mean-squared prediction error, IMSPE5 ), designs that optimize the criterion are preferred. Since the criterion will usually depend on the surrogate, which, in turn, depends on unknown parameters, what to use as a stand-in for the parameters before any data are collected is an issue. Extensive study of single-stage deterministic computer experiments resolved this dilemma by downplaying optimality and recommending readily computed "spacefilling" designs where no large region of input space is missed. Space-filling designs are readily computed, whereas optimizing IMSPE is complicated and without substantial advantage. For practical adoption, designs must be easy to produce as well as effective.

Multiple methods exist for obtaining space-filling designs, the most popular being Latin hypercube designs (LHDs;[START_REF] Mckay | Comparison of three methods for selecting values of input variables in the analysis of output from a computer code[END_REF]. 6 LHDs have proved adequate, especially when joined with an additional criterion, such as the maximin criterion, where one also maximizes the minimum distance between points in the design.7 Even a random LHD will often suffice. Sobol sequence designs [START_REF] Sobol | On the distribution of points in a cube and the approximate evaluation of integrals[END_REF] are equally effective for predicting the output of a deterministic simulator. The xs for Sobol designs are generated sequentially making it easy to retain the space-filling character when a multi-stage or sequential design strategy is used. 8 Pronzato and Müller (2011) have a lengthy discussion of these and other space-filling methods, some pertinent to non-rectangular geometries.

For stochastic simulators the picture is far less clear. The presence of intrinsic variability raises the complication of replication, not present in deterministic experiments. With the same inputs, a stochastic simulator can be run multiple times (replicated) providing different output values each time due to the intrinsic randomness. Replicates obviously have an effect on the estimation of the intrinsic variance, σ 2 v , and therefore on prediction (see Section 3.2), and so the number and location of replicates are important. A simple approach for a single-stage experiment is to use a space-filling design to establish the sites X n = (x 1 , . . . , x n ) of the experiment and then add replicates at each site. Determining the number, r i , of replicates at each site x i and how to apportion between replicates and sites, that is, how to choose the number of unique sites, n, given a total simulation budget N , is not well understood. In fact, there is limited theoretical evidence of the need for replicates altogether, although there is numerical evidence and wide belief that replicates can be advantageous, at least in appropriate contexts. For example, [START_REF] Wang | Controlling sources of inaccuracy in stochastic kriging[END_REF] produce designs by minimizing bounds on IMSPE. Their numerical results show no need for replicates unless σ 2 v (x) is large compared to σ 2 Z (a factor in measuring uncertainty in estimating the mean M ).

The presence of intrinsic variability suggests there is value in multi-stage designs where stage 1 is used to get information about σ 2 v (x) and later stages exploit this information to allocate replicates and select new inputs. Questions arise as to how inputs should be selected for stage 1, and also how to leverage the results from stage 1 to select new inputs and replicates in later stages. The two factors, replication and multiple stages (including fully sequential), are central to developing adequate design strategies. Attention is paid to both factors in the discussion below.

Single-Stage Design

A common approach in single-stage studies is to use space-filling designs for inputs, say n in number, and r replicates at each input, sometimes with no repeats i.e., r = 1. Predictions follow as described in Section 3 depending on the particular prediction model selected. Choices have to be made about the total number of runs and the number of replicates at each input site (N = nr). Often, N is a question of budget, but there is little insight into how r should be chosen except when meeting a specific surrogate model requirement, as in SK (Section 3.2).

For their single-stage study, [START_REF] Marrel | Global sensitivity analysis of stochastic computer models with joint metamodels[END_REF] use a standard LHD with no repeats to compare the performance of different statistical models. On the other hand, Plumlee and Tuo (2014) use a LHD with varying numbers of replicates r i at each x i . In their case, the number of replicates must be large, because the QK method (Section 3.3.1) depends on computing quantiles of the output y(x i ) at each input site of the design.

Two-Stage Design

The case for a two-stage design is largely to enable estimation of σ 2 v at stage 1 and use it for the second stage. [START_REF] Ankenman | Stochastic kriging for simulation metamodeling[END_REF] provide one solution in the context of SK. A first-stage design chooses the x i s via an LHD of size n 1 with a common number, r, of replicates at each of the inputs, resulting in a total number of N 1 = n 1 r runs at stage 1. The first-stage analysis uses the r replicates at each input to estimate σ 2 v (x i ) using the sample variances. As outlined in Section 3.2, a GP (working with log s 2 (x i )) is then used to produce a "plug-in" estimate of σ 2 v (x) for all x. A different GP uses that variance estimate to build a predictor for the mean output M .

For stage 2, n 2 additional unique input locations are chosen so that the combined set of design locations, X n = (x 1 , . . . , x n ), remains space-filling. The IMSPE is then calculated by integrating the MSPE all possible inputs X, using the GP model constructed in stage 1. Minimizing the IMSPE with respect to the number of replicates R n = (r 1 , . . . , r n ) provides the optimal number of replicates for the chosen X n . Details are in [START_REF] Ankenman | Stochastic kriging for simulation metamodeling[END_REF]. One difficulty is that the optimal R n might produce an r i for a first-stage site that is smaller than the r already used at stage 1. Some fix to the method would then be necessary.

In this setting, a Sobol sequence could be used to obtain a design that is space filling at both stage 1 and stage 2. This is not what is done in [START_REF] Ankenman | Stochastic kriging for simulation metamodeling[END_REF], but a Sobol sequence is easier to implement and likely to yield similar results. Choosing the unique inputs X n for stage 2 by optimizing the IMSPE could also be done, but adds to the computational burden. Suitable recommendations for the values of n 1 , n 2 , N and the replicates at each distinct input are lacking (in [START_REF] Ankenman | Stochastic kriging for simulation metamodeling[END_REF] the recommendations are ad hoc) and, as for one-stage experiments, open for study. A third-stage design (or indeed, any multi-stage design) can be constructed by repeating stage 2 in the above process.

Sequential Design

When the statistical design and resulting analysis are closely coordinated, it may be feasible to carry out a sequential process whereby, after the first stage, a run is chosen one-at-a-time. After each run all quantities of relevance can be updated in order to determine the next run. This addresses the issue of learning about σ 2 v and obtaining new runs without pre-specifying their allocations. An advantage of a sequential design is the possibility of stopping when a criterion is met before a budget constraint is reached. Another advantage is the increased likelihood of making useful runs of the simulator, replicates or otherwise. For some objectives, such as optimization (Section 6.3), a sequential design is usually essential. For global prediction, [START_REF] Binois | Replication or exploration? Sequential design for stochastic simulation experiments[END_REF] present an approach to sequential design, implemented in the previously mentioned hetGP package.

The strategy in [START_REF] Binois | Replication or exploration? Sequential design for stochastic simulation experiments[END_REF] begins at stage 1 with a space-filling design D 1 of n 1 inputs and an allocation of runs (r(x 1 ), . . . , r(x n 1 )). Using a GP for M and a latent GP prior on σ 2 v , as in Section 3.2, a MLE computation deals with all parameters, leads to predictors, and a calculable estimate of IMSPE(D 1 ). A new point z is considered, either as a new unique input x n 1 +1 or as a replicate of an existing input in D 1 . Selection z is added to the design D 1 if z minimises IMSPE(D 1 + z), yielding a new design D 2 . This myopic rule can be iterated and each time a new point is added the surrogate, including MLEs of its parameters is updated. The process stops when a criterion is met or the computational budget exhausted.

Computational viability is strained by the updating required after each run. On the other hand, the computational burden is eased by nature of it being "greedy": it only seeks the optimal data point for the very next simulator run, ignoring runs that may be better in the long run. This is not the only sequential design scheme available for global prediction problems. For example, the tree-generating processes used in TGP and BART (see Section 6.1) deliver specialized sequential design strategies. Details are available in [START_REF] Gramacy | Adaptive design and analysis of supercomputer experiments[END_REF] and [START_REF] Chipman | BART: Bayesian additive regression trees[END_REF].

Blurring the lines between multi-stage and sequential designs, it can sometimes be practical to run additional simulations in batches (e.g., as in making efficient use of a multi-core supercomputer). In such circumstances a "batch design" would be desirable. These have been developed for deterministic experiments (Loeppky et al., 2009a;[START_REF] Duan | Sliced full factorial-based latin hypercube designs as a framework for a batch sequential design algorithm[END_REF][START_REF] Erickson | Gradient based criteria for sequential design[END_REF], but not yet explicitly extended to stochastic cases.

When fully sequential methods are feasible the seqhetGP strategy sketched above is valuable. There are several aspects worth examining:

• The extensive use of a surrogate in the construction of the design requires scrutiny by diagnostics that assess the quality of the surrogate.

• The first stage of a sequential strategy must avoid a poor (e.g., too small) initial design lest a poor starting surrogate leads to poor choices thereafter.

• The utility of a sequential design depends on the relative cost of implementation compared to simulator runs. For challenging problems simulator runs are likely to be costly enough to make sequential design attractive.

• There may be modifications to a sequential design that reduce computational load without paying a significant cost in accuracy. For example, re-estimate parameters periodically rather than after each step.

Ocean Example. For the ocean model, we use an initial design of 50 sites, chosen by a maximin LHD in 2-d, each site with 5 replicates. The remaining 750 data points are then assigned via the sequential scheme. The resulting mean and standard deviation surfaces are in Figure 4. For the standard deviation surface the design sites are superimposed along with the number of replicates taken at the sites. The mean surface in the left panel is slightly different than for the non-sequential analyses (Figure 2, top row). The standard deviation surfaces look very different. For the design itself, new inputs are heavily replicated in regions where the standard deviation is large, and less so in regions where it is small. Additionally, the sequential design includes more unique sites than the fixed design, and more points on the boundaries of the input space.

Using the "truth" established in Section 3.2 we can compare the performances of the three methods. As discussed previously, the visual presence of heteroscedasticity is a deterrent to using homGP. Visually distinguishing between the performances of the hetGP and seqhetGP surrogates is more difficult: the means appear similar, and whilst some patterns in the true standard deviation appear to be captured by hetGP, imperfections are visible and the magnitude is not always correct. With the seqhetGP standard deviation, nuance seems lost. To properly compare the different methods, a numerical comparison can be more valuable.

Two useful numerical measures are root mean squared error, RMSE (the squareroot of the average squared difference between the surrogate's prediction of the mean and the "true" mean) and Score (the proper scoring rule from equation 27 in [START_REF] Gneiting | Strictly proper scoring rules, prediction, and estimation[END_REF]). RMSE measures the accuracy of the mean predictions and Score is an overall measure testing the accuracy of the combined mean and variance predictions. With a test set of inputs x 1 , . . . , x p and simulator outputs y 1 , . . . , y p , surrogate predictive means µ N (x 1 ), . . . , µ N (x p ) and variances σ 2 N (x 1 ), . . . , σ 2 N (x p ), Score is

1 p p i=1   - y i -µ N (x i ) σ 2 N (x i ) 2 -log(σ 2 N (x i ))   . (4.1)
Smaller RMSE is better while for Score, larger is better.

For the three methods, the RMSE for homGP, hetGP and seqhetGP are respectively 2.056, 1.985, and 1.567; and the Scores are respectively -3.999, -3.880, and -3.834. The RMSE results reveals that seqhetGP is best at predicting the mean, which was not obvious from the plots. The Scores for hetGP and seqhetGP are close but noticeably better than homGP, affirming the presence of heteroscedasticity.

The randomness in stochastic simulators as well as variability in design (there are many possible maximin LHDs) can induce a large degree of variability in specific results such as those just cited. It is therefore difficult to rely on a single result for making comparisons. As such, the above experiment is repeated 100 times and the resulting 100 RMSEs and Scores are summarized in boxplots in Figure 5. The boxplots confirm what was found with the single data set: heteroscedasticity is present and seqhetGP is preferred. Visual inspection of many of the standard devi-ation plots for the repeated experiments (as discussed in Section 3.2 and found in the supplement) reveals considerable variation and departure from the true standard deviation. The variance (the intrinsic variance and the GP uncertainty for the mean) can be hard to get right without an abundance of data, and the difficulty is compounded by the use of plug-in estimates whose uncertainty is not accounted for.

Designing for Statistical Model Parameter Estimation

Sections 4.2 and 4.3 construct designs that rely on a surrogate model based on stage 1 data in order to choose subsequent data points. The quality of the designs depends on the accuracy of the surrogate which, in turn, depends on the accuracy of its parameters.

An alternate approach to those used in Sections 4.2 and 4.3, is to construct an initial design with the express purpose of better estimating these parameters. Boukouvalas et al. (2014a) address the problem and focus on hetGP models, using a simple parametric function for the variance (σ 2 v (x) = exp(h(x)), where h is a simple function (e.g., a polynomial). They propose designs that maximize a criterion previously used for deterministic simulators by [START_REF] Abt | Fisher information and maximum-likelihood estimation of covariance parameters in Gaussian stochastic processes[END_REF]: the logarithm of the determinant of the Fisher information matrix, log |I|. Numerical results suggest this method gives improvements in estimating the parameters, but overall global prediction is no better, and sometimes worse, than using a space-filling design. When prediction is of prime importance, the question arises about how to make use of such designs for stage 1 in a multi-stage or sequential setting, where its impact on obtaining better initial surrogate models can be felt. For example, see [START_REF] Zhang | Distance-distributed design for Gaussian process surrogates[END_REF].

Calibration

Calibration is needed when there are inputs to the simulator that are neither known nor measurable, a common condition in practice. Transmissibility in the Ebola simulator and the diffusion coefficients in the ocean model are examples of such inputs. In order to infer (indirectly) values for these inputs and produce predictions, added information in the form of field data (experimental or otherwise) are necessary. Inclusion of field data and calibration parameters, labelled u C , leads to the observation model:

y F (x) = y S (x, u C ) + δ MD (x) + , (5.1)
where y F (x) are real-world field observations at controllable (or measurable) inputs x, y S is the simulator with additional unknown, non-measurable, inputs u C , is measurement error for the observations y F (x) (with variance σ 2 ), and δ MD (x) is an important term that accounts for the simulator not being a perfect representation of reality. y F "observes" reality with error ; reality = y S + δ MD .

Multiple competing methodologies and even philosophies exist for calibration. Several solutions to the calibration problem are outlined below. Despite the centrality of calibration in computer experiments, comprehensive comparisons are lacking.

Kennedy-O'Hagan Calibration (KOH)

The formulation in equation 5.1 was made by [START_REF] Kennedy | Bayesian calibration of computer models[END_REF] for deterministic simulators and is the basis for much of the calibration and related prediction work since. The strategy pursued by [START_REF] Kennedy | Bayesian calibration of computer models[END_REF], as implemented in [START_REF] Bayarri | A framework for validation of computer models[END_REF], obtains a surrogate for y S and models δ MD (x) with a GP (although other choices are possible). After replacing y S with the surrogate, posterior distributions for all unknowns can be obtained via a Bayesian analysis. In practice, the surrogate model is fit only using the simulator data, ignoring possible influences from the field data. Details and discussion of this modular approach can be found in [START_REF] Bayarri | A framework for validation of computer models[END_REF] and [START_REF] Liu | Modularization in Bayesian analysis, with emphasis on analysis of computer models[END_REF].

The KOH approach emphasizes the necessity to address calibration and model discrepancy together. Confounding between u C and δ MD (x) inevitably occurs because there are multiple combinations of u C and δ MD (x) that result in the same observed field data. Thus, u C is non-identifiable and its estimation is compromised, as is the discrepancy. Nonetheless, the resulting predictions for y and E(y) are sound, even if the individual estimates for u C and δ MD (x) aren't. For details and further discussion see [START_REF] Higdon | Combining field data and computer simulations for calibration and prediction[END_REF], [START_REF] Bayarri | A framework for validation of computer models[END_REF], [START_REF] Brynjarsdóttir | Learning about physical parameters: The importance of model discrepancy[END_REF], and [START_REF] Tuo | A theoretical framework for calibration in computer models: parametrization, estimation and convergence properties[END_REF].

Multiple attempts to circumvent confounding have surfaced. [START_REF] Tuo | Efficient calibration for imperfect computer models[END_REF] alleviates the ambiguity in u C by formally defining it as a least-squares quantity; [START_REF] Gu | Scaled Gaussian stochastic process for computer model calibration and prediction[END_REF] propose novel priors for the discrepancy that compromise between the [START_REF] Tuo | Efficient calibration for imperfect computer models[END_REF] strategy and KOH; and Plumlee (2017) introduces priors on the discrepancy that are orthogonal to the prior mean. In the stochastic simulator literature, [START_REF] Oakley | Calibration of stochastic computer simulators using likelihood emulation[END_REF] removes δ MD but compensates by inflating the variability in the prior distribution for u C . Ignoring δ MD altogether can be justified by strong evidence of the simulator being accurate, but such evidence is rare.

For stochastic problems, where reality is stochastic the discrepancy term δ MD (x) cannot be assumed deterministic. Modeling the discrepancy is likely be influenced by the model for the simulator while recognizing that discrepancy is often smoother. For example, if modeling y S calls for a hetGP with the Matern 5/2 correlation function then it is likely that a hetGP is needed for the discrepancy, perhaps with the smoother squared exponential correlation. A full Bayesian analysis in such circumstances may be prohibitively expensive and the above procedure would have to be modified. [START_REF] Sung | Calibration of computer models with heteroscedastic errors and application to plant relative growth rates[END_REF] use a hetGP for the discrepancy (but for a deterministic simulator), estimating parameters via maximum likelihood and following [START_REF] Tuo | Efficient calibration for imperfect computer models[END_REF] to avoid confounding.

Revisiting Ebola. The Ebola study [START_REF] Fadikar | Calibrating a stochastic, agent-based model using quantile-based emulation[END_REF] calibrates an ABM using the KOH framework. The simulator y S has 5 unknown, unmeasured inputs u C and the output is the log of the cumulative number of infected individuals up to week 1 and every week thereafter up to 57 weeks. The field data y F is a set of reported cumulative counts. For the statistical model, a QK strategy (Section 3.3.1) is followed by replicating each distinct simulation 100 times and then condensed into evenlyspaced quantiles at each time point (specifically, the 5% 27.5%, 50%, 72.5% and 95% quantiles). These quantile output trajectories are then reduced to a more manageable 5 dimensions using the principal component decomposition outlined in Section 3.4.

Underlying the approach is an assumption that the epidemic trajectories (actual and simulated) can be approximated by quantile trajectories (i.e., a realized epidemic that resembles the q th quantile at time 1 will also resemble the q th quantile at a later time). Accordingly, the quantile q is included as an input parameter (see Section 3.3.1) to allow KOH calibration to learn about the 5 calibration parameters as well as the (unknown) value of q for the observed epidemic. Because the difference between the simulator and reality quantile trajectories could not be noisy, the discrepancy is treated as deterministic (a smoothing spline is used rather than a GP for the discrepancy). Posterior distributions for unknown u C , δ MD , and q are obtained and used to make predictions of the cumulative counts and other quantities.

In the main analysis, which restricts the field data to only the first 20 weeks, the estimate of model discrepancy is almost zero. A subsequent analysis done using field data up to week 42 exposes some inaccuracy of the simulator (non-zero δ MD (x)) -the simulator continues to predict infections, even after the epidemic has died in reality.

Ocean Example. The previous ocean analyses fixed the two diffusion coefficients. Realistically, they are unknown and calibration is necessary. "Field" data are artificially created by averaging over 200 simulations at 150 different longitude-latitude coordinates, using the previously fixed values of the diffusion coefficients (K x = 700 and K y = 200). "True" values are obtained by adding a fake discrepancy, taken as a single realization from a GP with a squared-exponential correlation function, a variance of 1.64, and θ values of (1, 2) (equation 3.3). To these, normally distributed pretend "observation errors" with a variance of 4 are added, two such observations at each site. In real problems, the field data would be observed and not generated like this. Note that field data for this problem corresponds with the mean of the simulator, not individual draws from the simulator; a result of the simulator being a stochastic approximation.

With the diffusion coefficients now uncertain, the simulator has four inputs. A computer experiment is designed with runs at the 150 sites used for the field data and 500 unique selections of the calibration parameters K x and K y . This is done by combining copies of the 150 longitude and latitude sites with a size-500 maximin LHD for (K x , K y ), and then improving the combined design by maximizing the minimum distance between design points in the 4-dimensional space. Call this set of points D oc . The simulator experiment is carried out by taking 10 replicates at each point in D oc . Two distinct surrogates (a homGP and a hetGP) are fit with this fixed design. In addition, a seqhetGP surrogate is constructed, with an initial design of only 4 replicates of D oc and the remainder of the budget assigned following the strategy of [START_REF] Binois | Replication or exploration? Sequential design for stochastic simulation experiments[END_REF], described in Section 4.3.

For a KOH analysis done in modular fashion the surrogates are fit only using the simulated data. Because reality here is represented by the expectation of the simulator (rather than the simulator output itself), y S in equation 5.1 is replaced with E(y S ). Similarly, because reality is deterministic, δ MD is modeled as a standard GP. Of course, the simulated data are outputs from y S , not from E(y S ) -the surrogate is used to approximate the deterministic E(y S ). MCMC is then used to obtain posterior distributions for the remaining unknowns: the diffusion coefficients, K x and K y ; the variance and correlation parameters of the model discrepancy GP, σ 2 MD and θ MD ; and the observational error, σ 2 ). The posterior distributions for the key parameters are in Figure 6; their true values are K x = 700, K y = 200, σ 2 MD = 1.64 and σ 2 = 4. For all three surrogate models the posterior distributions for K x are fairly diffuse. The K y posteriors are highly concentrated, but not quite around the true value. The three posteriors for observational error are quite similar but all point to estimates closer to 5 rather than the true 4. The posteriors for the discrepancy variance are diffuse. These plots underline the dilemma of calibration: obtaining accurate values of calibration (and other) parameters in the presence of model discrepancy is problematic. Additionally, with noisy data, it is difficult to obtain precise estimates. However, KOH does yield useful posterior predictive distributions.

Table 1 compares predictions by KOH calibration with 3 other calibration approaches. The first one estimates K x and K y by ordinary least squares (OLS): (K x , K y ) is chosen such that the sum of the squared residual difference between the mean surrogate prediction and the observed data is minimized. New observations are then predicted by running the surrogate with the parameters (K x , K y ) replaced by the OLS estimates ( Kx , Ky ). The second approach follows a frequently adopted practice by guessing, or "judiciously selecting", specific values for K x and K y . Here, the choices K x = 600 and K y = 400 are made, and then predictions are made using the surrogate. Call this method SINGLE. The third method, NOCAL, generates predictions as if there were no field data and the distribution for (K x , K y ) is taken as their prior distribution, independent uniform priors on [100,1000]. In these alternative methods, the observational error variance is fixed at the true value, and 0-discrepancy is assumed (the former is overly generous and the latter is all too common in practice). For NO-CAL, a distribution for the oxygen concentration is obtained by sampling values of K x and K y from their prior distribution and plugging them into the surrogate, while for KOH, by sampling from the posterior distributions of all unknowns. 5; similarly for Score. ( Kx , Ky ) are least squares estimates for (K x , K y ), OLS presents the predictive results from least squares calibration, SINGLE the results from arbitrarily choosing (600, 400) for the diffusion coefficients, NOCAL the results from sampling the prior for (K x , K y ), and KOH the results from performing KOH calibration.

Although the differences in RMSE are negligible, the Scores indicate that KOH performs the best. It is also possible that the accuracy of OLS, SINGLE, and NOCAL is overstated, because the observational error variance is taken as known while in KOH it is estimated. That the least squares estimates ( Kx , Ky ) are not always close to the true values is unsurprising given the presence of discrepancy, along with possible imperfections and high variability in the surrogate. For similar reasons, scant differences appear among the three surrogates.

The similarity of RMSEs is a consequence of large variability in the surrogate, the presence of discrepancy, the dominance of the longitude and latitude inputs, and a weak effect from the calibration inputs. The first explains the magnitude of the RMSEs and the last explains why fairly inaccurate calibration inputs (in OLS and SINGLE) don't matter. Because KOH addresses discrepancy, its Score exceeds the others', showing that accounting for the discrepancy is necessary and can't be wished away Estimating Calibration Parameters The focus here, consistent with Section 4, has been on improving global prediction. If the problem was instead to provide good estimates for calibration parameters when model bias is absent, then different designs may be better suited. [START_REF] Damblin | Adaptive numerical designs for the calibration of computer codes[END_REF] address this in the context of deterministic simulators, but it is unclear how such methods extend to stochastic simulators. Additionally, whilst KOH facilitates capable predictions, the complexity and recorded pitfalls in KOH has led to competing calibration techniques that are also in common usage.

History Matching (HM)

History Matching (HM) is a common alternative to KOH calibration [START_REF] Craig | Pressure matching for hydrocarbon reservoirs: a case study in the use of bayes linear strategies for large computer experiments[END_REF][START_REF] Vernon | Galaxy formation: a Bayesian uncertainty analysis[END_REF][START_REF] Boukouvalas | Bayesian precalibration of a large stochastic microsimulation model[END_REF][START_REF] Andrianakis | Efficient history matching of a high dimensional individual-based HIV transmission model[END_REF]. HM searches for inputs where the simulator outputs closely match observed data, while recognizing the presence of the various uncertainties, including model discrepancy. The HM approach rules-out "implausible" inputs in a straightforward way, rather than attempting to find probable inputs. With an observation y F , and initially assuming u C makes up all the unspecified simulator inputs, u C is deemed implausible if:

|y F -µ N (u C )| σ 2 N (u C ) + σ 2 MD + σ 2 ) ≥ 3, (5.2)
where σ 2 N , σ 2 MD , and σ 2 are the variances of the surrogate, the model discrepancy, and the observational error respectively. In other words, an input is implausible if the difference between the observation and the simulator output using that input is sufficiently large relative to those uncertainties. The number 3 comes from [START_REF] Pukelsheim | The three sigma rule[END_REF] who shows that at least 95% of any unimodal distribution is contained within three standard deviations. When there are multiple outputs or additional, controllable inputs there are modifications to equation 5.2 [START_REF] Vernon | Galaxy formation: a Bayesian uncertainty analysis[END_REF].

The process can be repeated in so-called "waves", using non-implausible u C found at one wave to generate simulation runs for the next wave, sequentially reducing the space where u C could lie. With these waves HM aims to avoid regions of inputs where u C is unlikely to be and, in that regard, HM is a calibration design scheme. At any given wave, it is possible for all values of u C to be deemed implausible -the so-called terminal case [START_REF] Salter | Uncertainty quantification for computer models with spatial output using calibration-optimal bases[END_REF] -usually implying that σ 2 MD is set too low or that the simulator is not fit for purpose. [START_REF] Andrianakis | Bayesian history matching of complex infectious disease models using emulation: a tutorial and a case study on HIV in Uganda[END_REF] contains a thorough description of HM whilst applying it to a complex epidemiology model of HIV.

HM and KOH.

With KOH the estimation of u C is confounded with discrepancy, but predictions and their uncertainties are available. However, implementing KOH in complex problems may be burdensome if not intractable. Speculatively, a hybrid strategy may be to use HM to reduce the input space, confirm the absence of the terminal case, and then apply KOH in the narrowed space to get predictions and uncertainties. Complex models, unlike the fish and ocean examples in this review, would be ones for which this approach would be most appealing. Such hybrid strategies are a topic for further exploration.

Approximate Bayesian Computation (ABC)

Obtaining posterior distributions for calibration parameters u C and predictions (e.g., the KOH approach in Section 5.1) can be computationally challenging. ABC methods offer an alternative which have been found useful in moderately complex contexts [START_REF] Rutter | Microsimulation model calibration using incremental mixture approximate Bayesian computation[END_REF]; but less so in more ambitious settings [START_REF] Mckinley | Approximate Bayesian computation and simulation-based inference for complex stochastic epidemic models[END_REF].

ABC aims to produce samples from π(θ|Y F ), the posterior distribution of unknowns θ, given the field data Y F . For calibration, think of θ as u C . ABC does this by generating samples for the unknowns θ (s) and the output z (s) from π(Y F |θ)π(θ), that is, from the likelihood of the data given the unknowns, multiplied by the prior probability of the unknowns. For computer models, generating samples from the likelihood is equivalent to running the simulator. Such samples are only accepted if z (s) = Y F . For continuous settings, where exact equality cannot occur, acceptance is instead made if B(z (s) , Y F ) < τ , where B is a measure of distance and τ a level of tolerance. An approximated posterior distribution is then given by the collection of accepted θ (s) s. When there are multiple outputs (or there are other controllable inputs x, and so for any given θ (s) there are effectively multiple outputs), Y F and z (s) can be replaced with informative summary statistics. Finding a single statistic sufficient for all outputs is challenging, and a poorly chosen one can invalidate results.

The choice of the tolerance τ is important. If τ is small then it may take a very long time to generate a single sample that satisfies the inequality. If τ is not small then the approximation to the posterior is less reliable. For calibration, τ can be interpreted as a bound on the observational error and model discrepancy, leading to a "correct" posterior rather than an approximation [START_REF] Wilkinson | Approximate Bayesian computation (abc) gives exact results under the assumption of model error[END_REF]. This is then similar to HM with the subjective choice of bounds.

ABC can be done without the use of a surrogate, but many runs of the simulator itself to generate many θ (s) may be required. Otherwise, too few accepted θ will remain, or an overly high value of τ will be required. In either case accuracy can be compromised. Such computational barriers can be alleviated by the use of a surrogate.

Fish Example. Here we apply ABC to the fish simulator in order to estimate how many fish are in the population. Suppose that 25 fish are recaptured in the second round. A straightforward method to determine the total population size is to simulate many times from the NetLogo fish model, for many different values of the total fish population, and "accept" a simulation every time it leads to 25 fish being recaptured. This is exactly ABC, and is a fairly common practice with ABMs. Doing so 10,000 times, using a uniform prior on the integers between 200 and 4000, so each such population size has prior probability 1/3801, yields the results in the left panel of Figure 7. This direct use of the simulator produces only 52 accepted samples, which is a very small number, and this from 10,000 simulator runs. In comparison, a hetGP surrogate fit from only 400 runs, and from which 1,000,000 samples can be quickly drawn, yields 3811 accepted samples. This result, illustrated in the right panel of Figure 7, gives a less noisy histogram with the same overall shape. If the agent-based model is even marginally costly then a surrogate is unquestionably valuable for ABC computations. q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 1000 2000 3000 4000 0 20 40 60 Population Number of Marked in Recapture q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q qq q q q q q q q q q q q 0 500 1000 1500

Related Calibration Techniques

Bound-to-Bound [START_REF] Frenklach | Comparison of statistical and deterministic frameworks of uncertainty quantification[END_REF] is akin to HM, where an error bound that sweeps up all uncertainties is similarly defined and quadratic programming is then used to find feasible bounds for u C . Bayesian Melding [START_REF] Poole | Inference for deterministic simulation models: the Bayesian melding approach[END_REF][START_REF] Raftery | Inference from a deterministic population dynamics model for bowhead whales[END_REF] is a technique related to Bayesian calibration, used to reconcile differences between elicited prior distributions on inputs and outputs of a simulator. It has been applied in ecology, epidemiology, urban modeling, and pollution monitoring [START_REF] Ševčíková | Assessing uncertainty in urban simulations using Bayesian melding[END_REF][START_REF] Alkema | Probabilistic projections of HIV prevalence using Bayesian melding[END_REF][START_REF] Radtke | Bayesian melding of a forest ecosystem model with correlated inputs[END_REF][START_REF] Fuentes | Model evaluation and spatial interpolation by Bayesian combination of observations with outputs from numerical models[END_REF].

Other Methods and Objectives

Here we briefly outline other surrogate modeling and downstream tasks.

Regression Trees

In some situations the simulator mean M may have discontinuities or "regime changes", where a very different relationship between y and x exists in one part of input space compared to another part (i.e., non-stationarity). Regression Trees [START_REF] Breiman | Classification and Regression Trees[END_REF] form a class of methods that can be useful in these situations. They are also useful in contexts where some inputs are categorical rather than numerical. The problems are treated by dividing the input space into mutually exclusive regions within which independent surrogates (GPs or other regression methods) are fit. Two approaches: the treed GP (TGP [START_REF] Gramacy | Bayesian treed Gaussian process models with an application to computer modeling[END_REF] and Bayesian Additive Regression Trees (BART [START_REF] Chipman | BART: Bayesian additive regression trees[END_REF] have found wide application. Both use the data to automatically partition the input space, rely on Bayesian computation, and have public software: TGP in tgp on CRAN [START_REF] Gramacy | tgp: Bayesian Treed Gaussian Process Models[END_REF][START_REF] Gramacy | tgp: an R package for Bayesian nonstationary, semiparametric nonlinear regression and design by treed Gaussian process models[END_REF]; BART in several R packages, including BayesTree [START_REF] Chipman | BayesTree: Bayesian Additive Regression Trees[END_REF] and BART (McCulloch et al., 2019).

Other approaches by [START_REF] Rullière | Nested kriging predictions for datasets with a large number of observations[END_REF], and via Voronoi tessellations instead of trees (e.g., [START_REF] Kim | Analyzing nonstationary spatial data using piecewise Gaussian processes[END_REF][START_REF] Rushdi | VPS: Voronoi piecewise surrogate models for high-dimensional data fitting[END_REF][START_REF] Park | Patchwork kriging for large-scale Gaussian process regression[END_REF], have received less attention. [START_REF] Pratola | Heteroscedastic BART via multiplicative regression trees[END_REF] extends BART to heteroscedastic σ 2 v (HBART) by modeling M as a sum of Bayesian regression trees (as in BART) and the intrinsic variance σ 2 v (x) as a product of Bayesian regression trees, in a joint approach similar to that in Section 3.2.

Calibration methods capitalizing on the KOH approach and using regression trees as in TGP (Section 6.1) are explored in [START_REF] Konomi | Bayesian treed calibration: an application to carbon capture with AX sorbent[END_REF]. In each terminal node of the partition a GP with an independent constant intrinsic variance term is assumed for the computer model output. An independent GP is also deployed for the discrepancy term. Though σ 2 v is constant at each terminal node the constants can vary across the terminal nodes so heteroscedasticity is automatically incorporated.

Qualitative Inputs

Categorical (qualitative) variables are often present in stochastic simulators, especially those that incorporate characteristics of human behavior. While regression trees are capable of dealing with categorical inputs [START_REF] Broderick | Classification and categorical inputs with treed Gaussian process models[END_REF][START_REF] Gramacy | Categorical inputs, sensitivity analysis, optimization and importance tempering with tgp version 2, an R package for treed Gaussian process models[END_REF], GPs may be more effective as surrogates for smooth simulator output. [START_REF] Qian | Gaussian process models for computer experiments with qualitative and quantitative factors[END_REF], [START_REF] Zhou | A simple approach to emulation for computer models with qualitative and quantitative factors[END_REF][START_REF] Chen | Stochastic kriging with qualitative factors[END_REF] describe ways to extend the kernels used for numerical inputs to incorporate qualitative variables. Painting with a broad brush, their approaches take the correlation between two outputs y(x i ) and y(x j ) as the product of two correlation functions: C c (w i , w j ) dealing with the continuous inputs, w, and C q (z i , z j ) for the qualitative variables, z. A simple way of building C q takes C q (w i , w j ) = K k=1 τ k,w ik ,w jk (6.1)

where K is the number of qualitative variables and τ k,w ik ,w jk represents the correlation between w ik and w jk . One example of τ k,w ik ,w jk is:

τ j,w ik ,w jk = exp{-(φ ik + φ jk )I[w ik = w jk ]} (6.2)
where I is the indicator function (= 1 if its argument is true, = 0 if false), and φ > 0.

The cited references also provide other ways of modeling τ k,w ik ,w jk . Alternative methods exist e.g., [START_REF] Zhang | A latent variable approach to Gaussian process modeling with qualitative and quantitative factors[END_REF] make use of latent variables for qualitative models.

Optimization

A common experimental objective is to maximize an output of the simulator, i.e., to find an input x max that maximizes the output y(x). For minimization instead, replace y(x) by -y(x). Optimisation is usually a sequential process where successive xs are chosen to get closer and closer to the optimal x max -a sequential design problem (see Section 4). With stochastic simulators, y(x) is random, and optima are less concretely defined -the output is different every time the simulator is run at the same x. As a consequence, interest usually lies in maximizing a non-random quantity of interest, such as the mean, M , or possibly another scalar quantity such as the q th quantile.

For deterministic simulators Bayesian optimization [START_REF] Mockus | The application of Bayesian methods for seeking the extremum[END_REF][START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF] is a popular technique. An initial set of runs is used to build a GP surrogate and new runs are chosen by maximizing an "acquisition function" α(x). Iteratively choosing x new = arg max x α(x)) provides a progressively improved estimate for the maximum. A widely used choice for α(x) is the expected improvement (EI):

α EI (x) = E[max (y(x) -y max , 0)].

(6.3)

Maximizing EI chooses the input x new that maximizes the expected increase in the maximum value, y max , of already observed runs. With y modeled by a GP:

α EI (x) = (y max -µ N (x))Φ µ N (x) -y max σ N (x) + σ N (x)φ y max -µ N (x) σ N (x) (6.4)

where µ N (x) is the predictive mean of the GP, σ N (x) its standard deviation, φ is the standard normal density, and Φ the standard normal distribution function.

Alternative acquisition functions have generated extensive work on Bayesian optimization in recent years, mostly in the machine learning literature. The probability of improvement [START_REF] Kushner | A new method of locating the maximum point of an arbitrary multipeak curve in the presence of noise[END_REF] is an early example, and others, such as the GP upper confidence bound (GP-UCB) [START_REF] Srinivas | Gaussian process optimization in the bandit setting: No regret and experimental design[END_REF], consider homoscedastic simulator error. A recent summary can be found in [START_REF] Frazier | A tutorial on Bayesian optimization[END_REF]."

For stochastic simulators, the EI procedure can be extended by replacing y max , now a random variable, with the maximum estimated mean of currently run inputs, µ max = max i∈{1,...,N } µ N (x i ), see [START_REF] Vazquez | Global optimization based on noisy evaluations: an empirical study of two statistical approaches[END_REF]. Or, one can seek improvement over the maximum estimated mean of any possible input, max x µ N (x) [START_REF] Gramacy | Optimization under unknown constraints[END_REF]. In these cases, the σ N (x) term must exclude the σ 2 v (x) term that comes from say, a hetGP. Implementation of this method is provided in the hetGP package.

Alternative criteria for stochastic problems with constant intrinsic noise are discussed and compared in [START_REF] Picheny | A benchmark of kriging-based infill criteria for noisy optimization[END_REF]; with the above method is referred to as the "plugin" method. An R package for implementing several of these choices is available in DiceOptim [START_REF] Picheny | DiceOptim: Kriging-Based Optimization for Computer Experiments[END_REF]; [START_REF] Picheny | Noisy kriging-based optimization methods: a unified implementation within the diceoptim package[END_REF]). [START_REF] Jalali | Comparison of kriging-based methods for simulation optimization with heterogeneous noise[END_REF] also do a similar comparison for heteroscedastic noise.

The related goal of level set estimation to find regions where the output exceeds a threshold T can also be targeted with sequential criteria similar to EI. A simple criterion is maximum contour uncertainty (MCU), wherein new points are chosen according to a weighted sum of how close to T they are believed to be and the degree of uncertainty for that point. [START_REF] Lyu | Evaluating Gaussian process metamodels and sequential designs for noisy level set estimation[END_REF] provide some discussion here. This method is also implemented in hetGP.

Optimization using Gaussian processes, specifically in the presence of intrinsic variability that is potentially heteroscedastic (and potentially non-normal) is an interesting research question and possibly deserving of its own review. Nonetheless, the references provided here should provide a good introduction.

Sensitivity Analysis

Determining and measuring the effect of inputs on the output is usually part of any simulator experiment. Doing so assists scientific understanding of the system and enables screening out potentially superfluous variables. This goal has many related names: sensitivity analysis, screening, variable selection, etc., but the overall objective is generally the same -summarize and measure the influence of each input.

For deterministic simulators, Sobol indices [START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF] are widely used. Probabilistic distributions are assumed on the inputs of the simulator in order to represent their range of variations. Then, a functional Analysis of Variance (ANOVA) decomposition splits the variation of the simulator output into multiple components, each representing the individual contribution of an input variable x j or combination of input variables. A Sobol index is then computed as the percentage of the total simulator output variation explained by a component. Key Sobol indices include main effects (the percentage of variation explained by the individual x j s alone) and variation explained by interactive additive effects with other inputs. Computing the components takes large numbers of runs but the use of surrogate GPs make the calculations feasible [START_REF] Schonlau | Screening the input variables to a computer model via analysis of variance and visualization[END_REF][START_REF] Marrel | Calculations of Sobol indices for the Gaussian process metamodel[END_REF]. An enveloping discussion of sensitivity is provided by [START_REF] Oakley | Probabilistic sensitivity analysis of complex models: a Bayesian approach[END_REF].

Two extensions, by [START_REF] Marrel | Global sensitivity analysis of stochastic computer models with joint metamodels[END_REF] and [START_REF] Hart | Efficient computation of indices for stochastic models[END_REF], of Sobol indices for stochastic simulators yield the following expression for the stochastic simulator: y(x) = y(x, seed ) (6.5) where x is the set of controllable inputs. The input seed is responsible for output stochasticity, standing in for intrinsic variability, and is sometimes called a seed vari-able. As with a deterministic simulator, a probabilistic distribution (typically uniform) is assumed to represent the range of variation in controllable inputs.

In [START_REF] Marrel | Global sensitivity analysis of stochastic computer models with joint metamodels[END_REF], the total variation in the mean of the stochastic simulator is analysed through a functional ANOVA decomposition and Sobol indices are computed based on the percentage of the total simulator variation each component explains. The variation explained by the seed variable seed can also be computed, representing the total variation explained by the intrinsic variance. Additionally, a sensitivity analysis of the intrinsic variance σ 2 N (x) can be conducted separately to gather information on which input variables most impact the heteroscedasticity. [START_REF] Hart | Efficient computation of indices for stochastic models[END_REF] assumes the simulator can be run at different inputs x with the same seed seed . Rather than building a joint stochastic simulator surrogate for the mean and variance, as described in Section 3.2, they build a separate surrogate for a number of seeds. For each seed, they obtain a realization of each Sobol index, and by aggregating the realizations, they obtain distributions for the indices.

The extensive literature on model selection may have counterparts that can be effective for stochastic simulators. But a fully satisfactory approach even for deterministic simulators remains somewhat elusive.

Concluding Remarks

There are several key messages to be drawn from this review, each pointing to open or new research questions:

Gaussian Process Surrogates. GPs are discussed extensively because they provide a flexible way of allowing the data to inform about the shape of the underlying process. Moreover, they can be effective predictors and quantifiers of uncertainty. Diagnosing shortcomings in a GP for stochastic simulators (available in deterministic settings [START_REF] Bastos | Diagnostics for Gaussian process emulators[END_REF] is not yet well-established.

As noted in Section 3, neural network (deep learning) methods are in active use and under study, some of which may, in combination with GPs, offer promising research directions [START_REF] Schultz | Practical Bayesian optimization for transportation simulators[END_REF].

Additionally, it can be difficult to effectively capture non-normal variability. Doing so with as few simulations as possible, whilst also properly quantifying the various uncertainties, is likely to be an important research direction for stochastic simulator analysis. The wider quantile regression literature is likely a good starting point.

Design. Stochastic simulators differ from deterministic ones because they require much larger sample sizes and permit the use of replicates, whose treatment is generally ad hoc. This leads to the questions raised in Section 4, forming a direction of important research. Design size rules of thumb, useful even if imperfect, exist for deterministic simulators (Loeppky et al., 2009b), but are lacking for stochastic simulators.

Calibration. Accounting for model discrepancy in calibration is critical but there is no obvious "one-size-fits-all" method. A broad empirical comparison is needed with guidance about which strategies are effective under which conditions. Assessing the effectiveness of different methods can be challenging (see McKinley et al., 2018, for one comparison between ABC and HM), but sorely needed.

Simulator Complexity For complex stochastic simulators it may not be feasible to obtain enough runs. In some instances, the simulator can be replaced with a less complex one (e.g. [START_REF] Molina | Statistical inverse analysis for a network microsimulator[END_REF] that captures key features and permits adequate numbers of simulations. Another path, coupling stochastic simulators with deterministic simulators has been explored [START_REF] Baker | Predicting the output from a stochastic computer model when a deterministic approximation is available[END_REF] as a way to deal with low simulation budgets. Multi-fidelity modeling, where multiple simulators of varying complexity are coupled together [START_REF] Kennedy | Predicting the output from a complex computer code when fast approximations are available[END_REF][START_REF] Kennedy | Multilevel emulation for stochastic computer models with an application to large offshore windfarms[END_REF]) is a promising solution where possible.

In a similar vein, certain outputs may be less noisy than others, and the modeling of the less-noisy outputs can improve the modeling of the noisier ones. For example [START_REF] Wang | Enhancing response predictions with a joint Gaussian process model for stochastic simulation models[END_REF] use the expectation of a simulator to improve the estimation of noisier quantiles. This is related to the wider variance reduction literature, which has a long history [START_REF] Barton | History of improving statistical efficiency[END_REF]. Variance reduction has been applied in a number of examples but its use in ABMs is not apparent, perhaps due to the profusion of stochastic elements in an ABM. Fixing the initial seed in a stochastic simulator has played a role in sensitivity analysis (see Section 6.4), but leveraging information about the intrinsic randomness for wider purposes is an open problem. This review strives to raise awareness of existing tools and strategies for treating stochastic simulators and provide a starting point for practitioners interested in utilizing up-to-date statistical approaches. Despite the problems being pervasive and challenging there is a shortage of statistical research in this field. The problems pose computational and technical questions, as well as theoretical and philosophical ones. Current solutions are often capable, but there is a lack of comprehensive comparison between different solutions, and a lack of testing regarding their generalizability to complex situations (such as very large data sets). The hope is that the review provokes statistical researchers to engage the open questions discussed.

Figure 1 :

 1 Figure1: Fish example: 400 simulations consisting of 20 replicates at each of 20 population sizes (a maximin latin hypercube scaled to[150, 4000], rounded down to nearest integers). The left panels use homGP -the solid red line is the median of the predictive distribution and the dashed red-lines form the 95% uncertainty intervals. The right panels use hetGP. The upper panels include the data used to fit the surrogates; the lower panels omit the data but include the "true" values in black.

Figure 2 :

 2 Figure 2: Predictive mean and standard deviation surfaces for Ocean output using homGP and hetGP. Data are 1000 simulator runs consisting of 20 replicates at each of 50 input locations from a maximin Latin Hypercube Design (defined in Section 4) of size 50 in 2 dimensions. The top row provides predictive means, µ N , and the bottom gives standard deviations, σ N , of the predictive distribution of oxygen concentration. The left column uses homGP, the right uses hetGP.

Figure 3 :

 3 Figure3: Same setup as Figure1but with a QK surrogate. Mean predictions of 5 quantiles (5%, 25%, 50%, 75% and 95%) are provided on the left along with data, and mean predictions of 3 quantiles (5%, 50%, and 95%) are provided on the right along with the "true" values.

Figure 4 :

 4 Figure 4: Ocean Prediction with seqhetGP. Left plot is the mean, µ N , of the predictive distribution; right plot the standard deviation, σ N . Design sites with their replicates are superimposed on the right-hand plot.

Figure 5 :

 5 Figure 5: Performance results of the three ocean model surrogate model fits, repeated 100 times: the boxplots are for RMSE and Score computed, for each repetition, at the 500 test locations.

Figure 6 :

 6 Figure 6: Calibration results for the ocean model. The top row shows the posterior densities for the four parameters using homGP with the fixed design, the middle row uses hetGP with the fixed design, and the bottom row uses the hetGP surrogate with the sequential design. The budget for all three is 5000 runs. True values are superimposed as red vertical lines.

Figure 7 :

 7 Figure7: ABC fish calibration: directly (left) and via hetGP surrogate (right). The prior for the true population size is uniform on {200, . . . , 4000}, and 25 recaptured fish were observed. The left plot shows 10,000 simulations, highlights the 52 in agreement with the observation, and the histogram of accepted simulations. The histogram in the right panel is for the 3811 accepted draws out of 1,000,000 from hetGP.

Table 1 :

 1 Performance results of the three ocean model surrogates under KOH calibration. RMSE at the 500 test locations with the "true" values used for Figure

			( Kx , Ky )	OLS SINGLE NOCAL KOH
		homGP (824.9, 295.4) 9.16	9.16	9.22	9.14
	RMSE	hetGP	(754.9, 295.8) 9.16	9.16	9.20	9.15
		seqhetGP (496.3 276.0) 9.15	9.15	9.22	9.16
		homGP		-2.50	-2.66	-2.59	-2.32
	Score	hetGP		-2.55	-2.71	-2.62	-2.32
		seqhetGP		-2.55	-2.69	-2.61	-2.30

This terminology can have different meanings and connotations in different fields. In weather modeling, a stochastic simulator might refer specifically to a random weather generator[START_REF] Richardson | Stochastic simulation of daily precipitation, temperature, and solar radiation[END_REF][START_REF] Peleg | An advanced stochastic weather generator for simulating 2-d high-resolution climate variables[END_REF]

2017). In this work, we use the term to refer to any code that includes pseudorandom deviates in generating output.

If a large portion of the predictive distribution was negative in the transformed space, the untransformed intervals would be invalid, but this doesn't appear to be a problem in our example. Monotonic transforms exist to avoid this problem[START_REF] Johnson | Phenomenological forecasting of disease incidence using heteroskedastic Gaussian processes: a dengue case study[END_REF]. Predictions in the transformed space are also provided in the supplementary material.

The hetGP package also implements a Student-t variant[START_REF] Wang | Extended t-process regression models[END_REF][START_REF] Shah | Student-t processes as alternatives to Gaussian processes[END_REF][START_REF] Chung | Parameter and uncertainty estimation for dynamical systems using surrogate stochastic processes[END_REF].

With σ 2 N (x), the predictive variance, the IMSPE, of a design D is equal to x∈X σ 2 N (x) dx.

A Latin hypercube design is one where: on each dimension, the input space is divided into, usually, equal intervals and each interval is constrained to contain exactly one data point.

Such maximin LHDs are purportedly produced for example, by the maximinSLHD function of the R package SLHD[START_REF] Ba | SLHD: Maximin-Distance (Sliced) Latin Hypercube Designs[END_REF], or the lhs function from the Python package pyDOE[START_REF] Lee | pyDOE: The experimental design package for python[END_REF].

In R, the sobol function in the R package randtoolbox (Yohan[START_REF] Chalabi | randtoolbox: Toolbox for Pseudo and Quasi Random Number Generation and Random Generator Tests[END_REF]) can be used to generate Sobol sequences.

Acknowledgements

We gratefully acknowledge the support and funding provided by SAMSI during the Model Uncertainty: Mathematical and Statistical program, 2018-19. Moreover, Professor David Banks, Director of SAMSI, is thanked for suggesting and encouraging the writing of this review.

Robert Gramacy was supported in part by DOE LAB 17-1697 via a subaward from Argonne National Laboratory for SciDAC/DOE Office of Science ASCR and High Energy Physics, and by the National Science Foundation DMS-1821258.

Pierre Barbillon has received support from the Marie-Curie FP7 COFUND People Programme of the European Union, through the award of an AgreenSkills/AgreenSkills+ fellowship (under grant agreement n • 609398).

A Appendix -Ocean Truth

Throughout, reference to plots of the "truth" of the Ocean model is made. These plots are presented here, as well as in the supplementary material, for convenience.