Proceedings of the 2017 Winter Simulation Conference
W. K. V. Chan, A. D’Ambrogio, G. Zacharewicz, N. Mustafee, G. Wainer, and E. Page, eds.

IMPROVING PREDICTION FROM STOCHASTIC SIMULATION VIA MODELDISCREPANCY LEARNING
Henry Lam
Matthew Plumlee
Xinyu Zhang
Department of Industrial and Operations Engineering
University of Michigan
1205 Beal Ave.
Ann Arbor, MI 48109, USA

ABSTRACT
Stochastic simulation is an indispensable tool in operations and management applications. However,simulation models are only approximations to reality, and typically bear discrepancies with the generatingprocesses of real output data. We investigate a framework to statistically learn these discrepancies under thepresence of data on past implemented system configurations, which allows us to improve prediction using
simulation models. We focus on the case of general continuous output data that generalizes previous work.
Our approach utilizes (a combination of) regression analysis and optimization formulations constrained onsuitable summary statistics. We demonstrate our approach with a numerical example.
1

INTRODUCTION

Stochastic simulation is used ubiquitously in decision analytics in many operations and management
applications. It describes system dynamics under alternate configurations that can be used for analyticaltasks. While modelers often attempt to build models as realistic as possible, due to resource constraints
and ignorance of unobserved system features, these models are at best approximations and in almost allpractical cases bear discrepancies with reality. For example, in operational settings that are naturally castas “queues”, the serial structures of interarrival times (e.g., Livny et al. 1993), discretionary or strategic
behaviors (e.g., balking, reneging, queue selection; Pazgal and Radas 2008, Veeraraghavan and Debo 2009),
individualized and contextual behaviors of the customers and servers (e.g., Aksin et al. 2007), and structural
server changes over time (e.g., Brown et al. 2005) are all difficult to capture. Inadequate reflection of thesefeatures can degrade the predictive power of the simulation model and potentially affect the reliability ofdecision-making.
In this paper, we investigate a general framework that builds on our previous work Plumlee andLam (2016) to learn about the imperfectness of stochastic simulation models and subsequently use thisinformation to improve prediction. To describe our investigation, we first introduce some notations anddefine the notion of model discrepancy. Imagine that we are interested in an output variable, Y ∈ Y ,from a stochastic system which has a true distribution π(·). The variable can represent waiting time,queue length, reward, other responses or any collections of them that are relevant to decision-making. A
typical performance analysis consists of evaluating a quantity ρ(π), which in this paper we focus on the
expectation-type performance measure Eπ [ f (Y )] for some function f (·) : Y → R.
When real-world observations from π are abundant, the performance analyses can be accurately
approximated by replacing π with data, e.g., the empirical distribution. However, in decision analyses (e.g.,
optimality or feasibility tests), we are typically interested in system designs that are sparsely sampled oreven never adopted before. So, instead, model builders use “simulation models” that generate outputs with

Lam, Plumlee, and Zhang
distribution, say, π̃(·), as an approximation to π(·) that can be computed much more easily via simulation
replications. If π̃(·) 6= π(·), then ρ(π̃) may not equal ρ(π) and the outcomes of the analysis become
erroneous. Our interest is to combine both the observations on π and the simulation model π̃ to generateprediction for ρ(π) that is better than using either one only.
Our approach hinges on estimating the difference between π̃(·) and π(·) and then integrating this bias
estimate with the simulation model to predict ρ(π). It requires the existence of some control variable ordesign point x ∈ X that represents the set of system parameters, and that real system observations areavailable at some x. The value of a design point can include decision variables (e.g., number of servers,
routing policy) and important components of the system (e.g., arrival rate of customers). The availability ofobservations from different design points provides an opportunity to learn and correct for the imperfectnessof π̃, including at un- or under-observed points, by data-pooling.
We consider the class of mappings P : X × Y → R, such that if p ∈ P, then p(x, ·) is a valid
probability measure over Y . We expand the definition of the true model π and the simulation model π̃ tobe in P. In general, we have the pointwise relationship
π(x, y) = π̃(x, y) + δ (x, y)(1)

where we call δ (·, ·) the model discrepancy. When π̃ 6= π, we have δ 6= 0 and model discrepancy is present.
The relationship (1) can be viewed as a stochastic counterpart of a similar discrepancy notion in theliterature of deterministic computer experiments (where no y in (1) is present). The latter is under theframework of model calibration that refers generally to the refinement or correction of a hypotheticalsimulation model through real-world data to enhance the model prediction capability. Starting with theseminal statistical article of Kennedy and O’Hagan (2001), calibration for deterministic simulation has
been substantially investigated (e.g., Bayarri et al. 2012, Oakley and O’Hagan 2004, Higdon et al. 2004,
Plumlee 2016). In the stochastic simulation literature, other than parameter calibration (implemented
in some simulation software such as Anylogic), calibration on model structure is commonly conducted
together with model validation, through iterative checking (e.g., by statistical tests) and re-building of thesimulation model until a satisfactory model is obtained (Balci and Sargent 1982, Sargent 2005, Sargent 2013,
Kleijnen 1995). This paper detours from these common suggestions, and utilizes the idea of inferring model
discrepancies in deterministic computer experiments to improve the prediction of stochastic simulation.
Compared to the established practice, the avoidance of building increasingly sophisticated models could
save costs and time in refining the model, and prevent the unfortunate situation in case an ultimately
satisfactory model is absent.
We will consider two methods in this paper. The first one estimates ρ(δ (x, ·)) := ρ(π(x, ·)) − ρ(π̃(x, ·))directly by regression against the covariate x. The second method exploits further the fact that π and π̃
are probability distributions, so that they satisfy standard conditions for valid probability measures that
constrain the possible values of ρ(π). This second method can be thought of as inferring the difference ofdistributions δ (x, ·) = π(x, ·) − π̃(x, ·) and using it to calculate ρ(π). However, since δ (x, ·) can be a highor even infinite-dimensional object (in the case of continuous output distributions) that is challenging toestimate, our second approach instead operates on a collection of summary statistics on δ (x, ·), and posits
suitable optimization formulations to compute bounds that take into account the relation between thesesummary statistics induced by the probability distribution structure.
As mentioned before, this work builds on Plumlee and Lam (2016) and generalizes their investigation.
Besides the difference regarding the representation of model discrepancy (additive versus multiplicativeform), Plumlee and Lam (2016) focuses on finite-valued outputs in a Bayesian framework. In this paper, weare primarily interested in outputs that are continuous, which pose substantial additional challenges since they
are now infinite-dimensional objects. Our aforementioned second method novelly uses summary statistics
to facilitate tractable statistical estimation, but with infinite-dimensional optimization programs positedover the space of probability distributions to recover bounds incurred by the continuous distributions. With
well-chosen functions to define the summary statistics, these optimization programs can be reformulated by

Lam, Plumlee, and Zhang
duality as tractable finite-dimensional programs. These developments provide an implementable mechanism
to take into account distributional information that is infinite-dimensional in nature.
In the rest of this paper, we first describe the two proposed methods in Sections 2 and 3. Then we
provide an extensive numerical study and compare our performance with both simulation-only and data-only
predictions in Section 4. We conclude this paper in Section 5.
2

LEARNING MODEL DISCREPANCY VIA REGRESSION

Consider a collection of design point values, say {1, . . . , s}. For each value x = i, suppose there are i.i.d.
real-world output observations Yi j ∼ π(i, ·), j = 1, . . . , ni . The sample size ni can be possibly zero for somei (so that no Yi j exists).
We assume there is enough computational resource, so that the simulation output distribution π̃ canbe exactly obtained.
Suppose we are interested in ρ(π(x0 , ·)) = Eπ(x0 ,·) [ f (Y )] for some x0 . Using data alone, we will output
nx0
f (Yi j ). Clearly, if there is no data at x0 , this approach does not apply. On the otherȲx0 · = (1/nx0 ) ∑ j=1
hand, the simulation model gives ρ(π̃(x0 , ·)) = Eπ̃(x0 ,·) [ f (Y )] as the predicted value of ρ(π(x0 , ·)).
To combine both {Yi j } and π̃ in our prediction, we consider the estimation of ρ(δ (x0 , ·)) by the
regression problemf (Yi j ) − ρ(π̃(i, ·)) = h(i) + εi j(2)where h(·) is an unknown function of the design point x, and εi j , j = 1, . . . , ni are mean-zero i.i.d. noises
from the real data, for each i. The responses in the regression are f (Yi j ) − ρ(π̃(i, ·)) whose mean values
are the model discrepancies ρ(δ (i, ·)).
Regression (2) can be conducted by a variety of tools; we use local linear regression with a standard
R package in our experiments. This gives us a confidence interval (CI) for ρ(δ (x0 , ·)), say [L,U]. Our
prediction of ρ(π(x0 , ·)) is then expressed as a CI [ρ(π̃(x0 , ·)) + L, ρ(π̃(x0 , ·)) +U]. Clearly, if [L,U] indeed
covers the true ρ(δ (x0 , ·)) with probability 1 − α, then [ρ(π̃(x0 , ·)) + L, ρ(π̃(x0 , ·)) +U] covers ρ(π(x0 , ·))with the same probability.
The above scheme can be modified withf (Yi j ) = β ρ(π̃(i, ·)) + h(i) + εi j

(3)

The additional unknown scaling parameter β represents a multiplicative model discrepancy. The responsein (3) now has mean ρ(π(i, ·)). The CI for ρ(π(x0 , ·)) can be obtained directly from the regression (3).
3

COMBINING REGRESSION OUTCOMES WITH OPTIMIZATION

We consider an enhancement of the method in Section 2 by introducing the auxiliary statistics γk (π(x, ·)) =Eπ(x,·) [gk (Y )] for k = 1, . . . , m. One of these statistics can be taken as ρ(π(x, ·)), e.g., let g1 = f . Then,instead of using (2) only, we impose the multivariate-output regression
gk (Yi j ) − γk (π(i, ·)) = hk (i) + εi jk(4)where hk (·) is the individual regression function for the k-th statistic. Similar to Section 2, εi jk , j = 1, . . . , ni
are mean-zero i.i.d. noises, for each i, and gk (Yi j ) − γk (π(i, ·)) has mean equal to γk (δ (i, ·)).
The idea is that estimates of γk (δ (x0 , ·)), k = 1, . . . , m give more information than using ρ(δ (x0 , ·))alone. For convenience, we denote γ = (γk )k=1,...,m and g = (gk )k=1,...,m . To incorporate the new information,
we first need a confidence region (CR) for γ(δ (x0 , ·)). As in Section 2, this can be done by a variety ofstatistical tools, and here we use a basic scheme of running individual local regression at each k, making an
additional assumption that εi jk , k = 1, . . . , m are independent for each i, j. For each individual regression,
we can obtain the point estimate of γk (δ (x0 , ·)), say γ̂k , and its standard error, say sk . Under standardLam, Plumlee, and Zhang
approx.

assumption we have γ̂k ∼ N(γk (δ (x0 , ·)), σk2 ), independent among the k’s, where σk2 is the sampling
variance. This gives us a CR
(
)
m
2
(z
−
γ̂
)
k
k
2
U = (z1 , . . . , zm ) ∈ Rm : ∑≤ χ1−α,m
(5)
2
s
k=1
k
2
for γ(δ (x0 , ·)), where χ1−α,mis the (1 − α)-quantile of a χ 2 -distribution with degree of freedom m.
We then impose the optimization problems

andmax ρ(q) subject to γ(q) − γ(π̃(x0 , ·)) ∈ U(6)min ρ(q) subject to γ(q) − γ(π̃(x0 , ·)) ∈ U(7)where the decision variable q is a probability distribution on the response space Y . The constraints in(6) and (7) incorporate information on the possible values of γ(q) with probability 1 − α. Given thisinformation, the outcomes of (6) and (7) give the best upper and lower bounds on ρ(π(x0 , ·)). If the CR U
contains the true γ(π̃(x0 , ·)) with probability 1 − α, then (6) and (7) will form a CI for ρ(π(x0 , ·)) with level
at least 1 − α. This idea originates from the literature of data-driven distributionally robust optimization(e.g., Delage and Ye 2010, Ben-Tal et al. 2013, Goh and Sim 2010, Wiesemann et al. 2014), where U isoften known as the uncertainty set or the ambiguity set.
Note that (6) and (7) are generalized moment problems, where the optimization contains an objectivefunction and constraint functions that are all moments of a random variable. By using conic duality, (6)can be reformulated as the dual problemq2minκ,ν∈R,λ ∈Rm κ + (γ(π̃(x0 , ·)) + γ̂)T λ + χ1−α,mνT(8)subject toκ + g(y) λ − f (y) ≥ 0 for all y ∈ Yk(s)T λ k2 ≤ ν
where γ̂ = (γ̂k )k=1,...,m and s = (sk )k=1,...,m , and similarly for (7) (by considering − max{−ρ(π)}). Strong
duality holds if γ(q) − γ(π̃(x0 , ·)) lies in the interior of U and under appropriate topological assumptions on
Y (e.g., Shapiro 2001). Without the condition, weak duality still implies that (8) provides a conservative
approximation for (6).
Note that (8) still has an infinite number of constraints, but for specific choices of f and g one canreduce (8) to finite-dimensional optimization problems. We focus here on the setting that Y = [0, UB] ⊂ R
where UB is the largest possible value of Y , and f (y) and gk (y) are in the form yr where r is a rational
number. This allows us to reduce (8) to finite-dimensional semidefinite programs (SDP) by generalizing thetechnique in Bertsimas and Popescu (2005), which consider moment equalities only and does not containthe second order cone constraint in the dual formulation (8). For example, in the case that f (y) = y andgk (y) = yk/m̃ for k = 1, 2, · · · , m, (8) becomesq2minκ,ν∈R,λ ∈Rm κ + (γ(π̃(x0 , ·)) + γ̂)T λ + χ1−α,mνmsubject toκ + ∑ λi yi/m̃ − y ≥ 0 for all y ∈ [0, UB]
i=1
k(s)T λ k2 ≤ ν(9)Lam, Plumlee, and Zhang
By a change of variable y1/m̃ → y, we can transform the algebraic constraints into positive semidefinite
constraints, and getq2νminκ,ν∈R,λ ∈Rm ,U∈R(m+1)×(m+1) κ + (γ(π̃(x0 , ·)) + γ̂)T λ + χ1−α,msubject toU 0∑

ui j = 0, l = 1, · · · , mi, j:i+ j=2l+1

(10)

l

∑

 r/m̃
ui j = ∑ yr m̃−r
, l = 0, · · · , m
l−r UBi, j:i+ j=2l+2
k(s)T λ k2 ≤ ν

r=0

where y0 = κ, yi = λi 1i6=m̃ + (λi − 1)1i=m̃ , [U]i j = ui j . Similarly, for (7), the reformulation remains as (10)except that ym̃ is replaced by λm̃ + 1.
We choose the power in gk (·) as k/m̃ instead of integers so as to avoid a blow-up in the magnitudeof the moment as k increases. Note that there are other plausible choices of statistics, e.g., quantile-typestatistics in the form gk (y) = I(y ≤ bk ) for some bk , or a combination of these and the power moments,
which can also be similarly converted into SDP.
As in Section 2, an alternative is to consider the regression
gk (Yi j ) = βk γk (π̃(i, ·)) + hk (i) + εi jk ,(11)with responses gk (Yi j ) and additional scaling parameters βk . We can form a CR U for γ(π(x0 , ·)) in this caseas (5), where γ̂k is now the point estimate of γk (π(x0 , ·)) under the regression and sk is the corresponding
standard error. Then we use the optimization problems
max ρ(q) subject to γ(q) ∈ U
andmin ρ(q) subject to γ(q) ∈ U

The maximization can be dualized to (8) and in the special case considered above to (10) (without theγ(π̃(x0 , ·)) term). Similar treatments for the minimization as discussed above also hold.
We note that the accuracy of our proposed approach, here and in Section 2, depends on the complexityof the simulation model and the reality, in that it determines how challenging it is to find an acceptable
regression model for capturing the model discrepancy. In the next section, we will test our approach withsome numerical studies on a simple queueing system.
4

NUMERICAL STUDY

We use a queueing example to test the approaches we propose in Sections 2 and 3. We first consider asimulation model and the “real” system generated as follows: The simulation model is a M/M/x queuewith arrival rate λ = 5 and service rate µ = 0.05. The real system is a M/M/x queue with a random arrival
rate as an absolute value of a normally distributed random variable with mean 5 and standard deviation
5 · 0.1, and a random service rate as an absolute value of a normally distributed random variable withmean 0.045 and standard deviation 0.045 · 0.1 (Think of these as, e.g., some random daily characteristics).
In addition, each customer has a probability 0.2 to abandon the queue if the waiting time is larger thanan exponential random variable with mean 450
x where x is the number of servers. The output quantity ofinterest Y is the average waiting time among the first 10 customers.
The design point x corresponds to the number of servers, which ranges on x = 10, 11, · · · , 25. For thisexample, we can easily run plenty of simulation, say 105 replications, to get π̃ that can be considered asexactly obtained. For each x = 10, . . . , 24, we generate 10 observations on the real system. We deliberately
collect no observations for x = 25 as a test point.

Lam, Plumlee, and Zhang
We use f (y) = y and gk (y) = yk/m̃ , so that we can use the SDP (10) when adopting the optimizationenhanced approaches. We use a confidence level 1 − α = 95% in the experiment.
We consider five methods:•

Data-only method (Method D) only utilizes real data Yi j . For i = 1, . . . , 24, the CI is


σ̂i
σ̂i
,
Ȳi − z1−α/2 √ , Ȳi + z1−α/2 √
10
10
s
10

∑ Yi j
j=1

•
•
•
•
•

10

∑ (Yi j −Ȳi )2thwhere Ȳi = 10 , σ̂i =10−1 , and z1−α/2 is the (1 − α/2) quantile of a standard normal
distribution. For x = 25, there is no data so Method D does not apply.
Regression-only method 1 (Method R1) is based on the regression problem (2) in Section 2.
Optimization-enhanced method 1 (Method O1) is based on optimizing over the CR obtained from
(4) in Section 3.
Regression-only method 2 (Method R2) is based on (3) in Section 2.
Optimization-enhanced method 2 (Method O2) is based on optimizing over the CR obtained from(11) in Section 3.
Simulation-only method (Method S) uses the (accurate) point estimate of our simulation model.
i=1

Each method except Method S outputs a lower bound and a upper bound value for each experiment,
thereby forming a (1 − α)-level confidence interval. For Methods R1, O1, R2, and O2, we use the Rfunction “npreg” in the “np” package to make local linear regressions.
We first investigate the coverage probabilities of the CI generated by the methods. Table 1 shows theestimated coverage probabilities from repetition of 100 experiments with gk (y) = yk/10 , k = 1, · · · , 10. For
each entry in the table, we use a ± b to denote the CI of the estimated coverage probability from the 100
experiments, i.e., a is the average number of success in covering the true value and b is the half-width of
the associated CI for the coverage probability.
For Method D, the coverage probability ranges from 0.816 to 0.999. It is roughly around the theoretical
coverage probability 95% but there are mild fluctuations due to the small number of repetitions. That being
said, both regression and optimization methods obtain more fluctuated values as the number of servers
varies. For example, the coverage probability ranges from 0.372 to 0.999 for Method R1 and ranges from0.493 to 1 for Method O1. Therefore, the data-only method shows a more reliable performance overall.
Even though in some cases (e.g., when x = 20) the correction methods (i.e., Methods R1, O1, R2, O2) all
obtain reasonable coverages, the data-only method has uniformly better performances and is the best inother cases (e.g., x = 10). However, if one considers the case x = 25 (i.e., the extreme case where there isno data), then Method D does not apply, but the correction methods still give reasonable coverage (though
the performances vary). That the correction methods tend to have poorer coverage probabilities could be
due to the crudeness of the regression models, whose assumptions and post-hoc analyses have not been
carefully conducted in this experiment. In other words, the regression models used in this example could
deviate from the true behavior of ρ(δ (x, ·)) or ρ(π(x, ·)).

Lam, Plumlee, and Zhang
Table 1: Coverage probabilities for different methods
Numberservers
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

of

Method D

Method R1

Method O1

Method R2

Method O2

0.88±0.064
0.91±0.056
0.96±0.039
0.88±0.064
0.89±0.062
0.91±0.056
0.93±0.050
0.88±0.064
0.93±0.050
0.92±0.053
0.92±0.053
0.95±0.043
0.95±0.043
0.91±0.056
0.95±0.043-

0.47± 0.098
0.59± 0.097
0.67± 0.093
0.77± 0.083
0.76± 0.084
0.79± 0.080
0.85± 0.070
0.91± 0.056
0.93± 0.050
0.95± 0.043
0.96± 0.039
0.96± 0.039
0.90± 0.059
0.89± 0.062
0.89± 0.062
0.81± 0.077

0.59±0.097
0.63±0.095
0.78±0.082
0.83±0.074
0.79±0.080
0.78±0.082
0.91±0.056
0.95±0.043
0.91±0.056
0.97±0.034
0.96±0.039
0.97±0.034
0.89±0.062
0.88±0.064
0.92±0.053
0.70±0.090

0.57±0.098
0.67±0.093
0.81±0.077
0.90±0.059
0.90±0.059
0.98±0.028
0.94±0.047
0.98±0.028
0.95±0.043
0.90±0.059
0.99±0.020
0.95±0.043
0.99±0.020
0.99±0.020
0.87±0.066
0.87±0.066

0.66±0.093
0.77±0.083
0.88±0.064
0.93±0.050
0.90±0.059
0.96±0.039
0.96±0.039
0.97±0.034
0.93±0.050
0.94±0.047
0.97±0.034
0.94±0.047
0.95±0.043
0.95±0.043
0.87±0.066
0.73±0.087

Besides coverage probability, we are interested in gaining some understanding on how conservative
the methods are. This can be measured by the difference between the lower and upper confidence bounds
for each method. We call these differences the prediction gaps. Table 2 shows the statistics of this gap foreach of the methods, among the 100 experiments with gk (y) = yk/10 , k = 1, · · · , 10. For each entry in thetable except the ones in the last column, we use a ± b to denote the CI of the prediction gap from the 100
experiments (i.e., a is the average gap, and b is the CI half-width). The last column shows the absolutedifference between the point estimate of the simulation model and the truth, which can be viewed as themagnitude of the model error made in using the simulation.
Table 2: Prediction gaps for different methods
Numberservers
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

of

Method D

Method R1

Method O1

Method R2

Method O2

Method S

7.488 ± 0.355
6.005 ± 0.346
5.527 ± 0.275
4.751 ± 0.249
4.490 ± 0.220
3.878 ± 0.208
3.582 ± 0.177
2.987 ± 0.144
2.791 ± 0.138
2.533 ± 0.125
2.274 ± 0.108
2.112 ± 0.098
1.831 ± 0.086
1.673 ± 0.088
1.501 ± 0.077
1.422 ± 0.067

1.820 ± 0.263
1.561 ± 0.190
1.325 ± 0.150
1.315 ± 0.140
1.334 ± 0.134
1.172 ± 0.087
1.143 ± 0.087
1.038 ± 0.053
1.000 ± 0.057
0.981 ± 0.046
0.958 ± 0.043
0.906 ± 0.026
0.893 ± 0.025
0.885 ± 0.024
0.898 ± 0.019
0.999 ± 0.084

2.035 ± 0.242
1.670 ± 0.155
1.486 ± 0.127
1.467 ± 0.133
1.371 ± 0.103
1.247 ± 0.083
1.184 ± 0.086
1.060 ± 0.049
0.981 ± 0.050
0.973 ± 0.047
0.945 ± 0.045
0.886 ± 0.029
0.849 ± 0.028
0.832 ± 0.030
0.806 ± 0.026
0.871 ± 0.078

2.660 ± 0.205
2.631 ± 0.176
2.322 ± 0.125
2.127 ± 0.085
1.979 ± 0.040
1.890 ± 0.031
1.797 ± 0.030
1.713 ± 0.037
1.627 ± 0.051
1.527 ± 0.065
1.579 ± 0.067
1.489 ± 0.075
1.452 ± 0.081
1.522 ± 0.082
1.462 ± 0.084
1.514 ± 0.080

2.943 ± 0.171
2.884 ± 0.161
2.489 ± 0.094
2.360 ± 0.077
2.188 ± 0.050
2.035 ± 0.041
1.933 ± 0.030
1.808 ± 0.031
1.679 ± 0.038
1.557 ± 0.046
1.470 ± 0.047
1.383 ± 0.054
1.306 ± 0.045
1.272 ± 0.048
1.165 ± 0.050
1.157 ± 0.067

0.413
0.348
0.236
0.191
0.137
0.041
0.027
0.018
0.059
0.093
0.101
0.140
0.136
0.143
0.164
0.162

Lam, Plumlee, and Zhang
Table 2 shows that the prediction gaps using correction methods are significantly shorter than thatusing data-only method for every x. For example, at x = 18, the prediction gap ranges from 2.653 to 2.929
for Method D, but ranges only from 0.943 to 1.057 for Method R1 and from 0.937 to 1.037 for Method
O1. Moreover, at some particular x, both high coverage probability and small prediction gap are obtained
simultaneously for correction methods. For instance, the coverage probabilities at x = 18 for both data-only
and correction methods are roughly around 95% and the probabilities only differ by 5% ∼ 10% for many
choice of x such as 15, 16, 23, 24. This shows that, while the coverage probabilities are more reliable
for using only data, their predictions are generally much more conservative. The correction methods use
information on the simulation model to improve the conservativeness.
Considering the last column in Table 2, we see that the simulation model errors are in overall smaller
than the prediction gaps in this example. However, simulation method in this example is, in a sense, always
erroneous as the Monte Carlo errors of an inaccurate model are washed away with abundant simulation
runs, leaving point estimates that are different from the truth. We note that different simulation model will
give different model error, and there is no direct guarantee whether the model error as we define will besmaller or larger than the prediction gaps from the correction methods.
Next we repeat the experiment, but using a smaller number, but higher orders, of moments gk (y) =yk/2 , k = 1, · · · , 3. Tables 3 and 4 give the coverage probabilities and prediction gaps respectively under thisnew scheme. Table 3 shows that coverage probabilities for regression methods and optimization methods
are comparable to each other for every x. For example, for x = 16, the coverage probabilities differ by only
roughly 0.02 between regression methods and optimization methods. Under the new scheme, the prediction
gaps obtained from optimization methods in overall are shorter than that from regression methods. For
instance, in Table 4, we have prediction gap ranging from 1.106 to 1.34 for Method R1 and from 1.047 to
1.237 for Method O1 at x = 14. Also, the prediction gap ranges from 1.609 to 1.731 for Method R2 andranges from 1.417 to 1.539 for Method O2 at x = 19. Besides the shorter gaps, the optimization methods
give prediction gaps with smaller variability, e.g., the CI length of the prediction gap is shorter for almost
every x for Method O2 than that for Method R2.
Table 3: Coverage probabilities for different methods under a different set of moments
Numberservers
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

of

Method D

Method R1

Method O1

Method R2

Method O2

0.89 ± 0.06
0.93 ± 0.05
0.87 ± 0.07
0.91 ± 0.06
0.97 ± 0.03
0.90 ± 0.06
0.92 ± 0.05
0.92 ± 0.05
0.94 ± 0.05
0.90 ± 0.06
0.84 ± 0.07
0.93 ± 0.05
0.89 ± 0.06
0.87 ± 0.07
0.94 ± 0.05–

0.39 ± 0.10
0.59 ± 0.10
0.70 ± 0.09
0.79 ± 0.08
0.82 ± 0.08
0.82 ± 0.08
0.88 ± 0.06
0.89 ± 0.06
0.99 ± 0.02
0.98 ± 0.03
0.93 ± 0.05
0.97 ± 0.03
0.93 ± 0.05
0.93 ± 0.05
0.83 ± 0.07
0.74 ± 0.09

0.41 ± 0.10
0.59 ± 0.10
0.66 ± 0.09
0.75 ± 0.09
0.79 ± 0.08
0.79 ± 0.08
0.86 ± 0.07
0.85 ± 0.07
0.99 ± 0.02
0.98 ± 0.03
0.94 ± 0.05
0.97 ± 0.03
0.93 ± 0.05
0.93 ± 0.05
0.87 ± 0.07
0.68 ± 0.01

0.54 ± 0.10
0.76 ± 0.08
0.87 ± 0.07
0.92 ± 0.05
0.94 ± 0.05
0.95 ± 0.04
0.97 ± 0.03
0.95 ± 0.04
0.95 ± 0.04
0.99 ± 0.02
0.99 ± 0.02
0.95 ± 0.04
0.98 ± 0.03
0.98 ± 0.03
0.95 ± 0.04
0.77 ± 0.08

0.47 ± 0.10
0.71 ± 0.09
0.86 ± 0.07
0.90 ± 0.06
0.86 ± 0.07
0.93 ± 0.05
0.95 ± 0.04
0.94 ± 0.05
0.92 ± 0.05
0.97 ± 0.03
0.98 ± 0.03
0.94 ± 0.05
0.97 ± 0.03
0.97 ± 0.03
0.95 ± 0.04
0.78 ± 0.01

Lam, Plumlee, and Zhang
Table 4: Prediction gaps for different methods under a different set of moments
Numberservers
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

of

Method D

Method R1

Method O1

Method R2

Method O2

Method S

7.165 ± 0.390
6.156 ± 0.293
5.544 ± 0.285
4.850 ± 0.260
4.468 ± 0.208
3.729 ± 0.162
3.637 ± 0.176
3.145 ± 0.184
2.920 ± 0.138
2.563 ± 0.126
2.244 ± 0.120
2.053 ± 0.102
1.879 ± 0.089
1.690 ± 0.077
1.565 ± 0.072
1.388 ± 0.073

1.753 ± 0.238
1.498 ± 0.187
1.358 ± 0.157
1.295 ± 0.122
1.223 ± 0.117
1.229 ± 0.097
1.104 ± 0.082
1.082 ± 0.079
1.048 ± 0.066
0.977 ± 0.044
0.962 ± 0.039
0.917 ± 0.031
0.889 ± 0.025
0.867 ± 0.019
0.909 ± 0.026
1.017 ± 0.066

1.581 ± 0.190
1.345 ± 0.132
1.288 ± 0.129
1.202 ± 0.098
1.142 ± 0.095
1.150 ± 0.077
1.041 ± 0.064
1.018 ± 0.060
0.999 ± 0.050
0.941 ± 0.035
0.922 ± 0.029
0.904 ± 0.026
0.875 ± 0.024
0.873 ± 0.025
0.897 ± 0.028
0.966 ± 0.005

2.583 ± 0.176
2.314 ± 0.150
2.304 ± 0.115
2.146 ± 0.084
2.048 ± 0.056
1.885 ± 0.026
1.803 ± 0.030
1.754 ± 0.031
1.659 ± 0.047
1.670 ± 0.061
1.523 ± 0.065
1.557 ± 0.076
1.526 ± 0.083
1.523 ± 0.082
1.532 ± 0.078
1.519 ± 0.080

2.377 ± 0.149
2.079 ± 0.124
2.074 ± 0.097
1.950 ± 0.068
1.812 ± 0.041
1.665 ± 0.026
1.583 ± 0.034
1.585 ± 0.039
1.492 ± 0.051
1.478 ± 0.061
1.362 ± 0.064
1.383 ± 0.071
1.373 ± 0.075
1.376 ± 0.073
1.378 ± 0.071
1.345 ± 0.008

0.413
0.348
0.236
0.191
0.137
0.041
0.027
0.018
0.059
0.093
0.101
0.140
0.136
0.143
0.164
0.162

In comparing different regression models (2) and (3), we see that correction methods based on (3)generally give higher coverage probabilities compared with those based on (2). For example, at x = 12,
Method R1 obtains coverage probability 0.67 ± 0.093 and Method R2 obtains 0.81 ± 0.077 in Table 1,
and 0.70 ± 0.09 and 0.87 ± 0.07 respectively in Table 3. Likewise, for optimization methods, Method O1
obtains coverage probability 0.78 ± 0.082 and Method O2 obtains 0.77 ± 0.083 in Table 1, and 0.66 ± 0.09
and 0.86 ± 0.07 respectively in Table 3. However, higher coverage probabilities seem to come with the priceof wider prediction gaps. For instance, Method R1 gives prediction gaps 1.325 ± 0.015 and Method R2
gives 2.322 ± 0.0125 in Table 1, and 1.358 ± 0.157 and 2.304 ± 0.115 respectively in Table 3. Likewise, for
optimization methods, Method O1 gives prediction gaps 1.486±0.127 and Method O2 obtains 2.489±0.094
in Table 1, and 1.288 ± 0.129 and 2.074 ± 0.097 respectively in Table 3. From these observations, there
does not seem to be a concrete conclusion which regression model is better. One who prefers a higher
coverage probability at the cost of a conservative prediction gap may plausibly choose the second regression
model, and vice versa.
5

CONCLUSION

We investigate an approach based on a combination of regression and optimization to learn and correct for themodel discrepancy of a simulation model with real-world observations. Our approach targets at continuous
output variables that are not directly handleable by previous work due to their infinite-dimensional nature.
One version of our approach regresses the objective function of interest on the observed model discrepancies,or the observations themselves, against the design points. Another version further regresses the moments ofthe model discrepancies or observations and uses a moment optimization to compute bounds that accountfor the distributional nature of the outputs. We present the regression details and how to reformulate
and subsequently solve the optimization in the second version via duality and semidefinite programming.
Our experimental results show that the data-only approach could give better coverage probability in theavailability of data, but could be more conservative in terms of wider confidence interval length thanour correction methods. Moreover, the data-only approach is only applicable for historically observed
design points but trivially breaks down in the new design points. Our results also show that our correction
methods improve pure simulation models that are misspecified, in that the simulation-only estimates areLam, Plumlee, and Zhang
systematically different from the truth while our correction methods generate interval estimates that coveror are closer to it.
ACKNOWLEDGMENTS
We gratefully acknowledge support from the National Science Foundation under grants CMMI-1542020,
CMMI-1523453 and CAREER CMMI-1653339.
REFERENCES
Aksin, Z., M. Armony, and V. Mehrotra. 2007. “The modern call center: A multi-disciplinary perspectiveon operations management research”. Production and Operations Management 16 (6): 665–688.
Balci, O., and R. G. Sargent. 1982. “Some examples of simulation model validation using hypothesis
testing”. In Proceedings of the 14th conference on Winter Simulation-Volume 2, 621–629. Winter
Simulation Conference.
Bayarri, M. J., J. O. Berger, R. Paulo, J. Sacks, J. A. Cafeo, J. Cavendish, C.-H. Lin, and J. Tu. 2012. “Aframework for validation of computer models”. Technometrics.
Ben-Tal, A., D. Den Hertog, A. De Waegenaere, B. Melenberg, and G. Rennen. 2013. “Robust solutionsof optimization problems affected by uncertain probabilities”. Management Science 59 (2): 341–357.
Bertsimas, D., and I. Popescu. 2005. “Optimal inequalities in probability theory: A convex optimization
approach”. SIAM Journal on Optimization 15 (3): 780–804.
Brown, L., N. Gans, A. Mandelbaum, A. Sakov, H. Shen, S. Zeltyn, and L. Zhao. 2005. “Statistical analysis
of a telephone call center: A queueing-science perspective”. Journal of the American statisticalassociation 100 (469): 36–50.
Delage, E., and Y. Ye. 2010. “Distributionally robust optimization under moment uncertainty with application
to data-driven problems”. Operations research 58 (3): 595–612.
Goh, J., and M. Sim. 2010. “Distributionally robust optimization and its tractable approximations”. Operations research 58 (4-part-1): 902–917.
Higdon, D., M. Kennedy, J. C. Cavendish, J. A. Cafeo, and R. D. Ryne. 2004. “Combining field data andcomputer simulations for calibration and prediction”. SIAM Journal on Scientific Computing 26 (2):448–466.
Kennedy, M. C., and A. O’Hagan. 2001. “Bayesian calibration of computer models”. Journal of the Royal
Statistical Society: Series B (Statistical Methodology) 63 (3): 425–464.
Kleijnen, J. P. 1995. “Verification and validation of simulation models”. European journal of operational
research 82 (1): 145–162.
Livny, M., B. Melamed, and A. K. Tsiolis. 1993. “The impact of autocorrelation on queuing systems”.
Management Science 39 (3): 322–339.
Oakley, J. E., and A. O’Hagan. 2004. “Probabilistic sensitivity analysis of complex models: A Bayesian
approach”. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 66 (3): 751–769.
Pazgal, A. I., and S. Radas. 2008. “Comparison of customer balking and reneging behavior to queueing
theory predictions: An experimental study”. Computers & Operations Research 35 (8): 2537–2548.
Plumlee, M. 2016. “Bayesian calibration of inexact computer models”. Journal of the American Statistical
Association (just-accepted).
Plumlee, M., and H. Lam. 2016. “Learning stochastic model discrepancy”. In Winter Simulation Conference(WSC), 2016, 413–424. IEEE.
Sargent, G. R. 2013. “Verification and validation of simulation models”. Journal of Simulation 7 (1): 12–24.
Sargent, R. G. 2005. “Verification and validation of simulation models”. In Proceedings of the 37th
conference on Winter simulation, 130–143. winter simulation conference.
Shapiro, A. 2001. “On duality theory of conic linear problems”. In Semi-infinite programming, 135–165.
Springer.

Lam, Plumlee, and Zhang
Veeraraghavan, S., and L. Debo. 2009. “Joining Longer Queues: Information Externalities in Queue Choice”.
Manufacturing & Service Operations Management 11 (4): 543–562.
Wiesemann, W., D. Kuhn, and M. Sim. 2014. “Distributionally robust convex optimization”. Operations
Research 62 (6): 1358–1376.
AUTHOR BIOGRAPHIESHENRY LAM is an Assistant Professor in the Department of Industrial and Operations Engineering at the
University of Michigan (UM), Ann Arbor. He received his Ph.D. degree in statistics at Harvard University,and was an Assistant Professor in the Department of Mathematics and Statistics at Boston University before
joining UM. His research focuses on stochastic simulation, risk analysis, and simulation optimization. His
email address is khlam@umich.edu.
MATTHEW PLUMLEE is an Assistant Professor at the Department of Industrial and Operations Engineering at the University of Michigan. He received his Ph.D. degree from the School of Industrial andSystems Engineering at Georgia Tech. His interests lie in the interface of data and modeling, specificallyin methods for experimentation and uncertainty quantification for complex systems. His email address is
mplumlee@umich.edu.
XINYU ZHANG is a Ph.D. student in the Department of Industrial and Operations Engineering at the
University of Michigan (UM), Ann Arbor. He obtained his bachelor’s degree in UM, majoring in physics
and applied mathematics. His research interests include stochastic optimization and extreme event analysis.
His email address is xinyuzh@umich.edu.