Simulation Functions in Computer Experiments

• Simulation Functions are expensive deterministic mappings from inputs to outputs that are modeled with Gaussian processes to economize computer experiment runs.
• The approach employs expected improvement criteria to balance exploration and exploitation for objectives like global optimization and contour estimation.
• Empirical examples such as tidal power and volcano simulations illustrate effective surrogate modeling, uncertainty handling, and adaptive sampling.

In the design of computer experiments, a simulation function is the expensive simulator response $y(x)$ generated by a computer code at an input $x \in \chi \subset \mathbb{R}^d$. The modern treatment of such functions assumes that the simulator is deterministic but insufficiently sampled, so the unknown response surface is emulated statistically—primarily with a Gaussian process (GP)—and new simulator evaluations are chosen sequentially by maximizing an objective-specific expected improvement (EI) criterion. Within this framework, simulation functions are not studied only for interpolation; they are interrogated for features such as global optima, contours $S(a)=\{x:y(x)=a\}$, target-matching input sets, feasibility regions, robust designs, and related constrained or multi-fidelity objectives (Bingham et al., 2016).

1. Formal setting and scientific objectives

A computer experiment begins with a simulator or computer code that maps an input vector $x$ to an output $y(x)$, possibly one scalar output of interest extracted from a larger multivariate response. The central methodological assumption is that the simulator is expensive to run, so direct dense exploration of $\chi$ is infeasible, especially in moderate or high dimension. The practical problem is therefore not merely to approximate $y(x)$, but to allocate a small run budget toward a specific scientific goal (Bingham et al., 2016).

The canonical objectives are global optimization, contour or level-set estimation, and target matching. Global optimization seeks $\min_{x\in\chi} y(x)$ or $\max_{x\in\chi} y(x)$, with maximization typically handled by minimizing $-y(x)$. Contour estimation seeks the input set

$x \in \chi \subset \mathbb{R}^d$0

which is also the natural formalization of inverse design or “moving the output to a target.” The same framework extends to simultaneous multiple contour estimation, percentile contours, constrained optimization, robust design, noisy-simulator optimization, and multi-fidelity settings (Bingham et al., 2016).

Although the simulator is deterministic, the methodology models the unknown response surface as random away from sampled inputs. This stochasticity is epistemic rather than physical: rerunning the code at the same $x \in \chi \subset \mathbb{R}^d$1 yields the same output, but uncertainty persists because only finitely many inputs have been evaluated. A plausible implication is that the statistical model is not a replacement for the simulator; it is a device for deciding where the next simulator run is most informative.

2. Gaussian-process emulation of the simulation function

The default surrogate is a Gaussian process emulator,

$x \in \chi \subset \mathbb{R}^d$2

where $x \in \chi \subset \mathbb{R}^d$3 is a mean or regression function and $x \in \chi \subset \mathbb{R}^d$4 is a zero-mean GP with variance $x \in \chi \subset \mathbb{R}^d$5. In the applications discussed, the mean is taken to be constant and the correlation is stationary power-exponential,

$x \in \chi \subset \mathbb{R}^d$6

with $x \in \chi \subset \mathbb{R}^d$7 governing sensitivity to the $x \in \chi \subset \mathbb{R}^d$8-th input and $x \in \chi \subset \mathbb{R}^d$9 governing smoothness in that dimension (Bingham et al., 2016).

After $S(a)=\{x:y(x)=a\}$0 simulator runs, the predictive distribution at a new input is

$S(a)=\{x:y(x)=a\}$1

Under Gaussian assumptions, $S(a)=\{x:y(x)=a\}$2 and $S(a)=\{x:y(x)=a\}$3 are the conditional mean and conditional variance. Even without strict normality, $S(a)=\{x:y(x)=a\}$4 is the best linear unbiased predictor and $S(a)=\{x:y(x)=a\}$5 its associated mean squared prediction error. For deterministic simulators, the emulator interpolates: at any observed design point, $S(a)=\{x:y(x)=a\}$6 (Bingham et al., 2016).

The GP contains unknown mean, correlation, and variance parameters, usually estimated by maximum likelihood or Bayesian methods. The paper emphasizes that plug-in estimation makes the predictive distribution only approximately normal and may understate uncertainty if correlation-parameter uncertainty is ignored. This matters operationally because EI is computed from $S(a)=\{x:y(x)=a\}$7 and $S(a)=\{x:y(x)=a\}$8; if either is poorly estimated, the sequential design can misallocate expensive evaluations. Transformations of the response can therefore be essential. In the volcano example, the output $S(a)=\{x:y(x)=a\}$9 is transformed to $x$0 after cross-validation comparisons with $x$1 and $x$2 (Bingham et al., 2016).

3. Expected improvement as the sequential design principle

The core idea is to define an improvement function $x$3 that quantifies the value of running the simulator at a candidate input $x$4. Because the output at $x$5 is unknown before evaluation, $x$6 is random under the GP predictive distribution. The design criterion is its expectation,

$x$7

and the next run is chosen at

$x$8

This is the general EI template (Bingham et al., 2016).

For optimization, if $x$9 is the best observed value so far under minimization,

$y(x)$0

the improvement is

$y(x)$1

Under the GP predictive model,

$y(x)$2

with

$y(x)$3

The first term, $y(x)$4, is the exploration term; the second, $y(x)$5, is the exploitation term. The method therefore balances uncertain regions against currently promising regions without collapsing to either pure uncertainty sampling or pure surrogate optimization (Bingham et al., 2016).

This balance is the main reason EI is used for simulation functions rather than simply optimizing $y(x)$6. Near an already sampled optimum, $y(x)$7 is small or zero, so expected improvement can vanish even if the predictive mean is favorable. Conversely, highly uncertain regions are not automatically selected unless they have a plausible chance of improving the current objective.

4. Objective-specific EI criteria

For contour estimation or target matching at level $y(x)$8, the improvement function is defined by

$y(x)$9

with

$\chi$0

for some positive constant such as $\chi$1. Positive improvement occurs only when the output lies within an uncertainty-scaled tolerance band around the target. Under the GP model, the expected improvement is

$\chi$2

where

$\chi$3

The first term mainly favors points near the predicted contour; the last term can favor more distant but highly uncertain points; the second term is typically less influential (Bingham et al., 2016).

The same construction extends to multiple contours by replacing the squared deviation from a single target with the minimum squared deviation from several targets:

$\chi$4

Percentile contour estimation uses

$\chi$5

where $\chi$6 is the current estimate of the $\chi$7-th percentile of the simulator output and therefore changes during the sequential procedure (Bingham et al., 2016).

For constrained optimization, an expensive constraint function $\chi$8 can be emulated separately, feasibility $\chi$9 is represented through its predictive probability, and optimization EI is multiplied by that feasibility probability. For noisy or stochastic simulators, the paper cites a quantile-based objective using

$y(x)$0

with improvement defined relative to the smallest lower-quantile estimate observed so far. Robust design replaces a single deterministic response target with objectives involving $y(x)$1 or $y(x)$2 for controllable variables $y(x)$3, and multi-fidelity settings use augmented EI to trade off fidelity against evaluation cost (Bingham et al., 2016).

A plausible implication is that “expected improvement” is not one fixed acquisition function but a family of objective-dependent utilities. What remains invariant is the logic: evaluate next where the predictive distribution implies the greatest expected gain relative to the feature of the simulation function that matters scientifically.

5. Sequential workflow and empirical illustrations

The practical workflow begins with an initial space-filling design, typically a maximin Latin hypercube. A heuristic mentioned in the paper is $y(x)$4 initial runs for input dimension $y(x)$5. The simulator is evaluated at these initial points, a GP is fit—possibly after transforming the output—and EI is computed over candidate inputs. The next run is placed at the maximizer of EI, the surrogate is updated, and the process is repeated until the maximum EI becomes “small enough,” although no universal stopping rule is imposed (Bingham et al., 2016).

The Bay of Fundy example illustrates optimization. The simulator returns average extractable tidal power as a function of turbine location $y(x)$6, with $y(x)$7 along the tidal flow and $y(x)$8 perpendicular to it. Starting from a 20-point maximin Latin hypercube, the best observed value is $y(x)$9 MW. The EI surface then selects $\min_{x\in\chi} y(x)$0, where the expected improvement is 13.9 MW, but the actual simulator value is 159.7 MW. The sequential design ultimately finds the global maximum on a $\min_{x\in\chi} y(x)$1 grid with far fewer than 533 full-grid runs. The example also demonstrates that emulator uncertainty can be underestimated near a sharp peak, so realized improvement may greatly exceed predicted improvement (Bingham et al., 2016).

The volcano example illustrates contour estimation. Using TITAN2D, the inputs are pyroclastic flow volume and basal friction angle, the output is maximum flow height at a critical location, and the transformed response is $\min_{x\in\chi} y(x)$2. With 32 initial simulator runs, the objective is to estimate the contour $\min_{x\in\chi} y(x)$3. The contour EI criterion selects the next point near $\min_{x\in\chi} y(x)$4, a location close to the predicted contour but with relatively large predictive standard error. This embodies the intended behavior of contour EI: refine the scientifically relevant level set where the emulator is informative but still uncertain (Bingham et al., 2016).

6. Assumptions, limitations, and terminological scope

The methodology rests on a small number of consequential assumptions. The simulator is treated primarily as deterministic; the GP uncertainty is therefore epistemic. The standard closed-form EI expressions require the predictive distribution $\min_{x\in\chi} y(x)$5, but parameter estimation may make this only approximate. Transformations may be needed when GP assumptions are poor, and underestimation of predictive variance can mislead EI, as the tidal-power example shows. The framework is especially valuable in higher dimensions and small-budget settings, but it does not eliminate the need for model checking or careful surrogate specification (Bingham et al., 2016).

The term “SF” also requires disambiguation. In the arXiv literature, the same acronym denotes optical and mid-infrared ensemble structure functions in AGN variability (Li et al., 2023), scalar shrinkage functions in transform-domain image denoising (Hel-Or et al., 2020), atom-centered symmetry functions and polynomial symmetry functions in HDNNPs (Bircher et al., 2020), and scoring functions in knowledge graph embedding (Zhang et al., 2019). Within the computer-experiment context, however, the relevant meaning is the simulation function $\min_{x\in\chi} y(x)$6: the expensive black-box mapping interrogated by GP emulation and sequential EI design.

This terminological overlap suggests an important editorial caution. “Simulation Functions (SFs)” in design-of-experiments literature refers neither to variability statistics, denoising nonlinearities, atomistic descriptors, nor KGE plausibility maps, even though all of those are established uses of the same acronym. In the present sense, the concept is inseparable from expensive deterministic simulation, surrogate modeling, and adaptive sampling under objective-specific expected improvement.