--- title: Simulation Functions in Computer Experiments url: https://www.emergentmind.com/topics/simulation-functions-sfs type: topic --- # Simulation Functions in Computer Experiments 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 [1601.05887]. ## 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 [1601.05887]. 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 $$ S(a)=\{x:y(x)=a\}, $$ 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 [1601.05887]. 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\) 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, $$ Y(x)=\mu(x)+Z(x), $$ where \(\mu(x)\) is a mean or regression function and \(Z(x)\) is a zero-mean GP with variance \(\sigma^2\). In the applications discussed, the mean is taken to be constant and the correlation is stationary power-exponential, $$ R(x,x')=\exp\left(-\sum_{j=1}^d \theta_j |x_j-x_j'|^{p_j}\right), $$ with \(\theta_j \ge 0\) governing sensitivity to the \(j\)-th input and \(1\le p_j\le 2\) governing smoothness in that dimension [1601.05887]. After \(n\) simulator runs, the predictive distribution at a new input is $$ Y(x)\mid \text{data} \sim N(\hat y(x), s^2(x)). $$ Under Gaussian assumptions, \(\hat y(x)\) and \(s^2(x)\) are the conditional mean and conditional variance. Even without strict normality, \(\hat y(x)\) is the best linear unbiased predictor and \(s^2(x)\) its associated mean squared prediction error. For deterministic simulators, the emulator interpolates: at any observed design point, \(s(x)=0\) [1601.05887]. 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 \(\hat y(x)\) and \(s(x)\); 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 \(z\) is transformed to \(y=\sqrt z\) after cross-validation comparisons with \(z\) and \(\log(z+1)\) [1601.05887]. ## 3. Expected improvement as the sequential design principle The core idea is to define an improvement function \(I(x)\) that quantifies the value of running the simulator at a candidate input \(x\). Because the output at \(x\) is unknown before evaluation, \(I(x)\) is random under the GP predictive distribution. The design criterion is its expectation, $$ E[I(x)] = \int I(x)\, f(y\mid x)\, dy, $$ and the next run is chosen at $$ x_{n+1}=\arg\max_{x\in\chi} E[I(x)]. $$ This is the general EI template [1601.05887]. For optimization, if \(y_n\) is the best observed value so far under minimization, $$ y_n=\min\{y(x^{(1)}),\dots,y(x^{(n)})\}, $$ the improvement is $$ I(x)=\max\{y_n-y(x),0\}. $$ Under the GP predictive model, $$ E[I(x)] = s(x)\phi(u) + (y_n-\hat y(x))\Phi(u), $$ with $$ u=\frac{y_n-\hat y(x)}{s(x)}. $$ The first term, \(s(x)\phi(u)\), is the exploration term; the second, \((y_n-\hat y(x))\Phi(u)\), 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 [1601.05887]. This balance is the main reason EI is used for simulation functions rather than simply optimizing \(\hat y(x)\). Near an already sampled optimum, \(s(x)\) 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 \(a\), the improvement function is defined by $$ I(x)=\epsilon^2(x)-\min\{(y(x)-a)^2,\epsilon^2(x)\}, $$ with $$ \epsilon(x)=\alpha s(x), $$ for some positive constant such as \(\alpha=1.96\). 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 $$ \begin{aligned} E[I(x)] &= [\epsilon^2(x)-(\hat y(x)-a)^2]\left(\Phi(u_2)-\Phi(u_1)\right) \\ &\quad + s^2(x)\left[(u_2\phi(u_2)-u_1\phi(u_1))-(\Phi(u_2)-\Phi(u_1))\right] \\ &\quad + 2(\hat y(x)-a)s(x)\left(\phi(u_2)-\phi(u_1)\right), \end{aligned} $$ where $$ u_1=\frac{a-\hat y(x)-\epsilon(x)}{s(x)}, \qquad u_2=\frac{a-\hat y(x)+\epsilon(x)}{s(x)}. $$ 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 [1601.05887]. The same construction extends to multiple contours by replacing the squared deviation from a single target with the minimum squared deviation from several targets: $$ I(x)=\epsilon^2(x)-\min\left\{(y(x)-a_1)^2,\ldots,(y(x)-a_k)^2,\epsilon^2(x)\right\}. $$ Percentile contour estimation uses $$ I^g(x)=\epsilon^g(x)-\min\{(y(x)-\hat\nu_p)^g,\epsilon^g(x)\}, $$ where \(\hat\nu_p\) is the current estimate of the \(100p\)-th percentile of the simulator output and therefore changes during the sequential procedure [1601.05887]. For constrained optimization, an expensive constraint function \(c(x)\) can be emulated separately, feasibility \(a