Proceedings of the 2016 Winter Simulation Conference
T. M. K. Roeder, P. I. Frazier, R. Szechtman, E. Zhou, T. Huschka, and S. E. Chick, eds.

LEARNING STOCHASTIC MODEL DISCREPANCY
Matthew Plumlee
Henry Lam
Department of Industrial and Operations Engineering
University of Michigan
1205 Beal Ave
Ann Arbor, MI 48109, USA

ABSTRACTThe vast majority of stochastic simulation models are imperfect in that they fail to fully emulate the entiretyof real dynamics. Despite this, these imperfect models are still useful in practice, so long as one knows
how the model is inexact. This inexactness is measured by a discrepancy between the proposed stochastic
model and a true stochastic distribution across multiple values of some decision variables. In this paper, we
propose a method to learn the discrepancy of a stochastic simulation using data collected from the system
of interest. Our approach is a novel Bayesian framework that addresses the requirements for estimation ofprobability measures.
1

INTRODUCTION

Stochastic, or simulation, models are only approximations to the reality. A conjectured model may notalign with the true system because of unobserved complexity. Moreover, some highly accurate models,even if formulable, may not be implementable due to computational barriers and time constraints, in whichcase a simpler surrogate model is adopted. In all these cases, there is a discrepancy between the modeland the reality.
This important issue is traditionally handled through model validation and calibration in the simulation
literature. Model validation refers to the task of checking whether a developed stochastic model sufficiently
reflects the reality, often by comparing simulation outputs with real output data. A conventional method
compares the relevant outputs from these two sources statistically (e.g. Schruben 1980, Balci and Sargent1982). If the stochastic model is found to be inaccurate, ad-hoc calibration is performed to improve themodel recursively until the model is sufficiently consistent with the reality (Sargent 2005). While intuitive,
this conventional method may require building increasingly sophisticated models and it lacks a mechanism
to capture persistent discrepancies outside the structure conjectured by the model builders.
More broadly, the term calibration refers to the class of methods that can adjust the model to bring it
closer to reality. Model calibration using data has been well studied in the deterministic model case (e.g.,
Kennedy and O’Hagan 2001). General methodologies to handle the stochastic model case, however, have notgone much beyond the proposal to re-calibrate given an insufficient coherence between simulation and real
data revealed by testing. A unique property of the stochastic case is that we must account for the discrepancy
that considers the properties of probability distributions, which constitute the typical representation of theobserved outputs. Most existing statistical calibration methods in the deterministic domain, e.g. Higdon
et al. (2004), Higdon et al. (2008), Storlie et al. (2015), are thus not directly applicable. In the stochastic
domain, Goeva et al. (2014) studied calibration for the distribution of the simulation input model, but
assuming the simulation logic is correctly specified. A stream of work in the queueing inference literaturestudied the calibration of input processes and system performances (e.g., Whitt 1981, Basawa et al. 1996,
Pickands III and Stine 1997, Bingham and Pitts 1999, Hall and Park 2004, Larson 1990), given congestion

Plumlee and Lamor transaction data. These work are based on specific queueing structures that could be approximated either
analytically or via diffusion limits, and, moreover, that are not subject to misspecification errors.
This paper studies a general framework to learn, and correct, for discrepancies between stochastic
outputs from simulation models and the reality. In contrast to the past literature, we do not restrict thesimulation logic or presume the logic matches reality. Instead, we focus on situations where responses fromthe real system are recorded for multiple scenarios. Taking into account the discrepancies between the real
data and the model outputs at these different scenarios, we seek to give an estimate for the performancemeasure of interest at a described scenario. Our key methodological contribution is to build a Bayesianlearning framework that operates on the probabilistic representation inherent in simulation models by using
likelihood ratios as our main inference objects.
Pedagogical example To fix the conceptual ideas in this paper, we borrow from the call center data
originally analyzed in Brown et al. (2005). This dataset is associated with a call center where a customer isplaced into a queue once a call is received until one of n servers is available. From these data, the sample
mean of the waiting time (from entry to service for a customer) from 9:00 to 10:00 am is calculated. In thisnarrow time period, the arrival rate and service rate, which is time inhomogenous according to Brown et al.
(2005), should be approximately homogenous. The average arrival rate in this data was 1.76 customers
per minute and the average service time was 3.72 minutes per customer. In this paper, we also account for
the number of servers operating in the system at any given time, which appears to differ between days. To
our reading, this subset of the data was by-and-large ignored in Brown et al. (2005)’s original analysis.
Our model for this call center will be an n-server first-come-first-serve queue. Following practice,both the interarrival and service times will be modelled as exponentially distributed. After a warm up
period, the sample average of the waiting time is measured over the course of a one hour window. This,in principal, agrees with Brown et al. (2005). Moreover, this modeling framework likely emulates many
common practices. In the spirit of ad-hoc calibration, two additional features were added: (i) the arrival
rate is randomly generated each day from a log normal distribution with associated mean parameter 1.8
and variance 0.4 and (ii) with probability 2/3, a customer will abandon the queue if the wait time is longer
than an exponential random variable with mean 1.5, otherwise they will remain in the queue until served.
Adding both of these features resulted in a simulation model that was closer to the observed data.
When this example is revisited in this article, it will become clear that the simulation model output does
not exactly agree with the collected responses. Even though the parameters of the model are well-estimatedand the model has been deliberately calibrated, the model is not perfect thus it does not agree with thecollected responses from the real system. This creates the problem in practice: how does one go about
calibrating this model in a precise way using these observations?
2

GENERAL SETTING

In this paper, the system of interest outputs a discrete (or discretized) response over a space Y with cardinality
q. This response depends on a vector of design variables, denoted x, and each potential combination areindexed from i = 1, . . . , m, where m is the total feasible combinations of these variables. The variable x canbe broadly defined to include input variables that can not necessarily be “controlled”. An example of the
response is the waiting time in call centers (which is the running pedagogical example) or hospitals (Helmand Van Oyen 2014). In the first example, design variables could be the number of servers, the system
capacity, and the arrival rate. In the second example, the design variables could be the rate of elective
admissions. The probability measure πi describes the distribution of the response of the real system. The
term design point will be understood to mean the point in the space corresponding to a combination ofdesign variables.
The objective is to make decisions about the distribution of the response’s distribution πi for all ibetween 1 and m. In reality, a user rarely has direct access to the true distribution functions. Perhaps thesimplest tool to avoid this problem is by directly observing many independent responses from the system,

Plumlee and Lam
possibly from a designed experiment (Li et al. 2015). Using this data,πbi (γ) =1 Ni∑ I (yik = γ)Ni k=1(1)

is an estimate of the probability mass function for the true response. Here, I is the indicator function andyik is the kth recorded response on the ith setting. For each of the m design points, there are Ni recorded
responses corresponding to the ith decision variable. As Ni grows, πbi (γ) will converge to the true massfunction, πi (γ), under reasonable conditions.
In practice, this purely empirical estimate is not sufficient to make decisions. A user needs many
samples in order to directly estimate this quantity to sufficient precision with πbi (γ) alone. This problem
is compounded because doing this for all possible xi is almost never a reasonable option. A user may
not have access to the some of the design points. For example, this can happen when a potential design
point has been proposed but has not yet been implemented. So for many i, Ni is 0, which makes directly
estimating πi purely from observed responses impossible.
Instead of attempting to directly conduct inference on the unknown πi (·), operation researchers oftenchoose a surrogate model to use instead. These models are described using state-of-the-art understanding ofa given system, perhaps simplified for computational reasons. The example mentioned in the introductionis the use of first-come first-serve multi-server queues with a Poisson arrival process and an exponentiallydistributed service times. Despite being a reasonable model, it shows obvious discrepancies with the realdata.
Let π̃i be the surrogate model’s distribution associated with the design point xi . For simplicity, assume
that the surrogate model is exactly provided in the sense that simulation is cheap relative to the cost ofphysical experimentation, and that all parameters in this model are fixed. If the surrogate model is assumed
exactly equal to the true model, then πi (γ) can be found precisely with the relation πi (γ) = π̃i (γ), andno data is needed. But, as discussed before, there is almost always a discrepancy between πi and π̃i inpractice. Thus, in this paper, we will operate from the assumption that the surrogate model is inexact,
and that data is collected from the system. In general, this means that πi (γ) 6= π̃i (γ) for at least someγ ∈ Y and some i ∈ {1, . . . , m}, where πi needs to be statistically inferred. As of yet, there is no systematic
way to address this uncertainty and learn the inexactness to the authors’ knowledge. This motivated ourapproach to automatically correct the deficiencies using data, by viewing the model as a black box. The
next section will explain the definition of discrepancy in simulation models so that the proposed correction
and uncertainty quantification scheme is well defined. Before that, we revisit our running example:
Pedagogical example (continued) The response is discretized into the four categories < 1, 1 − 2, 2 − 3,
> 3 and thus q = 4. In this example, 5 − 9 servers were considered and thus m = 5. We also have π̃ to areasonable approximation by simulating from the described model 1000 times for each design point. Figure
1 shows the frequency diagram of the data along with the surrogate model. The data we use here does
not include 5 servers and 9 servers. Using this simulation model, we can conjecture what would happen
at these staffing levels. Going by the simulation model alone, one would expect that when 5 servers arepresent, the average waiting time should often be larger than 3, and when 9 servers are present, the average
waiting time will almost never be above 1.
But when looking at the collected responses next to the simulation model, the picture is more muddled.
According to the simulation model, the chance of the average waiting time between 1 and 2 minutes
is around 8% when 8 servers are used. But in our collected responses, about 25% of the 23 recorded
observations are in this range. This indicates that the model is inexact. This knowledge should affectour prediction of the probability of the response being between 1 and 2 when 9 servers are present. Our
methodology presented in the next sections will explicitly bring such knowledge into consideration.

Plumlee and Lam
5 servers, n1 = 0

Probability

6 servers, n2 = 64

7 servers, n3 = 38

0.8

0.5

8 servers, n4 = 23
1

0.8
0.6

0.4

0.6
0.3

0.4

0.4

0.2
0.1
0
0

1

2

3

Average waiting time

0.8

0.8

0.6

0.6

0.4

0.4
0.2

0.2

0.2

0.2

0

0

0

0

1

2

3

Average waiting time

0

1

2

3

9 servers, n5 = 0
1

00

Average waiting time

1

2

3

Average waiting time

0

1

2

Figure 1: Frequency diagram of the real data (the left bars) and the probabilities computed from thesurrogate model (the right bars) for the pedagogical example.
3

DISCREPANCIES AND VALID DISCREPANCIESThis article suggests directly of modeling the discrepancy between πi and π̃i . The discrepancy δi is definedasπi (γ)
δi (γ) =, for all γ ∈ Y .
(2)π̃i (γ)i.e., the ratio between πi and π̃i . A value of δi (γ) that is exactly 1 for all γ ∈ Y would represent a surrogate
model with no discrepancy. The δi reflects the the relative probability of the response between the true
response and simulated response. Model inexactness can then be cast as δi (γ) 6= 1 for some γ ∈ Y and i
between 1 and m.
There are direct analogies between δi and the likelihood ratio (or adjustment factor) in the context ofimportance sampling (e.g., Owen (2013), Chapter 9; Asmussen and Glynn (2007), Chapter IV; Glasserman(2003), Chapter 4). However, whereas importance sampling is used to speed up Monte Carlo schemes
under a known stochastic model, and the associated πi and δi in that situation are therefore known andare designed by the user, ours are not fully observable. In the uncertainty quantification context, δi also
appears in worst-case optimizations formulated to quantify model errors. These optimizations computeconservative estimates of performance measures subject to beliefs or uncertainty on the true model relative
to the surrogate (often known as the baseline model). For instance, Hu et al. (2012) considered Gaussianmodels with mean and covariance uncertainty represented by linear matrix inequalities. Glasserman andXu (2014) considered nonparametric uncertainty measured by Kullback-Leibler divergence. A notabledifference between these lines of work and ours is the role of δi . While the former is a decision variable
in an optimization program that computes a worst-case bound, ours is an object to be inferred from data
to improve the prediction accuracy.
For simplicity, we assume that our surrogate probabilities π̃i (γ) are all strictly positive, thus this valueof δi (γ) will naturally be ≥ 0 and finite. Any π̃i (γ) for which this is not naturally true can be slightly
adjusted by adding small probability for all γ ∈ Y and renormalizing the distribution such that the masses
sum to one.
In addition to the fact that 0 ≤ δi (y) < ∞, there is another important constraint on δi that relates to thefact that πi is a probability distribution. Because∑ πi (γ) = 1, then ∑ δi (γ)π̃i (γ) = 1.
γ∈Y

γ∈Y

Thus the L2 inner product between δi and π̃i over Y must be exactly equal to 1 for all i from 1 to m.
From the proceeding discussion,
Definition Say π̃i is a probability distribution over the discrete sample space Y with π̃i (γ) > 0 for allγ ∈ Y . Any mapping h : Y → R is a valid discrepancy with respect to π̃i if(i)

0 ≤ h(γ) < ∞ for all γ ∈ Y and3

Average waiting time

Plumlee and Lam(ii)∑γ∈Y h(γ)π̃i (γ) = 1.

All δi as defined in (2) will meet the criteria outlined above, regardless of the true distribution (provided
the true distribution is indeed a proper probability distribution over Y ).
These two conditions are consistent with the traditional assumptions on likelihood ratios. We highlightthem here because they characterize the type of objects that we conduct inference on. Moreover, theyshow a distinction between inferring stochastic and deterministic model discrepancies, in that probabilistic
structure plays a role in the former but not the latter (e.g. compare to Kennedy and O’Hagan (2001)).
We also comment at this point that the finite support assumption of the output distribution, hence thediscrepancy, avoids the technical difficulties in both the statistical inference and the computation involving
infinite-dimensional parameters. For continuous distributions, one approach is to project them onto a finite
number of statistics, such as moments. Such investigation is out of the scope of this work.
4

LEARNING DISCREPANCY

Calibration problems are often located in data-poor environments (in terms of observed responses) which
are also rich in prior knowledge. Beyond the notion of a valid discrepancy provided in Section 3, there
is often other information available. Even though observed responses do not exist for every design point,
similar discrepancies might be anticipated for similar design points. This similarity is measured by the
distance in terms of the design variables associated with each point. Recall the example of waiting times
in call centers; the discrepancy when the number of servers is 5 is expected to be more similar to thediscrepancy when the number of servers is 6 relative to when the number of servers is 8.
To incorporate this type of prior information in a data-poor environment, this paper will use a Bayesianlearning approach. Under this framework, the observed responses have a certain likelihood of occurring
given the knowledge of the discrepancy. Combining the likelihood of these observed responses with theprior information creates a posterior for the discrepancy.
More concretely, we let p generally represent a probability density and p(a|b) implies the conditional
probability density of a given b. Let a vector d represent a vectorized version of the discrepancy:d = [δ1 (γ1 ), . . . , δ1 (γq ), · · · , δm (γ1 ), . . . , δm (γq )]T .

(3)

The term data substitutes for the event of all observed responses from the true system:
data = Yi j = yi j , i = 1, . . . , m, j = 1, . . . , Ni .
Our objective is thus to conjecture onp (d |data ) ,which would define the posterior of d. Using Bayes rule,p (d |data ) ∝ p(data|d)p(d),(4)

where ∝ stands for equality up to a constant multiplier. The computational methods that we will propose
later require only this proportional value.
The first element in the right hand side of (4) is referred to as the likelihood and the second element in
the right hand side of (4) is referred to as the prior distribution. The likelihood is straightfowardly defined
based on π. Our prior distribution will be such thatp (d = g + Gθ ) ∝ exp ε1T log(g + Gθ ) − λ θ T θ ,where ε is a small positive constant, λ is a tuning parameter left to the user, θ ∈ R pq relates to d via anaffine map d = g + Gθ with g and G designed such that we meet the criteria of a valid discrepancy. BasedPlumlee and Lamon this prior distribution, our posterior on θ isp (θ |data ) ∝ exp (N + ε)T log(g + Gθ ) − λ θ T θ ,(5)and we can recover d with d = g + Gθ . Note that N is the length mq vector of counts
"
N=

n1

n1

nm

nm

#T

∑ I(y1 j = γ1 ), . . . , ∑ I(y1 j = γq ), · · · , ∑ I(ym j = γ1 ), . . . , ∑ I(ym j = γq )

j=1

j=1

j=1

.

j=1

We explain how the chosen prior distribution agrees with the concept of a valid discrepancy. Basically, it
can be interpreted as a multivariate Normal prior on the discrepancy subject to the probabilistic constraints.
The vector g and matrix G are derived from a positive definite matrix R and a mean vector µ. These R andµ can be considered like the mean and covariance for the prior on the bias. For example, it is reasonable
to use a vector of 1s as the prior mean for µ. But a prior distribution on d of Normal with a mean vectorof 1s and a covariance matrix that is proportional to R will not satisfy the probabilistic definition of a valid
discrepancy with probability one in general. To resolve this, we let g and G correspond to−1
g = µ + RΠ̃T Π̃RΠ̃T(1 − Π̃µ)
and

q
−1
Π̃R.
G = R − RΠ̃T Π̃RΠ̃T

These can be thought of as the conditional mean and square root of the covariance after conditioning on(ii) in the definition of a valid discrepancy. To the authors’ knowledge, this method of conditioning normal
distributions to resolve the probabilistic definition of a valid discrepancy has not been used elsewhere in
the Bayesian statistical literature.
We also need to consider condition (i) when assigning a prior distribution, which is taken care of by
the inclusion of the ε in our prior. This can be considered as replacing a sharp constraint with what are
known in optimization as barrier functions (Ben-Tal and Nemirovski 2001, pp 274). That is, I(δi (γ j )) ≥ 0)is replaced with the barrierexp(ε log(δi (γ j ))),
where ε > 0 is some small constant.
Thus we conclude the following result based on the described prior distribution:
Theorem 1 Say that π̃i (γ) > 0 for all i ∈ {1, . . . , m} and γ ∈ Y . There exists a constant C such that
C exp ε1T log(g + Gθ ) − λ θ T θ ,is a valid probability density on Rmq . Moreover, if θ is drawn from this distribution on Rmq , then withprobability one(i) g + Gθ > 0 and (ii) Π̃ (g + Gθ ) = 1.
Pedagogical example (continued) Consider our pedagogical example described in Sections 1 and 2. Figure
2 shows the results of the Bayesian analysis of the data. The top panels demonstrate what would happen
if the probabilistic constraints that control the notion of a valid discrepancy were ignored in our analysis.
Clearly, the non-sensical results with probabilities larger than 1 would not meet any criteria for decision
making. While not present here, sometimes these bounds can dip below zero, inducing severe inferential
problems. Using the proposed approach that takes into account probabilistic constraints, these issues appear
resolved as illustrated in the bottom panels.

Plumlee and Lam
Calibrated, invalid
Probability

5 servers, n1 = 0
3

7 servers, n3 = 38

8 servers, n4 = 2344

9 servers, n5 = 0

6

4

4

3

6
2

2
1

2

0

0
0

Calibrated, valid
Probability

6 servers, n2 = 64
6

1

2

3

0

1

2

2

2

0

0

3

0

1

2

3

0.8

1
0
0

1

2

3

1

0.6

0

1

2

3

1

2

3

1

0.8
0.6
0.6

0.4
0.4
0.2
0
1

2

3

Average waiting time

0.5

0.2

0
0

0.5

0.4

0.2

0
0

1

2

3

Average waiting time

0
0

1

2

3

Average waiting time

0
0

1

2

3

0

Average waiting time

Average waiting time

Figure 2: Frequency diagram of the observed responses and the probabilities computed from the surrogate
model for the pedagogical example with predictive intervals marked by the rectangles for π based on theobserved responses, the surrogate model, and our prior information. The top panels use prior information that
ignores validity and the bottom panels use prior information that includes the notion of a valid discrepancy.
A few issues were brought up in the discussion of this example in Section 2. Specifically, when there
were 9 servers, how likely was it that the average waiting time would be between 1 and 2 minutes? Thereis no data, but the simulation model combined with the observed responses and our prior information gives
us an estimate of somewhere between 2.5% and 8%. This accounts for the discrepancy that we observed
based on the recorded responses at 8 servers. Overall, this range appears to agree with both the data andthe surrogate outputs. We are not confident, for example, that staffing 5 servers will produce the same
results as the simulation model, which has average waiting times over 3 minutes about 15% of the time.
Based on the recorded responses, this could be 22%, but it could also be as low as about 2%. More data
would produce tighter bounds; see the result with 64 data points when the number of servers is 6. The
confidence intervals are much thinner due to more available data.
5

POSTERIOR COMPUTATIONSThis section discusses Bayesian computation methods for a quantity of interest, described by a functionζ (·, ·) : X × Y → R, from the posterior. That is, we find upper and lower bounds formm∑ ∑ ζ (xi , γ)πi (γ) = ∑ ∑ ζ (xi , γ)π̃i (γ)δi (γ).

i=1 γ∈Y

i=1 γ∈YLetz̃ = [ζ (x1 , γ1 )π̃1 (γ1 ), . . . , ζ (x1 , γq )π̃1 (γq ), · · · , ζ (xm , γ1 )π̃m (γ1 ), . . . , ζ (xm , γq )π̃m (γq )] .
5.1 Sampling-based approach
Given w samples from the length mq vectors of the discrepancy parameter θ (1) , . . . , θ (w) , samples from theposterior of the quantity of interest are{z̃T g + z̃T Gθ (k) ,k = 1, . . . , w}.

The 95% credible interval is then the gap between the 2.5% and 97.5% empirical quantiles.
Because of the nonregularity of the posterior (and the fact that there is an unknown proportionality
constant), direct Monte Carlo sampling is not possible. Sophisticated Markov chain Monte Carlo samplers

Plumlee and Lam
were designed explicitly for this purpose (Gelman et al. 2014, Chapters 11 and 12). Perhaps the simplest
mechanism for sampling from these nonregular type is a simple Metropolis Hastings algorithm with asymmetric proposal (Gelman et al. 2014, pp 278-280). Here, some initial θ (0) is chosen, typically as the
maximizer of the posterior in (5). Then an arbitrary symmetric distribution is selected such that the next
point is θ (0) + ρ (0) , where ρ (0) is drawn from this symmetric distribution. Starting with i = 0, the next
sample from step i is


(
(i)
(i)
(i) + ρ (i) , with probability min 1, p(θ +ρ |data)θ(i)p(θ |data)θ (i+1) =θ (i) , otherwise.
In the above ρ (i) are independent draws from this symmetric distribution.
Another approach to consider, which has grown in popularity, is the idea of using Hamiltonian dynamics
to modify the Metropolis Hastings algorithm with a non-symmetric proposal distribution. The key benefitof this approach is that it considers the derivative of the posterior when constructing proposal samples. The
leap-frog algorithm, taken from Neal (2011), pp 121) and modified here for our specific density function,is that given θ (i) the next sample will be:
Algorithm 1 Hamiltonian Monte Carlo Sampler
θ prop = θ (i)Draw v0 as a mq vector of independent standard
Gaussian random variables,

1
T
prop
−1
Let v = v0 + 2 (N + ε1) (g + Gθ )G − λ θ propfor L times doθ prop = θ prop + ev
v = v + (N + ε1)T (g + Gθ prop )−1 G− 2λ θ prop .
1θ (prop) , with probability min 1, p(θ (prop) |data) exp(kv+ 2 ((N+ε)T (g+Gθ prop )−1 )G−λ θ prop k2 /2)
exp(kv0 k2 /2)p(θ (i) |data)θ (i+1) = (i)θ , otherwise.
The numerical constant e > 0 and the integer L > 0 are parameters of this algorithm that can be adjusted
to reach a distribution more inline with the results in the problem.
Lastly, consider slice sampling (Neal 2003). This sampler, similar to principal to the Hit and Run
algorithm (Bélisle et al. 1993) and the approach of Chen and Schmeiser (1998). Because of the uniqueproperties of our posterior (namely, convexity of the level sets), the reflective slice sampling has particularappeal. The basic idea is to sample uniformly from level-sets of the density by choosing a direction andreflecting off the barrier of the level set. Start with some positive constant α and let θ ∗ be the maximizerof the posterior. The algorithm is that given θ (i) , α > 0 the next sample will be:
Algorithm 2 Slice Sampler
Set θ (i+1) = θ (i) .
Draw random direction hfor L times do
Find largest step size ∆ along h such that p(θ (i+1) + ∆h|data) ≥ p(θ ∗ |data) exp(−α).
Set θ (i+1) = θ (i+1) + ∆h.

TLet v = (N + ε1)T (g + Gθ (i+1) )−1 G − 2λ θ (i+1) and update h = h − 2v vvThv .
Find largest step size ∆ along h such that p(θ (i+1) + ∆h|data) ≥ p(θ ∗ |data) exp(−α).
Let U be a uniform random variable between 0 and 1 and θ (i+1) = θ (i+1) +U∆h.
Let α = log p(θ ∗ |data) − log p(θ (i+1) |data) − log(U).

Plumlee and Lam
The numerical constant L > 0 is a parameter of this algorithm that can be adjusted. There is no needfor an acceptance/rejection step. The major computational the burden of this sampler lies in our ability tofind the maximum step size quickly and effectively.
Pedagogical example (continued) We attempted the samplers from Section 5.1 on the posterior describedin the running example. The results are located in Table 1. Attempts were made to optimize all threemethods and balance the computational time by using different total amounts of samples collected. Half
the samples for each sampler were discarded as a warm-up period. Though these results are limited, they
indicate that the choice of sampler will play a major role in being able to efficiently sample from theposterior of the quantity of interest. The method with the best effective sample size on a subsample of
200 was the slice sampling approach. This is despite that the number of iterations is low compared to theother two methods.
Table 1: The results for the sampling approaches using the queueing data as described in Section 5.1.
Total time (seconds)
Recorded samples
ESS% rejectedMH
0.51
10,00048/20044.6 %

HMC
0.61
800
107/200
18.0 %

SS
1.05
400
198/200
0.0%

5.2 Optimization-based approach
Because of the complexity in the setup and execution of the above samples, we also propose a new method
for posterior inference using convex optimization to develop confidence bounds for the quantity of interest.
This approach is inspired from the literature of robust optimization (e.g., Ben-Tal and Nemirovski 2002,
Ben-Tal et al. 2009, Bertsimas et al. 2011) that represents uncertain or ambiguous parameters in an
optimization problem by imposing constraints, forming the so-called ambiguity or uncertainty sets. Whenthe uncertain parameters follow probability distributions, the uncertainty sets are constructed as a predictive
set that contains the truth with a prescribed confidence. This is particularly useful in approximating chanceconstrained programs (e.g., Ben-Tal and Nemirovski 2002, Chapter 2), and performance measures fromcomplex stochastic models (e.g., Bandi and Bertsimas 2012, Bandi and Bertsimas 2014). Recently, Gupta(2015) considered obtaining Bayesian statistical guarantees for uncertainty sets, particularly for the case of
distributionally robust optimization where the uncertain parameters are themselves probability distributions.
Beyond that, to our best knowledge, there has been no attempt to use robust optimization as a tool forcomputing Bayesian posterior. Thanks to the probabilistic structure and the high dimensionality of theparameter space (qm) induced in our setting, we advocate robust optimization as an efficient computational
tool for inferring on the posterior distribution.
To begin, consider an uncertainty set of d based on the posterior. This uncertainty set replaces theprobabilistic posterior with a deterministic set-based statement of the feasibility of the discrepancy. Thatis, given some constant c that controls the “size” of the uncertainty set, let
U = {d = g + Gθ |p(θ |data) ≥ c, θ ∈ Rqm } .
Clearly, as c gets smaller, U is nondecreasing in size. Motivated by some theoretical results excluded
from this article, a reasonable approach would be to set
1
c = Φ−1 (w)2 + maxqm p(θ |data).
θ ∈R
2
The benefits of this choice are two fold. First, this guarantees that U will be nonempty. Second, this choiceleads to asymptotic bounds very close to 95% confidence bounds, which are popular in diverse areas ofscience and business.

Plumlee and Lam
Our optimization problem to find the bounds is stated as
max or min z̃T d,
subject to
d∈U.

(6)

which can be solved with most convex programming approaches if the feasible region is convex (such as
penalty methods and interior point methods; see Bertsekas (1999)). Because of the structure of the problem,
Proposition U is a strictly convex set on the interior of Rqm+ for all ε > 0 and λ > 0.
Note that the assumption of the prior plays a critical role. If ε was not included in the prior, we would
not be guaranteed that U is a bounded set on the interior of valid discrepancies (by the coercivity ofp(·|data)). The prior for θ was chosen to be log-concave. If a different choice was made for the prior onθ , U may not be convex.
Pedagogical example (continued) We now offer a graphical explanation of the differences in the methods
described in Sections 5.1 and 5.2. Say a user is concerned with staffing exactly 6 servers. What is anupper bound on the probability that the average wait time is below 2 minutes? Using the discretization
introduced for this example, this would correspond toπ2 (γ1 ) + π2 (γ2 ) = π̃2 (γ1 )δ2 (γ1 ) + π̃2 (γ2 )δ2 (γ2 ) = 0.251δ2 (γ1 ) + 0.294δ2 (γ2 ).
i.e., ζ (x2 , γ1 ) = ζ (x2 , γ2 ) = 1, and all other ζ ’s are 0.
Consider the sampling based approach. If 200 outputs are sampled from θ (1) , . . . , θ (200) , the upper
credible interval would be determined by finding the hyperplane parallel to 0.251δ2 (γ1 ) + 0.294δ2 (γ2 ) suchthat exactly 5 samples are to the right and above it. The robust optimization approach works similarly,
with the major difference that the hyperplane is maximized over the convex set U described in Section
5.2. This is illustrated in Figure 3, where it should be emphasized that these are both the samples and thetrue uncertainty set computed from the data, not merely a cartoon. Here, the convexity (and nonellipsoidal)
shape of the projection of the uncertainty set is verified visually. The samples appear to roughly agree withthe uncertainty set despite them being generated from an approximate sampler.
6

Conclusion

We have studied inference on model discrepancies between surrogate stochastic or simulation models andreal response data under the availability of multiple design points. This probabilistic constraints of this caseis unique to stochastic settings and deviates from previous work in the deterministic setting where no such
constraint is considered, and in other stochastic calibration problems where specific model structures are
assumed known or perfectly specified. This articles has also compared two ways of computing posterior
predictive intervals, typical samplers and a new optimization approach. We anticipate that learning stochastic
model discrepancy can be extended to higher-dimensional problems (such as continuous, and thus infinitedimensional, distributions). The major complexity here is that the number of positivity constraints is on
the order of the space space, thus any proposed method must be careful designed to be well-defined andcomputable.
ACKNOWLEDGMENTS
We thank Ilan Guedj for the data organization and for Avi Mandelbaum for continuing to place the data on
the website http://ie.technion.ac.il/serveng/. We also gratefully acknowledge support from the the National
Science Foundation under grants CMMI-1400391/1542020 and CMMI-1436247/1523453.
REFERENCES
Asmussen, S., and P. W. Glynn. 2007. Stochastic Simulation: Algorithms and Analysis. Springer Science& Business Media.

discrepancy @ [1, 2)

Plumlee and Lam

2

O pt i

1.8

97.5%

1.6

mi z a

tion

maxi

quan
tile

2.5% Uquan
tile
O pt imi z a
tion
mi

1.4
1.2
0.65

0.7

0.75

0.8

0.85

0.9

mum

nimu

0.95

m
1

1.05

1.1

discrepancy @ [0, 1)
Figure 3: Graphical description of the differences between optimization- and sampling-based approaches forthe stochastic discrepancy problem. The bottom axis corresponds to δ2 (γ1 ) and the vertical axis corresponds
to δ2 (γ2 ). The light arrows describe the gradient of the objective function. The 200 dots are samples froma well-tuned slice sampler. The upper limit of the 2.5% and 97.5% quantile is indicated by the dashedlines line. The region labeled U is the projection of the convex set U onto this two-dimensional plane.
The maximum of the optimization is determined by the pictured line.
Balci, O., and R. G. Sargent. 1982. “Some examples of simulation model validation using hypothesis
testing”. In Proceedings of the 1982 Winter Simulation Conference, edited by Highland, Chao, andMadrigal, 621–629. Piscataway, NJ: IEEE.
Bandi, C., and D. Bertsimas. 2012. “Tractable Stochastic analysis in high dimensions via robust optimization”.
Mathematical Programming 134 (1): 23–70.
Bandi, C., and D. Bertsimas. 2014. “Robust option pricing”. European Journal of Operational Research 239(3): 842–853.
Basawa, I. V., U. N. Bhat, and R. Lund. 1996. “Maximum likelihood estimation for single server queues
from waiting time data”. Queueing Systems 24 (1-4): 155–167.
Bélisle, C. J., H. E. Romeijn, and R. L. Smith. 1993. “Hit-and-run algorithms for generating multivariate
distributions”. Mathematics of Operations Research 18 (2): 255–266.
Ben-Tal, A., L. El Ghaoui, and A. Nemirovski. 2009. Robust Optimization. Princeton University Press.
Ben-Tal, A., and A. Nemirovski. 2001. Lectures on Modern Convex Optimization: Analysis, Algorithms,
and Engineering Applications. SIAM.
Ben-Tal, A., and A. Nemirovski. 2002. “Robust optimization–methodology and applications”. Mathematical
Programming 92 (3): 453–480.
Bertsekas, D. P. 1999. Nonlinear Programming. Athena Scientific Belmont.
Bertsimas, D., D. B. Brown, and C. Caramanis. 2011. “Theory and applications of robust optimization”.
SIAM Review 53 (3): 464–501.
Bingham, N., and S. M. Pitts. 1999. “Non-parametric estimation for the M/G/∞ queue”. Annals of theInstitute of Statistical Mathematics 51 (1): 71–97.
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 Statistical
Association 100 (469): 36–50.
Chen, M.-H., and B. Schmeiser. 1998. “Toward black-box sampling: A random-direction interior-point
Markov chain approach”. Journal of Computational and Graphical Statistics 7 (1): 1–22.

Plumlee and Lam
Gelman, A., J. B. Carlin, H. S. Stern, and D. B. Rubin. 2014. Bayesian Data Analysis. 3 ed. Taylor &Francis.
Glasserman, P. 2003. Monte Carlo Methods in Financial Engineering. Springer Science & Business Media.
Glasserman, P., and X. Xu. 2014. “Robust risk measurement and model risk”. Quantitative Finance 14 (1):29–58.
Goeva, A., H. Lam, and B. Zhang. 2014. “Reconstructing input models via simulation optimization”. In
Proceedings of the 2014 Winter Simulation Conference, edited by A. Tolk, S. D. Diallo, I. O. Ryzhov,
L. Yilmaz, S. Buckley, and J. A. Miller, 698–709. Piscataway, NJ: IEEE.
Gupta, V. 2015. “Near-Optimal Ambiguity Sets for Distributionally Robust Optimization”. Technical report.
Hall, P., and J. Park. 2004. “Nonparametric inference about service time distribution from indirect measurements”. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 66 (4): 861–875.
Helm, J. E., and M. P. Van Oyen. 2014. “Design and optimization methods for elective hospital admissions”.
Operations Research 62 (6): 1265–1282.
Higdon, D., J. Gattiker, B. Williams, and M. Rightley. 2008. “Computer Model Calibration Using HighDimensional Output”. Journal of the American Statistical Association 103 (482): 570–583.
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.
Hu, Z., J. Cao, and L. J. Hong. 2012. “Robust simulation of global warming policies using the DICE
model”. Management Science 58 (12): 2190–2206.
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.
Larson, R. C. 1990. “The queue inference engine: Deducing queue statistics from transactional data”.
Management Science 36 (5): 586–601.
Li, J. Q., P. Rusmevichientong, D. Simester, J. N. Tsitsiklis, and S. I. Zoumpoulis. 2015. “The value offield experiments”. Management Science 61 (7): 1722–1740.
Neal, R. M. 2003. “Slice sampling”. Ann. Statist. 31 (3): 705–767.
Neal, R. M. 2011. “MCMC using Hamiltonian dynamics”. Handbook of Markov Chain Monte Carlo:113–
162.
Owen, A. B. 2013. Monte Carlo Theory, Methods and Examples.
Pickands III, J., and R. A. Stine. 1997. “Estimation for an M/G/∞ queue with incomplete information”.
Biometrika 84 (2): 295–308.
Sargent, R. G. 2005. “Verification and validation of simulation models”. In Proceedings of the 2005 Winter
Simulation Conference, edited by M. E. Kuhl, N. M. Steiger, F. B. Armstrong, and J. A. Joines, 130–143.
Piscataway, NJ: IEEE.
Schruben, L. W. 1980. “Establishing the credibility of simulations”. Simulation 34 (3): 101–105.
Storlie, C. B., W. A. Lane, E. M. Ryan, J. R. Gattiker, and D. M. Higdon. 2015. “Calibration of computationalmodels with categorical parameters and correlated outputs via Bayesian smoothing spline ANOVA”.
Journal of the American Statistical Association 110 (509): 68–82.
Whitt, W. 1981. “Approximating a point process by a renewal process: The view through a queue, anindirect approach”. Management Science 27 (6): 619–636.
AUTHOR BIOGRAPHIES
MATTHEW PLUMLEE (mplumlee@umich.edu) is an Assistant Professor researching simulation experiments and calibration at the Department of Industrial and Operations Engineering at the University ofMichigan.
HENRY LAM (khlam@umich.edu) is an Assistant Professor researching stochastic simulation, risk analysis,
and simulation optimization at the Department of Industrial and Operations Engineering at the Universityof Michigan.