<figure><img src="image_1.png" /></figure>

HAL
open science

# Gearing Gaussian process modeling and sequential design towards stochastic simulators

Mickael Binois, Arindam Fadikar, Abby Stevens

<figure><img src="image_2.png" /></figure>

## To cite this version:

Mickael Binois, Arindam Fadikar, Abby Stevens. Gearing Gaussian process modeling and sequential design towards stochastic simulators. 2024. hal-04824305

HAL Id: hal-04824305

https://inria.hal.science/hal-04824305v1

Preprint submitted on 6 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Gearing Gaussian process modeling and sequential design towards noisy simulators

Mickaël Binois*

Arindam Fadikar†

Abby Stevens†

## Abstract

This chapter presents specific aspects of Gaussian process modeling in the presence of complex noise. Starting from the standard homoscedastic model, various generalizations from the literature are presented: input varying noise variance, non-Gaussian noise, or quantile modeling. These approaches are compared in terms of goal, data availability and inference procedure. A distinction is made between methods depending on their handling of repeated observations at the same location, also called replication. The chapter concludes with the corresponding adaptations of the sequential design procedures. These are illustrated in an example from epidemiology.

## 1 Introduction

Accurately reproducing real-world dynamics often requires stochastic simulators, particularly in fields like epidemiology, operations research, and hyperparameter tuning. In these contexts it becomes important to distinguish between aleatoric uncertainty – arising from noise in observations, from epistemic uncertainty – stemming from uncertainty in the model. The former is sometimes called intrinsic uncertainty while the latter is referred to as extrinsic uncertainty, see e.g., Ankenman et al. (2010).

Gaussian process (GP) based surrogate methods (see, e.g., Rasmussen and Williams (2006); Gramacy (2020)) can be easily adapted from deterministic to noisy settings while maintaining strong predictive power, computational efficiency, and analytical tractability. Even in the deterministic setup, it is common to add a small diagonal nugget (also known as a jitter) term to the covariance matrix of the GP equations to ease its numerical inversion. It is also interpreted as a regularization term, especially in the reproducing kernel Hilbert space (RKHS) context, see, e.g., Kanagawa et al. (2018). This can be contrasted to the use of pseudo-inverses, which reverts to interpolation, see for instance the discussion by Mohammadi et al. (2016). Here we will prefer the term noise variance to relate it to intrinsic uncertainty, and also because the nugget effect has a different meaning in the kriging literature (see e.g., Roustant et al. (2012)).

<footnote>*Inria, Université Côte d'Azur, CNRS, LJAD, Sophia Antipolis, France mickael.binois@inria.fr</footnote>

<footnote>†Argonne National Laboratory, Decision and Infrastructure Sciences, U.S.A.</footnote>
Surrogates are often developed to replace calls to expensive simulators, and are hence used in a variety of contexts: uncertainty propagation, sensitivity analysis, Bayesian optimization, calibration, or inversion, see e.g., Gramacy (2020). Typically, from an initial set of observations, additional ones are added sequentially based on a suitable criterion. While relying on deterministic tools for these is possible, ignoring the noise variance may be detrimental to understanding complex dynamics and further adaptations are generally needed. For example, in both modeling and sequential design contexts, experiments are often replicated (i.e., different outputs are repeatedly generated from the same inputs due to the stochasticity). If this should be prevented in the noiseless context, it is an efficient tool to learn the noise features. Savings in total evaluation budget and computational speed may be obtained by adaptively balancing replication in the input space.

The goal of this chapter is to provide detailed treatments for situations when the intrinsic uncertainty becomes significant, resulting in low signal-to-noise ratios. We address both the modeling aspects and then the sequential design strategies under these conditions. In particular, we focus on alternatives for modeling the corresponding noise process, with various degrees of complexity and depending on data availability, with or without replication. While detailing only a few of the possible models and associated references, we refer for instance to Kleijnen (2015); Gramacy (2020); Baker et al. (2022) for complementary and more exhaustive viewpoints on the options. Sequential design procedures to enrich an initial data set are detailed, again taking into account the option to replicate.

This chapter is organized as follows: Section 2 details the simple adaptations for noisy modeling, disambiguating different contexts while keeping analytical tractability. Section 3 discusses further refinements requiring approximations. Section 4 details some modifications towards sequential designs. Finally, Section 5 summarizes the main concepts and opens to remaining challenges.

## 1.1 Related frameworks

For numerical simulators, the aleatoric behavior is generally controlled by a random seed argument. Given this seed information, seen as a categorical variable, the simulator becomes deterministic. As a result, dedicated covariance functions can be defined, see for instance Chen et al. (2012); Pearce et al. (2022) for common random numbers (CRN) oriented ones, or Zhou et al. (2011); Roustant et al. (2020); Zhang et al. (2020); Deshwal et al. (2021) for the more general options.

A similar context is when additional environmental variables are present: they correspond to variables that may not be controlled but whose values are known. Given the control and environmental variables, the output value remains deterministic. An additional related instance is when realizations of the random noise process are available, which may be used to estimate the behavior of the environmental variables, see e.g., El Amri et al. (2023). This pathwise or trajectory interpretation is further explored and exploited, e.g., by Wilson et al. (2021); Fadikar et al. (2023).

Considering input noise is an alternative to replace the environmental variables effect, with a given probability distribution on the inputs. For further analysis, a robust
GP is defined which integrates out the environmental variables or input noise. Without further assumptions, Monte Carlo integration is required to do so, while analytical expressions are available for simpler setups, typically with Gaussian noise, see e.g., Girard et al. (2002); McHutchon and Rasmussen (2011); Janusevskis and Le Riche (2013); Qing et al. (2022). More generally, GPs may also be defined directly on distributions, see e.g., Bachoc et al. (2020).

Notice that these differ from the i.i.d. noise assumption that we make with noisy GPs, with no extra information to exploit (irreducible noise). These setups, however, may be combined, e.g., having both input and output noise.

# 2 Noisy GP modeling

While the modification to account for noise in Gaussian processes is relatively straightforward, inference for an unknown input-varying noise variance function comes with some challenges. Options are discussed and illustrated on a data set.

## 2.1 Extending the constant noise framework

Denoting $Y$ the GP model of the simulator, recall the classical GP predictive equations for $N$ observations with i.i.d. additive Gaussian noise, $Y_i(\mathbf{x}_i) = f(\mathbf{x}_i) + \varepsilon_i$, $1 \le i \le N$, where $\varepsilon_i \sim N(0, r(\mathbf{x}_i))$ with $r(\cdot)$ is the (positive) variance function:

$$m_N(\mathbf{x}) = \mathrm{E}[Y(\mathbf{x})|y_1, \dots, y_N] = \mathbf{c}(\mathbf{x})^{\top}(\mathbf{C}_N + \mathbf{\Lambda}_N)^{-1}\mathbf{y} \quad (1)$$

$$s_N^2(\mathbf{x}) = \mathrm{Var}[Y(\mathbf{x})|y_1, \dots, y_N] = r(\mathbf{x}) + \sigma^2 (c(\mathbf{x}, \mathbf{x}) - c(\mathbf{x})^{\top}(\mathbf{C}_N + \mathbf{\Lambda}_N)^{-1}c(\mathbf{x})) \quad (2)$$

with $\mathbf{\Lambda}_N = \sigma^{-2}\mathrm{Diag}(r(\mathbf{x}_1), \dots, r(\mathbf{x}_N))$.

Compared to the standard homoscedastic GP model, with constant noise, that is $\mathbf{\Lambda}_N = \nu I$, here the noise variance can vary with $\mathbf{x}$ while keeping the analytical tractability of the GP prediction. Note that the predictive variance includes $r(\mathbf{x})$, removing which corresponds to $\mathrm{Var}[f(\mathbf{x})|y_1, \dots, y_N]$, that is the model for the underlying deterministic process $f$. Interestingly, as shown by van Beek et al. (2021), the predictive variance can also be decomposed into interpolation error and correction due to noise:

$$s_N^2(\mathbf{x}) = r(\mathbf{x}) + \sigma^2 \left( c(\mathbf{x}, \mathbf{x}) - c(\mathbf{x})^{\top} \mathbf{C}_N^{-1} c(\mathbf{x}) + \sum_{i=1}^n c(\mathbf{x})^{\top} \zeta_i c(\mathbf{x})^{\top} \right),$$

with details on the computation of $\zeta_i$ in the same article. In terms of theoretical analysis for the prediction error and convergence, some existing results with constant noise variance can be found, e.g., in Lederer et al. (2019); Garnett (2023).

So far this assumes that the noise function $r$ is known, which is seldom the case in practice. One exception is with Monte Carlo error, see e.g., Picheny et al. (2013a). Another appears in the context of hyperparameter tuning, where the accuracy of the output is controlled by an additional parameter. This may be the size of the data set on which a given machine learning model is trained, see e.g., Klein et al. (2015).
## 2.2 Inference for heteroscedastic GP modeling

With the variance process $r(\mathbf{x})$ seldom known, a variety of options have been proposed to learn it from the data. Increasing slightly in complexity, a parametric model may be used, such as $r(\mathbf{x}) = \exp(h(\mathbf{x}))$ with $h(\mathbf{x})$ from a parametric family, for example, simple polynomials, see e.g., Boukouvalas et al. (2014) or exponential family Le et al. (2005).

For more flexibility, many works entertain a second GP to model the log variance, ensuring positive variance at the same time, starting from Goldberg et al. (1998). The main caveat is that these log variances are unobserved, thus being latent variables that must be estimated as well. Goldberg et al. (1998) rely on a fully Bayesian approach and Markov chain Monte Carlo (MCMC) methods to do so. Later works took approximate inference methods, expectation-maximization (EM) Kersting et al. (2007); Boukouvalas and Cornford (2009), variational approximation Titsias and Lázaro-Gredilla (2011) or maximum likelihood Binois et al. (2018). Even if the learning is simplified by the underlying smoothing effect of using a GP on the log noise, the number of latent variables grows with the number of observations, unless sparse GPs are used, see e.g., Snelson and Ghahramani (2006). They rely on inducing points to summarize the entire distribution.

A perhaps simpler idea is to learn $r(\mathbf{x})$ from empirical estimates of the variance, which can only be computed if replication is present in the data. Denote $y_i^{(j)}, 1 \le j \le a_i$, $a_i$ repeated observations at a given $\mathbf{x}_i$, with mean $\bar{y}_i$, then $\hat{\sigma}_i^2 = \frac{1}{a_i-1} \sum_{j=1}^{a_i} (y_i^{(j)} - \bar{y}_i)^2$. The GP on the log noise variance can then be trained directly on these variance estimates, with the so-called stochastic kriging model Ankenman et al. (2010). With replication, among the total $N$ observations, there are only $n$ unique designs $\bar{\mathbf{x}}_i, N = \sum_{i=1}^n a_i$. An additional benefit is that the predictive equations (2) are equivalent, see e.g., Ankenman et al. (2010); Picheny et al. (2013a); Binois et al. (2018), to:

$$m_n(\mathbf{x}) = \mathbf{c}(\mathbf{x})^{\top}(\mathbf{C}_n + \mathbf{\Lambda}_n\mathbf{A}_n^{-1})^{-1}\bar{\mathbf{y}}, \quad (3)$$

$$s_n^2(\mathbf{x}) = r(\mathbf{x}) + \sigma^2 (c(\mathbf{x}, \mathbf{x}) - \mathbf{c}(\mathbf{x})^{\top}(\mathbf{C}_n + \mathbf{\Lambda}_n\mathbf{A}_n^{-1})^{-1}\mathbf{c}(\mathbf{x})) \quad (4)$$

where $\bar{\mathbf{y}} = (\bar{y}_1, \dots, \bar{y}_n)$, $\mathbf{\Lambda}_n = \sigma^{-2}\text{Diag}(r(\bar{\mathbf{x}}_1), \dots, r(\bar{\mathbf{x}}_n))$ and $\mathbf{A}_n = \text{Diag}(a_1, \dots, a_n)$. A direct benefit of these expressions is that the computational complexity is much reduced, depending on $n$ rather than $N$. In terms of error analysis, results are available for stochastic kriging, e.g., by Wang and Haaland (2019).

In between the latent variance and empirical approaches, hybrids have been proposed to leverage the reduced computational complexity coming with replicates while allowing fewer (even no) replicates per unique design, see e.g., Boukouvalas and Cornford (2009); Binois et al. (2018). This also raises the question of how to allocate replicates better than fixing the same amount at every location, which will be discussed in Section 4.
## 2.3 Illustration

We provide an illustrative example of a 1d SIR (Susceptible–Infected–Recovered) problem in Figure 1. For one data set, we use 25 unique designs each with 100 replicates while there are 2500 unique locations for the second data set. We use 10000 replicates at 51 unique locations for computing a reference. Already from this simple data set, it appears that the empirical estimation is not precise, even using 100 replicates. This even shows for the reference on the variance and skewness estimation. Without replication, estimating the variance or higher moments requires learning latent variables. Trajectories corresponding to the CRN setup, that is fixing the seed values, are also added for illustration of the difference with the iid setup.

<figure><img src="image_4.png" /><figcaption>Figure 1: SIR simulator data sets. Top left: 2500 observations with replications, allowing to compute the estimated mean and quantiles. Top right: 2500 observations without replications. Five realizations with the same seed value across the x-axis are also depicted. Bottom: the empirical (points) and reference (dashed lines) variance and skewness for the data set 1.</figcaption></figure>

Then we fit the first GP models on these two data sets, using the hetGP Binois and Gramacy (2021) R R Core Team (2024) package, with results in Figure 2. The predictive means on the 5 trajectories are also shown. From this relatively simple example, one can observe the ability of a heteroscedastic model to better represent the black-box at hand, in terms of the first two moments. As for the 90% predictive intervals, since the predictions remain Gaussian, such a model cannot capture higher moments accurately.
As will be illustrated throughout the entire chapter, capturing finer features requires more data. Even a GP with constant noise is sufficient to estimate the mean response, then heteroscedastic GP (hetGP) models also captures the noise variance, but even more complex models are required to go beyond. The difference between the outcomes for data sets 1 and 2 are slight, but the training times are very different: 0.08 vs 20.2s for the homoscedastic model, 0.093 vs 396s for the heteroscedastic one. A stochastic kriging model would need fitting two GPs, and is only possible on data set 1. From the estimation of variances, a benefit of the stochastic kriging model is the availability of an estimate of the predictive variance on the log variance. This is also available depending on the inference procedure with hetGP models. The use of replicates is beneficial in terms of variance estimation, it better captures the initial values around $x \approx 0$ while hetGP models are in the right ballpark.

# 3 Extensions beyond the Gaussian case

For more complex noise structures, more complex models are needed, at the cost of less analytical tractability. Besides transforming the outputs with a parameterized monotonous function as proposed e.g., by De Oliveira et al. (1997); Snelson et al. (2004), we detail some of the variations. A discussion from the RKHS regularization viewpoint is proposed e.g., by Aravkin et al. (2015) to improve robustness to outliers.

## 3.1 Some options

### 3.1.1 Student-t variations and non-Gaussian likelihoods

A first modulation of the Gaussian noise assumption is with Student-t noise, with larger tails. Using Student-t processes, Shah et al. (2014); Wang et al. (2017), analytical expressions may be obtained but can revert back to GP ones when the number of observations grows, limiting the interest to smaller data sets.

A more popular assumption is to use a Gaussian process coupled with Student-t noise, see e.g., Vanhatalo et al. (2009); Hartmann and Vanhatalo (2019); Lyu et al. (2021). Nevertheless, the likelihood is no-longer Gaussian and the posterior distribution cannot be obtained in closed form. It can still be approximated, using for instance the Laplace approximation of the posterior, as detailed e.g., by Rasmussen and Williams (2006). This later approximation is also used with Poisson distribution for count data, e.g., by Flaxman et al. (2015). Another typical use if for classification with Gaussian processes, with the Bernoulli distribution.

A wider range of likelihood functions is possible when using MCMC for approximating the posterior function, made computationally feasible by relying on variational inducing points, as proposed, e.g., by Flaxman et al. (2015). Coupling variational inducing points and variational inference for non-Gaussian likelihoods is the approach proposed e.g., by Hensman et al. (2015b,a) and available for instance in the GPflow Matthews et al. (2017), GPyTorch Gardner et al. (2018) or GPJax Pinder and Dodd (2022) Python packages. The point of the variational framework is to approximate the intractable posterior distribution by a simpler parameterized one. These parameters
<figure><img src="image_5.png" /><figcaption>Constant noise model, data set 2</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Noise varying model, data set 2</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Data set 1</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Data set 2</figcaption></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>Figure 2: Illustration on a 1d SIR simulator. Homoscedastic versus heteroscedastic models predictions on the two data sets. Bottom: corresponding variance predictions.</figcaption></figure>
are jointly optimized with the other model parameters by minimizing the Kullback-Leibler divergence between approximate and true posterior via the so-called evidence lower bound (ELBO). Combined with the latent variables idea from Section 2, one can model non-Gaussian heteroscedastic models, as proposed by Saul et al. (2016) with Student-t or even beta distributed noises. In terms of flexibility, the generalized lambda distribution, with its four parameters, can also model a variety of non-Gaussian distributions and be cast within the same variational framework. This distribution has also been used with polynomial chaos expansions latent functions by Zhu and Sudret (2021). They offer a peek at the underlying phenomenon with the ability to learn the latent variables.

It remains that these recent generalization capabilities are enabled by several layers of approximations. Taking a step back, if the goal is really to learn the tails of unknown noise distributions, then perhaps looking directly at quantiles (or expectiles) is sufficient, and it as been proposed, e.g., by Picheny et al. (2022) with two latent GPs.

### 3.1.2 Quantile GP modeling

Instead of assuming i.i.d. additive noise on $ε_i$, Plumlee and Tuo (2014) proposed reconstructing replicate variability by estimating quantiles as a function of the inputs. This approach relies on two key principles: (1) with a sufficient number of quantiles, the empirical CDF can be represented as a mixture of point masses, allowing for fast sampling; and (2) quantiles are assumed to be continuous, making them suitable for modeling with a Gaussian Process (GP). The idea is to fit a GP directly to the empirical quantiles, accounting for estimation error, without assuming Gaussian noise. Alternatively, the training data can be modified by estimating multiple quantiles from noisy simulations and augmenting the corresponding quantile levels into the input. Applying a standard GP to this modified data provides similar kriging estimates while avoiding distributional assumptions on the noise. The same formulations in (2) can be modified as:

$$Y_i(\mathbf{x}_i, \alpha) = f(\mathbf{x}_i, \alpha) + \varepsilon_i, \quad \varepsilon_i \sim N(0, \sigma^2), \quad 1 \le i \le N, \qquad (5)$$

where, $0 < \alpha < 1$ is the quantile level. For instance, Picheny et al. (2013a) also considered the GP predictive quantiles, which are available in closed form. If the level of the quantile is not decided beforehand, then it can be an additional parameter to the approach, as entertained by Fadikar et al. (2018). Taking a loss function tailored for quantiles and latent input noise variables, Quadrianto et al. (2009) propose a quantile regression framework based on GPs. It is further extended by Picheny et al. (2022) relying on a variational framework.

### 3.1.3 Deep GPs

Noisy simulations can alternatively be viewed through the lens of non-stationarity. By applying non-linear transformations to a GP, we can create a stochastic process with input-dependent noise variance, albeit without explicit control over the noise structure. We focus on warped GPs, where input and/or output are transformed to model complex data. Warping functions, such as the monotonic neural-net sum
of tanh functions in Snelson et al. (2004), add flexibility while maintaining the full probabilistic framework for inference. This model can be further improved by placing a non-parametric GP prior on the warping function itself, as described in (Lázaro-Gredilla, 2012), where the inference is performed using variational methods. Deep GPs are generalizations of the warped GPs. A deep GP is composed of multiple hidden layers of latent variables, where each node serves as the output for the layer above and the input for the layer below. The observed outputs are positioned at the bottom layer, while Gaussian processes control the mappings between layers. A single-layer deep GP (also known as GP-Latent Variable Model or GP-LVM) was first proposed by Lawrence (2003) and a hierarchical version of GP-LVM was introduced in Lawrence and Moore (2007). Damianou and Lawrence (2013) later developed the deep GP framework and the variational framework to fully infer the model parameters for any number of layers. Currently, a number of Python and R packages exist for fitting a deep GP model, like for instance GPyTorch Gardner et al. (2018), GPflux Dutordoir et al. (2021) or deepgp Booth (2024).

## 3.2 Illustration

Here we focus on directly modeling the quantiles of the data set, relying either on replicates or latent variables. The outcomes are represented in Figure 3. For the latent quantile GP, we rely on Picheny et al. (2022). On data set 1, the results of both approaches are similar, with broader prediction intervals for the GPs on empirical quantiles. Without replication in data set 2, only the latent quantile model can be applied, which manages to fit the actual quantiles well, except at the origin. Replicates with no or small variance can force the model to take it into account. We fit three layers deep GP models to both data sets using deepgp package in R. However, the model predictions are not significantly different between the two data sets. Figure 3 shows the model fit on the data set with replication.

After presenting broad options of modeling capabilities, the specific choice for a given application would take into account various aspects. First, the amount of available data will drive the choice, the most complex models requiring more data. Using an appropriate noise model is also important to accurately represent the data, especially if this noise structure is known. Finally, the structure of the design of experiments may also restrict options: with replication, estimating moments is facilitated for modeling them directly. Otherwise, inducing points and variational inference of latent quantities may be necessary to keep a reasonable computational cost.

Next we will review the options to improve the model for specific targets.

# 4 Adapted sequential design strategies

While noise modeling brings more complexity, it also introduces a regularization component which may actually alleviate numerical issues (covariance matrix conditioning) and sequential design procedures (the acquisition function value is smoother as it is not zero at evaluated designs).
<figure><img src="image_11.png" /><figcaption>Figure 3: Quantile GP models: directly on empirical quantiles (top-left) and using latent variables (right). Reference quantiles are in shaded blue. Deep GP model on data set 1 (bottom-left).</figcaption></figure>

Before delving into specific methods for various targets, some general techniques are possible, see for instance Jalali et al. (2017); Rojas-Gonzalez and Van Nieuwenhuyse (2020). A first naive one is simply to ignore the noise and do exactly as in the noiseless case. This is probably fine as long as the signal to noise ratio is high. Increasing this latter is achievable by fixing a number of replicates for each evaluation. A more careful approach is to filter out the noise, by using the prediction by the noisy GP model as inputs for a deterministic GP employed in the sequential approach. This is called reinterpolation, see, e.g., Forrester et al. (2008).

We continue by detailing several typical applications: global accuracy and optimization. Estimating a probability of failure, see e.g., Hao et al. (2021), or calibration Baker et al. (2022), or a level set, see e.g., Lyu et al. (2021), are other related tasks, and so is sensitivity analysis Mazo (2021). We refer to e.g., Kleijnen (2015); Santner et al. (2018); Ginsbourger (2018) and references therein for more general aspects.
## 4.1 Global accuracy

Broadly speaking, getting a globally accurate surrogate model amounts to have good space-filling properties, see e.g., Gramacy (2020). There is a reasonable criterion to adopt is the integrated mean square prediction error (IMSPE): $\alpha_{IMSPE}(\mathbf{x}) = \int_{\mathcal{X}} s_N^2(\mathbf{x}) d\mathbf{x}$. Optimizing it sequentially attempts to reduce the predictive variance everywhere, without putting too much effort on the frontier of the domain. Alternatively, targeting model estimation is studied for instance by Boukouvalas et al. (2014) with the Fisher information.

The heteroscedastic noise context brings a twist, to spread the effort between areas where the noise variance is more or less large. When the input set $\mathcal{X}$ is discrete, the optimal allocation of replicates over designs is given by Ankenman et al. (2010). For continuous inputs, Binois et al. (2019) show that replication may be preferred to adding new designs depending on $r(\mathbf{x})$ and propose further strategies to promote replication. A key idea is to look ahead, meaning taking into account the effect of future decisions. That is, reducing the number of alternatives to keep tractability, to assess whether or not evaluating a new design now and then replicating is better than waiting and replicate a few times before adding a new design. Since IMSPE only depends on the predictive variance, which does not depend on the observations, this is a favorable context for these less myopic options. Batch extensions are detailed e.g., in Zhang et al. (2022).

## 4.2 Optimization

Global accuracy criteria are usually related to the predictive variance or covariance matrix only. This is less simple for other goal oriented approaches such as those geared towards optimization – over the mean in this section – that require adaptations. Notable exceptions where the noisy setup is transparent include the probability of improvement, knowledge gradient or entropy based criteria, albeit these two latter families require approximations to be computed and an estimate of $r(\mathbf{x})$.

The difficulty for the popular expected improvement (EI) criterion is to define the current best value $T$: $\alpha_{EI}(\mathbf{x}) = E[\max(0, T - Y(\mathbf{x})|y_1, \dots, y_N)]$. Indeed, noise realizations may trick the surrogate model, leading to an overly optimistic estimation of this target. It will then require more samples for correction. This effect may be more pronounced for multi-objective optimization, where a set of Pareto optimal solutions are needed. Properly identifying the current best estimate from the surrogate model is called the identification step by Jalali et al. (2017), which is required at least at the final iteration to return a solution to the user.

For EI, a simple adaptation is to plug-in the predictive mean instead of function values, at the evaluated designs Vazquez et al. (2008) or over $\mathcal{X}$ Gramacy and Lee (2011). An offset may be added for more robust estimation of this plugin, combined to a penalty for the noise variance at $\mathbf{x}_{n+1}$ by Huang et al. (2006). With additional noisy constraints, a Bayesian version of EI is proposed by Letham et al. (2019). Thorough discussions available in Frazier (2018); Gramacy (2020); Garnett (2023) details relationships between EI, KG or entropy infill criteria, and extensions to related topics: non myopic, batch or multi-fidelity criteria. Last but not least, Astudillo and Frazier
(2019) insists on modeling the actual output directly, rather than modeling a sum, an integral or even a composition.

Criteria dedicated to the noisy context include the expected quantile improvement Picheny et al. (2010), looking at the improvement over a quantile. Notice that the upper confidence bound (UCB) criterion can be seen as directly optimizing a GP quantile for a given level. Benchmark of the above criteria are presented for instance by Picheny et al. (2013b); Jalali et al. (2017). Jalali et al. (2017) insist on the importance of a replication strategy, interleaving an additional stage to the usual procedure and on an appropriate identification strategy to actually return a good solution at the end.

So far only standard GPs, including the heteroscedastic version have been consid-
ered. In the case of Student-t processes, dedicated criteria have also been obtained
in closed form, for multi- van der Herten et al. (2016) or single-objective Xie and
Chen (2017) optimization. In the context of CRNs, we refer to Pearce et al. (2022)
for exploiting the relationship between seeds to obtain the best mean results and to
Fadikar et al. (2023) for finding good individual realizations (irrespective of the mean
behavior).

Only a few criteria exist to directly choose the number of replicates when evaluating, see e.g., Jalali et al. (2017); van Beek et al. (2021); Binois et al. (2021). But two stage approaches have been proposed in the operation research literature, first new points are evaluated, then the remaining budget is used to replicate at the most promising locations. Some rules have been derived for selecting and ranking solutions, sometimes under the term optimal computing budget allocation(OCBA) Chen et al. (2000); Quan et al. (2013), or for multi-objective optimization counterparts, e.g., Hunter et al. (2017); Rojas-Gonzalez and Van Nieuwenhuyse (2020); Gonzalez et al. (2020).

Another way to see this is to consider the cost associated with a noise variance
reduction, hence cost-aware strategies can be employed, see e.g., Garnett (2023). We
mention another simple strategy when the signal to noise ratio is low: enforcing a
minimal variance reduction for a new evaluation. This helps ensuring that there will
be sufficient information to actually modify the prediction significantly, hence reducing
the time for BO.

4.3 Robust optimization

A related topic when dealing with uncertainty is to edge against it by the formulation of the problem: rather than considering simply the expected behavior, a threshold on the variance may be used (or estimating a mean-variance Pareto front), a worst case approach may be adopted, optimizing a given quantile or taking a chance constrained framework. Nevertheless these are generally tackled considering environmental variables (necessary, e.g., for worst-case frameworks).

Using noisy surrogates, optimization on quantiles is tackled e.g., by Browne et al. (2016); Picheny et al. (2022), with a model of the quantiles and an EI-like infill criterion, or by Makarova et al. (2021) for a weighted sum of the mean plus variance. Additional criteria looking at the (conditional) value at risk are studied for instance by Chen and Kim (2016) within a stochastic kriging framework, when Picheny et al. (2022) uses latent GPs for directly modeling the quantiles with entropy and Thompson sampling
acquisition functions. Finally, a multi-objective version is described, e.g., by Rivier and Congedo (2022).

## 4.4 Illustration

We show two strategies: one selecting individual evaluation at every iteration, but targeting replicates, one choosing simultaneously the next design with its replication budget, corresponding to a variance reduction of at least 10%. Which one is preferable depends on the level of concurrency for evaluations, the signal to noise ratio or the time for single evaluations.

A first task is to reach global accuracy, hence we follow the approach relying on looking ahead with IMSPE to lightly enforce replication. For the second task, we consider a level set exercise to identify where $f = 0.5$. The contour stepwise uncertainty reduction (SUR) criterion proposed by Lyu et al. (2021) is used. Similar to the knowledge gradient idea, it focuses on the effect of a new observation on the estimation of the contour. Results for these two setups are presented in Figure 4. For IMSPE, a total of 900 points are selected sequentially, one per iteration, resulting in ≈ 120 unique designs. With the contour SUR, only ≈ 30 iterations are necessary to reach a total of 1000 evaluations. Both succeed in improving the initial estimate towards the goal, either placing replicates in a more or less space-filling way for IMSPE, except on the left where the noise variance is smaller while for the level set the replicates are on both sides of the crossing.

# 5 Conclusion and perspectives

In this chapter, we reviewed a variety of methods for modeling noisy simulators. Some benefit from or require replication, while others make no such assumption. The former may be simpler to interpret and re-implement, while the latter are available with advanced implementations and inference frameworks. Both frameworks have been adapted for sequential designs.

We focused mostly on GP models, but other random fields model may be relevant as well, see e.g., Kozubowski et al. (2013). Similar in structure to their neural networks counterparts, deep GPs use Gaussian processes for their layers, see e.g., Damianou and Lawrence (2013); Dunlop et al. (2018); Sauer et al. (2023a), and have been adapted to sequential learning as well Hebbal et al. (2021); Yazdi (2022); Sauer et al. (2023b) but, as discussed in Section 3.1.3, without dedicated treatment for the noise. Further studies of their performance in the noisy setup would be needed, perhaps with additional focus on the inference procedure.

In terms of remaining road blocks, most are shared with regular GPs, such as coping with non-stationarity or high-dimensional inputs. More specific ones in the noisy context include model selection depending on the amount of available data or more fine-tuned handling of replicates.
<figure><img src="image_12.png" /><figcaption>Figure 4: Sequential learning illustration. Left: initial state, right: state after using IMSPE (global accuracy, top) with empirical means in red when there are replicates. Bottom: state for 0.5 level set estimation. Gray vertical bars represent the amount of replication per design.</figcaption></figure>

References

Ankenman, B., Nelson, B. L., and Staum, J. (2010). Stochastic kriging for simulation metamodeling. *Operations research*, 58(2):371–382.

Aravkin, A. Y., Bell, B. M., Burke, J. V., and Pillonetto, G. (2015). The connection between Bayesian estimation of a Gaussian random field and RKHS. *IEEE transactions on neural networks and learning systems*, 26(7):1518–1524.

Astudillo, R. and Frazier, P. (2019). Bayesian optimization of composite functions. In *International Conference on Machine Learning*, pages 354–363. PMLR.

Bachoc, F., Suvorikova, A., Ginsbourger, D., Loubes, J.-M., and Spokoiny, V. (2020). Gaussian processes with multidimensional distribution inputs via optimal transport and Hilbertian embedding. *Electronic Journal of Statistics*, 14:2742–2772.

Baker, E., Barbillon, P., Fadikar, A., Gramacy, R. B., Herbei, R., Higdon, D., Huang,
J., Johnson, L. R., Ma, P., Mondal, A., et al. (2022). Analyzing stochastic computer models: A review with opportunities. *Statistical Science*, 37(1):64–89.

Binois, M., Collier, N., and Ozik, J. (2021). A portfolio approach to massively parallel Bayesian optimization. *arXiv preprint arXiv:2110.09334*.

Binois, M. and Gramacy, R. B. (2021). hetgp: Heteroskedastic Gaussian process modeling and sequential design in R. *Journal of Statistical Software*, 98(13):1–44.

Binois, M., Gramacy, R. B., and Ludkovski, M. (2018). Practical heteroscedastic Gaussian process modeling for large simulation experiments. *Journal of Computational and Graphical Statistics*, 27(4):808–821.

Binois, M., Huang, J., Gramacy, R. B., and Ludkovski, M. (2019). Replication or exploration? Sequential design for stochastic simulation experiments. *Technometrics*, 61(1):7–23.

Booth, A. S. (2024). *deepgp: Bayesian Deep Gaussian Processes using MCMC*. R package version 1.1.3.

Boukouvalas, A. and Cornford, D. (2009). Learning heteroscedastic Gaussian processes for complex datasets. *Group*, 44(0).

Boukouvalas, A., Cornford, D., and Stehlík, M. (2014). Optimal design for correlated processes with input-dependent noise. *Computational Statistics & Data Analysis*, 71:1088–1102.

Browne, T., Iooss, B., Gratiet, L. L., Lonchampt, J., and Remy, E. (2016). Stochastic simulators based optimization by Gaussian process metamodels–application to maintenance investments planning issues. *Quality and Reliability Engineering International*, 32(6):2067–2080.

Chen, C.-H., Lin, J., Yücesan, E., and Chick, S. E. (2000). Simulation budget allocation for further enhancing the efficiency of ordinal optimization. *Discrete Event Dynamic Systems*, 10:251–270.

Chen, X., Ankenman, B. E., and Nelson, B. L. (2012). The effects of common random numbers on stochastic kriging metamodels. *ACM Transactions on Modeling and Computer Simulation (TOMACS)*, 22(2):1–20.

Chen, X. and Kim, K.-K. (2016). Efficient VaR and CVaR measurement via stochastic kriging. *INFORMS Journal on Computing*, 28(4):629–644.

Damianou, A. and Lawrence, N. D. (2013). Deep Gaussian processes. In *Artificial intelligence and statistics*, pages 207–215. PMLR.

De Oliveira, V., Kedem, B., and Short, D. A. (1997). Bayesian prediction of transformed Gaussian random fields. *Journal of the American Statistical Association*, 92(440):1422–1433.
Deshwal, A., Belakaria, S., and Doppa, J. R. (2021). Bayesian optimization over hybrid spaces. In *International Conference on Machine Learning*, pages 2632–2643. PMLR.

Dunlop, M. M., Girolami, M. A., Stuart, A. M., and Teckentrup, A. L. (2018). How deep are deep Gaussian processes? *Journal of Machine Learning Research*, 19(54):1–46.

Dutordoir, V., Salimbeni, H., Hambro, E., McLeod, J., Leibfried, F., Artemev, A., van der Wilk, M., Hensman, J., Deisenroth, M. P., and John, S. (2021). GPflux: A library for deep Gaussian processes. *arXiv preprint arXiv:2104.05674*.

El Amri, M. R., Helbert, C., Zuniga, M. M., Prieur, C., and Sinoquet, D. (2023). Feasible set estimation under functional uncertainty by Gaussian process modelling. *Physica D: Nonlinear Phenomena*, 455:133893.

Fadikar, A., Collier, N., Stevens, A., Ozik, J., Binois, M., and Toh, K. B. (2023). Trajectory-oriented optimization of stochastic epidemiological models. In *2023 Winter Simulation Conference (WSC)*, pages 1244–1255. IEEE.

Fadikar, A., Higdon, D., Chen, J., Lewis, B., Venkatramanan, S., and Marathe, M. (2018). Calibrating a stochastic, agent-based model using quantile-based emulation. *SIAM/ASA Journal on Uncertainty Quantification*, 6(4):1685–1706.

Flaxman, S., Wilson, A., Neill, D., Nickisch, H., and Smola, A. (2015). Fast Kronecker inference in Gaussian processes with non-Gaussian likelihoods. In *International conference on machine learning*, pages 607–616. PMLR.

Forrester, A., Sobester, A., and Keane, A. (2008). *Engineering design via surrogate modelling: a practical guide*. John Wiley & Sons.

Frazier, P. I. (2018). Bayesian optimization. In *Recent advances in optimization and modeling of contemporary problems*, pages 255–278. Informs.

Gardner, J., Pleiss, G., Weinberger, K. Q., Bindel, D., and Wilson, A. G. (2018). Gpytorch: Blackbox matrix-matrix Gaussian process inference with gpu acceleration. *Advances in neural information processing systems*, 31.

Garnett, R. (2023). *Bayesian Optimization*. Cambridge University Press.

Ginsbourger, D. (2018). *Sequential Design of Computer Experiments*, pages 1–9. American Cancer Society.

Girard, A., Rasmussen, C., Candela, J. Q., and Murray-Smith, R. (2002). Gaussian process priors with uncertain inputs application to multiple-step ahead time series forecasting. *Advances in neural information processing systems*, 15.

Goldberg, P. W., Williams, C. K., and Bishop, C. M. (1998). Regression with input-dependent noise: A Gaussian process treatment. *Advances in neural information processing systems*, pages 493–499.
Gonzalez, S. R., Jalali, H., and Van Nieuwenhuyse, I. (2020). A multiobjective stochastic simulation optimization algorithm. *European Journal of Operational Research*, 284(1):212–226.

Gramacy, R. B. (2020). *Surrogates: Gaussian Process Modeling, Design, and Optimization for the Applied Sciences*. CRC Press.

Gramacy, R. B. and Lee, H. K. (2011). Optimization under unknown constraints. *Bayesian Statistics*, 9(9):229.

Hao, P., Feng, S., Liu, H., Wang, Y., Wang, B., and Wang, B. (2021). A novel nested stochastic kriging model for response noise quantification and reliability analysis. *Computer Methods in Applied Mechanics and Engineering*, 384:113941.

Hartmann, M. and Vanhatalo, J. (2019). Laplace approximation and natural gradient for Gaussian process regression with heteroscedastic student-t model. *Statistics and Computing*, 29(4):753–773.

Hebbal, A., Brevault, L., Balesdent, M., Talbi, E.-G., and Melab, N. (2021). Bayesian optimization using deep Gaussian processes with applications to aerospace system design. *Optimization and Engineering*, 22:321–361.

Hensman, J., Matthews, A., and Ghahramani, Z. (2015a). Scalable variational Gaussian process classification. In *Artificial Intelligence and Statistics*, pages 351–360. PMLR.

Hensman, J., Matthews, A. G., Filippone, M., and Ghahramani, Z. (2015b). MCMC for variationally sparse Gaussian processes. *Advances in Neural Information Processing Systems*, 28.

Huang, D., Allen, T. T., Notz, W. I., and Zeng, N. (2006). Global optimization of stochastic black-box systems via sequential kriging meta-models. *Journal of global optimization*, 34:441–466.

Hunter, S. R., Applegate, E. A., Arora, V., Chong, B., Cooper, K., Rincón-Guevara, O., and Vivas-Valencia, C. (2017). An introduction to multi-objective simulation optimization. *Optimization Online*.

Jalali, H., Van Nieuwenhuyse, I., and Picheny, V. (2017). Comparison of kriging-based methods for simulation optimization with heterogeneous noise. *European Journal of Operational Research*, 261(1):279–301.

Janusevskis, J. and Le Riche, R. (2013). Simultaneous kriging-based estimation and optimization of mean response. *Journal of Global Optimization*, 55:313–336.

Kanagawa, M., Hennig, P., Sejdinovic, D., and Sriperumbudur, B. K. (2018). Gaussian processes and kernel methods: A review on connections and equivalences. *arXiv preprint arXiv:1807.02582*.
Kersting, K., Plagemann, C., Pfaff, P., and Burgard, W. (2007). Most likely heteroscedastic Gaussian process regression. In *Proceedings of the 24th international conference on Machine learning*, pages 393–400. ACM.

Kleijnen, J. P. (2015). *Design and Analysis of Simulation Experiments*, volume 230. Springer.

Klein, A., Bartels, S., Falkner, S., Hennig, P., and Hutter, F. (2015). Towards efficient Bayesian optimization for big data. In *NIPS 2015 Bayesian Optimization Workshop*.

Kozubowski, T. J., Podgórski, K., and Rychlik, I. (2013). Multivariate generalized Laplace distribution and related random fields. *Journal of Multivariate Analysis*, 113:59–72.

Lawrence, N. (2003). Gaussian process latent variable models for visualisation of high dimensional data. *Advances in neural information processing systems*, 16.

Lawrence, N. D. and Moore, A. J. (2007). Hierarchical Gaussian process latent variable models. In *Proceedings of the 24th International Conference on Machine Learning*, ICML '07, page 481–488, New York, NY, USA. Association for Computing Machinery.

Lázaro-Gredilla, M. (2012). Bayesian warped Gaussian processes. *Advances in Neural Information Processing Systems*, 25.

Le, Q. V., Smola, A. J., and Canu, S. (2005). Heteroscedastic Gaussian process regression. In *Proceedings of the 22nd international conference on Machine learning*, pages 489–496. ACM.

Lederer, A., Umlauft, J., and Hirche, S. (2019). Uniform error bounds for Gaussian process regression with application to safe control. *Advances in Neural Information Processing Systems*, 32.

Letham, B., Karrer, B., Ottoni, G., and Bakshy, E. (2019). Constrained Bayesian optimization with noisy experiments. *Bayesian Analysis*, 14(2):495.

Lyu, X., Binois, M., and Ludkovski, M. (2021). Evaluating Gaussian process metamodels and sequential designs for noisy level set estimation. *Statistics and Computing*, 31(4):1–21.

Makarova, A., Usmanova, I., Bogunovic, I., and Krause, A. (2021). Risk-averse heteroscedastic Bayesian optimization. *Advances in Neural Information Processing Systems*, 34:17235–17245.

Matthews, A. G. d. G., Van Der Wilk, M., Nickson, T., Fujii, K., Boukouvalas, A., León-Villagrá, P., Ghahramani, Z., and Hensman, J. (2017). GPflow: A Gaussian process library using TensorFlow. *J. Mach. Learn. Res.*, 18(40):1–6.
Mazo, G. (2021). A trade-off between explorations and repetitions for estimators of two global sensitivity indices in stochastic models induced by probability measures. *SIAM/ASA Journal on Uncertainty Quantification*, 9(4):1673-1713.

McHutchon, A. and Rasmussen, C. E. (2011). Gaussian process training with input noise. In *Advances in Neural Information Processing Systems*, pages 1341-1349.

Mohammadi, H., Riche, R. L., Durrande, N., Touboul, E., and Bay, X. (2016). An analytic comparison of regularization methods for Gaussian processes. *arXiv preprint arXiv:1602.00853*.

Pearce, M. A. L., Poloczek, M., and Branke, J. (2022). Bayesian optimization allowing for common random numbers. *Operations Research*, 70(6):3457-3472.

Picheny, V., Ginsbourger, D., and Richet, Y. (2010). Noisy expected improvement and on-line computation time allocation for the optimization of simulators with tunable fidelity. In *2nd International Conference on Engineering Optimization, September 6-9, 2010, Lisbon, Portugal*.

Picheny, V., Ginsbourger, D., Richet, Y., and Caplin, G. (2013a). Quantile-based optimization of noisy computer experiments with tunable precision. *Technometrics*, 55(1):2-13.

Picheny, V., Moss, H., Torossian, L., and Durrande, N. (2022). Bayesian quantile and expectile optimisation. In *Uncertainty in Artificial Intelligence*, pages 1623-1633. PMLR.

Picheny, V., Wagner, T., and Ginsbourger, D. (2013b). A benchmark of kriging-based infill criteria for noisy optimization. *Structural and multidisciplinary optimization*, 48:607-626.

Pinder, T. and Dodd, D. (2022). GPJax: A Gaussian process framework in JAX. *Journal of Open Source Software*, 7(75):4455.

Plumlee, M. and Tuo, R. (2014). Building accurate emulators for stochastic simulations via quantile kriging. *Technometrics*, 56(4):466-473.

Qing, J., Dhaene, T., and Couckuyt, I. (2022). Spectral representation of robustness measures for optimization under input uncertainty. In *38th International Conference on Machine Learning (ICML)*.

Quadrianto, N., Kersting, K., Reid, M. D., Caetano, T. S., and Buntine, W. L. (2009). Kernel conditional quantile estimation via reduction revisited. In *Data Mining, 2009. ICDM'09. Ninth IEEE International Conference on*, pages 938-943. IEEE.

Quan, N., Yin, J., Ng, S. H., and Lee, L. H. (2013). Simulation optimization via kriging: a sequential search using expected improvement with computing budget constraints. *Iie Transactions*, 45(7):763-780.
R Core Team (2024). *R: A Language and Environment for Statistical Computing*. R Foundation for Statistical Computing, Vienna, Austria.

Rasmussen, C. E. and Williams, C. (2006). *Gaussian Processes for Machine Learning*. MIT Press.

Rivier, M. and Congedo, P. M. (2022). Surrogate-assisted bounding-box approach applied to constrained multi-objective optimisation under uncertainty. *Reliability Engineering & System Safety*, 217:108039.

Rojas-Gonzalez, S. and Van Nieuwenhuyse, I. (2020). A survey on kriging-based infill algorithms for multiobjective simulation optimization. *Computers & Operations Research*, 116:104869.

Roustant, O., Ginsbourger, D., and Deville, Y. (2012). DiceKriging, DiceOptim: Two R packages for the analysis of computer experiments by kriging-based metamodeling and optimization. *Journal of Statistical Software*, 51(1):1–55.

Roustant, O., Padonou, E., Deville, Y., Clément, A., Perrin, G., Giorla, J., and Wynn, H. (2020). Group kernels for Gaussian process metamodels with categorical inputs. *SIAM/ASA Journal on Uncertainty Quantification*, 8(2):775–806.

Santner, T. J., Williams, B. J., and Notz, W. I. (2018). Some criterion-based experimental designs. In *The Design and Analysis of Computer Experiments*, pages 201–246. Springer.

Sauer, A., Cooper, A., and Gramacy, R. B. (2023a). Vecchia-approximated deep Gaussian processes for computer experiments. *Journal of Computational and Graphical Statistics*, 32(3):824–837.

Sauer, A., Gramacy, R. B., and Higdon, D. (2023b). Active learning for deep Gaussian process surrogates. *Technometrics*, 65(1):4–18.

Saul, A. D., Hensman, J., Vehtari, A., and Lawrence, N. D. (2016). Chained Gaussian processes. In *Artificial intelligence and statistics*, pages 1431–1440. PMLR.

Shah, A., Wilson, A., and Ghahramani, Z. (2014). Student-t processes as alternatives to Gaussian processes. In *Artificial intelligence and statistics*, pages 877–885. PMLR.

Snelson, E. and Ghahramani, Z. (2006). Variable noise and dimensionality reduction for sparse Gaussian processes. In *Proceedings of the Twenty-Second Conference on Uncertainty in Artificial Intelligence*, pages 461–468.

Snelson, E., Rasmussen, C. E., and Ghahramani, Z. (2004). Warped Gaussian processes. *Advances in neural information processing systems*, 16:337–344.

Titsias, M. K. and Lázaro-Gredilla, M. (2011). Variational heteroscedastic Gaussian process regression. In *Proceedings of the 28th International Conference on Machine Learning (ICML-11)*, pages 841–848.
van Beek, A., Ghumman, U. F., Munshi, J., Tao, S., Chien, T., Balasubramanian, G., Plumlee, M., Apley, D., and Chen, W. (2021). Scalable adaptive batch sampling in simulation-based design with heteroscedastic noise. *Journal of Mechanical Design*, 143(3).

van der Herten, J., Couckuyt, I., and Dhaene, T. (2016). Hypervolume-based multi-objective Bayesian optimization with student-t processes. *arXiv preprint arXiv:1612.00393*.

Vanhatalo, J., Jylänki, P., and Vehtari, A. (2009). Gaussian process regression with student-t likelihood. *Advances in neural information processing systems*, 22.

Vazquez, E., Villemonteix, J., Sidorkiewicz, M., and Walter, E. (2008). Global optimization based on noisy evaluations: An empirical study of two statistical approaches. *Journal of Physics: Conference Series*, 135:012100.

Wang, W. and Haaland, B. (2019). Controlling sources of inaccuracy in stochastic kriging. *Technometrics*.

Wang, Z., Shi, J. Q., and Lee, Y. (2017). Extended t-process regression models. *Journal of Statistical Planning and Inference*, 189:38-60.

Wilson, J. T., Borovitskiy, V., Terenin, A., Mostowsky, P., and Deisenroth, M. P. (2021). Pathwise conditioning of Gaussian processes. *Journal of Machine Learning Research*, 22(105):1-47.

Xie, G. and Chen, X. (2017). A heteroscedastic t-process simulation metamodeling approach and its application in inventory control and optimization. In *2017 Winter Simulation Conference (WSC)*, pages 3242-3253. IEEE.

Yazdi, F. (2022). *Fast deep Gaussian process modeling and design for large complex computer experiments*. PhD thesis, Simon Fraser University.

Zhang, B., Gramacy, R. B., Johnson, L. R., Rose, K. A., and Smith, E. (2022). Batch-sequential design and heteroskedastic surrogate modeling for delta smelt conservation. *The Annals of Applied Statistics*, 16(2):816-842.

Zhang, Y., Tao, S., Chen, W., and Apley, D. W. (2020). A latent variable approach to Gaussian process modeling with qualitative and quantitative factors. *Technometrics*, 62(3):291-302.

Zhou, Q., Qian, P. Z., and Zhou, S. (2011). A simple approach to emulation for computer models with qualitative and quantitative factors. *Technometrics*, 53(3).

Zhu, X. and Sudret, B. (2021). Emulation of stochastic simulators using generalized lambda models. *SIAM/ASA Journal on Uncertainty Quantification*, 9(4):1345-1380.