Proceedings of the 2020 Winter Simulation Conference
K.-H. Bae, B. Feng, S. Kim, S. Lazarova-Molnar, Z. Zheng, T. Roeder, and R. Thiesing, eds.

CALIBRATING INPUT PARAMETERS VIA ELIGIBILITY SETS
Yuanlu Bai
Henry Lam

Industrial Engineering and Operations Research
Columbia University
500 West 120th Street
New York, NY 10027, USA

ABSTRACT
Reliable simulation analysis requires accurately calibrating input model parameters. While there has beena sizable literature on parameter calibration that utilizes directly observed data, much less attention has
been paid to the situation where only output-level data are available to justify input parameter choices.
This latter problem, which is known as the inverse problem and relates to the model validation literature,
involves several new challenges, one of which is the non-identifiability issue. In this paper we introduce
the concept of eligibility set to bypass non-identifiability, by relaxing the need of consistent estimation to
obtaining bounds on the input parameter values. We reason this concept from the worst-case notion inrobust optimization, and demonstrate how to compute eligibility set via empirical matching between thesimulated and the real outputs. We substantiate our procedure with theoretical error analysis and validate
its effectiveness via numerical experiments.
1

INTRODUCTION

Stochastic simulation is widely-used in system performance evaluation and decision-making. It operates
by generating random variates, fed through the system logic, to obtain random outputs that are utilized for
subsequent analysis. To ensure the reliability of the analysis, the input parameters in the simulation need
to be accurately calibrated. These parameters can arise in the probability distributions of the input randomvariates and the system logic, and are typically obtained via data or domain knowledge. The literature
on input modeling and input uncertainty, which has grown substantially in recent years, precisely handles
these issues (see, for example, the surveys Henderson (2003), Chick (2006), Barton (2012), Song et al.
(2014), Lam (2016)).
While calibrating input parameters from direct data has been actively studied, much less attention has
been paid to situations where only output-level data, instead of direct input-level data, are available toassist input parameter selection. This nonetheless happens ubiquitously in practice. As a simple example,
models of inter-arrival times and service times are needed to simulate generic service queues. Some service
centers, however, may not have advanced electronic systems to track each customer at all times, and may
only record the times a customer arrives or leaves the building. It is thus unclear how the sojourn time iscontributed between the waiting time and the service time, yet the latter is needed to calibrate the model
for use. One could imagine analogous issues can get substantially more complicated as the system size
grows.
In this paper, we study input parameter calibration using only output-level data. In scientific modeling,
the task of inferring input from output is often referred to as an “inverse problem” (Tarantola 2005). Related
investigations include model error inference or uncertainty quantification, typically addressed via Bayesian
approaches (Kennedy and O’Hagan 2001; Currin et al. 1991). In stochastic simulation, these issues areBai and Lam
closely related to so-called model calibration and validation (Sargent 2010; Kleijnen 1995). Upon building
a simulation model, the modeler checks the differences of the simulated outputs with real-world outputs(“validation”), and makes refinements to the simulation model to close the gap (“calibration”). This process
is conducted iteratively, each time possibly with expansive effort to collect more real data, and using expert(or otherwise ad hoc) judgements to decide how to validate and what data to collect. Some statistical tests
have been suggested for validation, such as the two-sample mean-difference tests (Balci and Sargent 1982)
and the Schruben-Turing test (Schruben 1980). However, despite these works and the importance of model
calibration and validation, there appears relatively few studies on concrete general-purpose methodologies
for these tasks (Nelson 2016).
A challenge in inverse problem or model calibration is the so-called non-identifiability (Tarantola 2005).
That is, there exist more than one input parameter value that exactly matches the output behavior. This is
an intrinsic barrier in inverse model calibration, in the sense that this cannot be resolved even if one has an
infinite amount of data. It stems, instead, from the possibility that the input-output map is non-injective.
Compared to standard statistical inference, one could also view non-identifiability as arising because thedata is only partially observed. Non-identifiability appears generally in inverse problems and is not only
restricted to the setting of stochastic simulation. In the latter, however, one additional challenge is that theoutput is represented by a probability distribution. Thus, when speaking of a “match” between the modeland real system behaviors, one needs to think of a close distance between probability distributions.
Our main goal of this paper is to introduce a notion of eligibility set to fundamentally bypass the
non-identifiabliity issue. The idea is that, instead of insisting on consistent parameter estimation which
could be intrinsically impossible as discussed, we relax our “target” to obtaining tight bounds or regions
on all the possible parameter values given the data. Intuitively, when this “eligibility set” is big, this means
the problem is very “ill-posed” and it is difficult to pinpoint exactly the parameter location. On the other
hand, if the set is small, this means we are capable of accurately inferring the true location. Regardlessof the set size, however, the key of this concept is that it is statistically valid in all circumstances. In
other words, the eligibility set contains the true parameter with certain given confidence level, and thus
is correct in a rigorous statistical sense. The set size dictates only the conservativeness level rather thancorrectness. Constructing a valid eligibility set needs to handle both the input-output map and the statisticalnoise from the simulation and real data. We will demonstrate how to do so, and also derive some theoreticalresults regarding our calibration errors. Finally, we validate the performances of our method via numericalexperiments.
We discuss some related works. First, our construction of eligibility set resembles concepts in robust
optimization (RO) (Bertsimas et al. 2011; Ben-Tal et al. 2009), and relatedly distributionally robust
optimization (DRO) (Delage and Ye 2010; Wiesemann et al. 2014). This literature advocates decisionmaking under the worst-case scenario, in situations where the optimization problem contains unknown or
uncertain parameters. The worst-case scenario is typically obtained among a set that likely contains thetrue parameter value, which is often known as the uncertainty set or ambiguity set. In DRO in particular,
the unknown parameter is the underlying probability distribution in a stochastic problem. Our eligibility
set utilizes, in a sense, the concept of uncertainty set that attempts to capture the truth with high likelihood.
Second, our work is related to a couple of recent works on calibration in the simulation literature. Likeus, Goeva et al. (2019) investigate the calibration of input models from output-level data. They focus
on a nonparametric setting in which the input model is represented as a general probability distribution,and employ DRO to calibrate the input model by providing bounds to certain target input-dependent
quantities. However, without an explicit input parameter in their framework, the concept of eligibility
set is not clearly defined there. As a result, their construction of bounds is different from our presented
methodology. Moreover, despite the study of asymptotic effectiveness, they do not provide instructions
on how to choose the simulation size corresponding to the data size as we do here. Other related works
include Lam et al. (2017), Zhang and Zou (2016), Plumlee and Lam (2017), which consider the inference
on the model error between the simulation model and the real world using Bayesian approaches or machine

Bai and Lam
learning methods. Lastly, there is a line of work using maximum entropy to calibrate distributions, inpricing financial derivatives (Avellaneda et al. 2001; Glasserman and Yu 2005), and also in simulation andprobabilistic analysis (Kraan and Bedford 2005; Goeva et al. 2014).
The remainder of this paper is as follows. Section 2 formulates the problem setting and introduces
notations. Sections 3 and 4 describe in detail our concepts and procedure respectively. Section 5 presents
theoretical statistical guarantees. Section 6 demonstrates the numerical results. Section 7 concludes thispaper.
2

PROBLEM SETTING AND NOTATIONS

We consider a family of output probability distributions Pθ where θ denotes the finite-dimensional parameterthat we want to calibrate. We suppose that the true value of θ , denoted by θ0 , is known to lie in theparameter space Θ. We assume that real-world output data, denoted X1 , X2 , · · · , XN ∼ Pθ0 , are available tous. On the other hand, we also have the capability to simulate a random sample Y1θ ,Y2θ , · · · ,Ynθ ∼ Pθ , given
any input value θ ∈ Θ. The latter condition is a basic paradigm in stochastic simulation that is generally
satisfied and essential to our developments. In practice, besides parameter uncertainty, the risk of model
misspecification is also a common issue, but in this paper we focus on solving the former problem assuming
the correctness of the parametric model.
It is noteworthy that we make no assumptions on the relationship between θ and the output distributionPθ . All we need is the simulability of samples from Pθ given θ , and the underlying map can be extremelycomplicated. A generic example would be a general queuing system or network where the inter-arrivaland service time distributions are known to be in certain parametric families. Suppose we have data fromsome output quantities such as the average sojourn times, in which case Pθ0 can denote their distribution.
Given the parameter value, the system dynamics can be readily simulated to generate random samples,
yet it is too complicated to write down analytical expression. In other cases, we could encounter a “blackbox” as the simulation model developed by a third party, in which case our framework would still apply.
To fix some notations, the real output data set X1 , · · · , XN and the simulated output data set Y1θ , · · · ,Ynθ
each determines an empirical probability distribution, which we respectively write as PNθ0 and Pnθ . Moreconcretely, they are given by
1 N
PNθ0 (·) = ∑ δXk (·)
N k=1
and
Pnθ (·) =1 n∑ δY jθ (·)
n j=1

where δXk (·) and δY θ (·) denote the Dirac measures at Xk and Y jθ .
jWe also use F θ to denote the cumulative distribution function of Pθ , so thatF θ0 (x) = P(X1 ≤ x), x ∈ RandF θ (x) = P(Y1θ ≤ x), x ∈ R.
Correspondingly, we denote the empirical distribution functions of X1 , · · · , XN and Y1θ , · · · ,Ynθ respectively
as FNθ0 and Fnθ . That is,
1 N
FNθ0 (x) = ∑ I(Xk ≤ x), x ∈ R
N k=1
and
Fnθ (x) =1 n∑ I(Y jθ ≤ x), x ∈ R
n j=1

Bai and Lam
where I(·) is the indicator function.
3

ELIGIBILITY SET: BASIC CONCEPTS

Our goal is to estimate θ0 from the real output data X1 , · · · , XN and using the capability to simulate samples
given θ , namely Y1θ , · · · ,Ynθ for any given θ of our choice, potentially repeatedly. The standard statistical
approach for this task would be to obtain an estimate as close to θ0 as possible, i.e., a consistent estimatorwith error shrinking to zero as the data size increases.
To understand some challenge of the above and how we can handle it, note that our real observations
only inform us on Pθ0 , and our model allows us to observe Pθ . In other words, any estimation needs tobe done at the output level, by matching these two distributions in some sense. To make it more precise,consider a statistical distance d(·, ·) between two probability distributions. The best scenario one could
hope for is to obtain θ such that d(Pθ , Pθ0 ) = 0. However, because of the intricacy of the input-output map,
there could be other θ 6= θ0 such that the distance is 0. This issue is known as non-identifiability. Moreover,even in the case that there is a unique θ0 giving rise to the zero distance, in a finite sample situation there
could still be many θ that gives small d(Pθ , Pθ0 ), making the problem ill-posed. In scientific computing,these are usually handled via Bayesian approaches or some regularization to identify good candidates ofθ , but since we are looking at the match of distributions in the stochastic simulation setting, it is not clear
these approaches can apply easily.
Our idea is to relax the notion of consistent estimation, or to look for some “best” parameter, into a
region that likely contains the true parameter value. We call this region the eligibility set of θ . To explain,suppose we have an infinite amount of real data, so that Pθ0 is fully known. The eligibility set in this casewould be(1){θ ∈ Θ : d(Pθ , Pθ0 ) = 0}
It is trivial to see that (1) contains θ0 . Note that the choice of statistical distance d(·, ·) does not affect thecorrectness of this claim, but it can affect the conservativeness, i.e., some choices of d(·, ·) can give rise
to a larger-volume set than others. Typically, the smaller the set size, the more successful it is consideredin our calibration task.
To proceed, in practice we only have a finite real data size N, which incurs statistical error in informingPθ0 . In this situation, we considerΘ̃ = {θ ∈ Θ : d(Pθ , PNθ0 ) ≤ η}(2)where η ∈ R+ is a suitable constant. Suppose that {Q : d(Q, PNθ0 ) ≤ η} is a (1 − α)-level confidence region
for the true output distribution Pθ0 , that is,P(d(Pθ0 , PNθ0 ) ≤ η) ≥ 1 − α(3)where P refers to the probability with respect to the real data X1 , . . . , XN . Thus we haveP(θ0 ∈ Θ̃) = P(d(Pθ0 , PNθ0 ) ≤ η) ≥ 1 − α.
In other words, Θ̃ is a 1 − α confidence region for θ0 . We summarize this as:
Theorem 1 Suppose {Q : d(Q, PNθ0 ) ≤ η} is a (1 − α)-level confidence region for the true output distribution
Pθ0 . Then the set Θ̃ in (2) is a (1 − α)-level confidence region for θ0 . Similarly, if the confidence guarantee
for {Q : d(Q, PNθ0 ) ≤ η} holds asymptotically as N increases, then a corresponding asymptotic guarantee
holds for Θ̃.
Note that in general, one can replace {Q : d(Q, PNθ0 ) ≤ η} by any other confidence region of Pθ0 .
However, the form d(Q, PNθ0 ) ≤ η is especially handy since d(Q, PNθ0 ) can represent a goodness-of-fit

Bai and Lam
statistic for Q, and the resulting form of confidence region is obtained via the standard duality from ahypothesis test to a confidence region. Moreover, the form depicted by Θ̃ in (2) has a simple interpretation:
If for some θ ∈ Θ, the distribution Pθ differs significantly from the true distribution Pθ0 , as indicated bya large d(Pθ , PNθ0 ), then we do not accept this choice of θ , and vice versa.
The way we construct Θ̃ above bears resemblance to the RO and DRO literature. This literature considers
decision-making under the worst-case scenario for problems with uncertain or unknown parameters, where
the worst case is taken over a set that likely contains the true parameter value. This latter set is often
known as the uncertainty set or ambiguity set. DRO, in particular, focuses on stochastic problems where
the uncertainty is on some underlying probability distributions. Furthermore, data-driven RO or DRO refers
to situations where the uncertainty set is constructed from the data, typically as a confidence region for the
unknowns. Then, via an optimization problem, one translates the confidence guarantee from the uncertainty
set to confidence bounds on the performance measure or objective of interest. In our discussion, the sets
{Q : d(Q, PNθ0 ) ≤ η} and {θ ∈ Θ : d(Pθ , PNθ0 ) ≤ η} can be viewed as uncertainty sets. Suppose we areinterested in a quantity ψ(Pθ0 ) where ψ is a deterministic function, and then we may evaluate
maximize/minimize ψ(Pθ )subject to d(Pθ , PNθ0 ) ≤ η (i.e., θ ∈ Θ̃).

(4)

Similar to the discussion above, if Θ̃ is a (1 − α)-level confidence region, or equivalently that (3) holds,
then the optimal values of the two optimization problems in (4) provide (1 − α)-level confidence upperand lower bounds to the true value ψ(Pθ0 ). Formulation (4), which aims to compute worst-case bounds
under a confidence region, thus resembles closely the framework of data-driven RO (or DRO). Although
we focus on the calibration of parameter θ in this paper, one may keep in mind that our framework canbe generalized to tackle relevant problems in the form above.
4

ELIGIBILITY SET: CONSTRUCTION PROCEDURE

To construct an eligibility set with the guarantee presented in Theorem 1, we need to resolve two questions.
First, how should we pick a proper discrepancy measure d(·, ·) and the associated η to achieve a given
confidence level 1 − α? Second, how can we determine whether d(Pθ , PNθ0 ) ≤ η or not for any θ ∈ Θ,
especially when the relationship between θ and Pθ is quite complicated?
For the first question, we acknowledge that there are many different ways to measure the discrepancy
between two probability distributions, and also many of them could potentially be used here. Indeed,as long as the confidence property (3) is satisfied for some η > 0, then the pair (d, η) should satisfy
our requirements. However, a good choice d should not only well capture the discrepancy between two
distributions, but also efficient to compute (This is even more so if we were to tackle the optimization
problem (4)). Similarly, η should be easy to calibrate to attain the confidence level 1 − α. Based on thesetwo considerations, a natural choice isd(P1 , P2 ) = sup |F1 (x) − F2 (x)|x∈Rwhere P1 , P2 are probability distributions and F1 , F2 are respectively their cumulative distribution functions.
Then we immediately get thatd(Pθ0 , PNθ0 ) = sup |F θ0 (x) − FNθ0 (x)|,x∈Rwhich is exactly the the Kolmogorov-Smirnov(KS) statistic for the goodness-of-fit for F θ0 . Correspondingly,√
we may choose η = q1−α / N, where, q1−α is the (1 − α)-quantile of supt∈[0,1] |BB(t)| and BB(·) denotes
a standard Brownian bridge. This choice satisfies the asymptotic version of the confidence property (3).

Bai and Lam
Now we address the second question on how to determine whether d(Pθ , PNθ0 ) ≤ η or not. We haveassumed that Pθ comes potentially from a black-box, so that we should resort to sampling Y1θ , · · · ,Ynθ ∼ Pθ
to get Pnθ to approximate Pθ . Correspondingly, we use Fnθ to approximate F θ , and in turn d(Pnθ , PNθ0 ) toapproximate d(Pθ , PNθ0 ). This approximation is utilizable if the simulation size n is much larger than thedata size N, which we will justify in detail in Section 5. In other words, for this approach to work, weneed enough computational capacity measured by the amount of simulation size we can obtain.
To summarize, we now choose Θ̃ defined in (2) more explicitly as


√
θ0
θ
Θ̂ = θ ∈ Θ : sup |Fn (x) − FN (x)| ≤ q1−α / N(5)x∈Rwhere q1−α is the (1 − α)-quantile of supt∈[0,1] |BB(t)|.
We provide a more computationally friendly equivalent to (5), which is quite standard but we show here
for user’s convenience. The empirical distribution functions Fnθ and FNθ0 are both non-decreasing staircase
functions with limx→−∞ Fnθ (x) = limx→−∞ FNθ0 (x) = 0 and limx→∞ Fnθ (x) = limx→∞ FNθ0 (x) = 1. So it is easy
to verify that√sup |Fnθ (x) − FNθ0 (x)| ≤ q1−α / N
x∈R√√⇔|Fnθ (Xk ) − FNθ0 (Xk )| ≤ q1−α / N, |Fnθ (Xk −) − FNθ0 (Xk −)| ≤ q1−α / N, ∀k = 1, · · · , N
whereFnθ (u−) = lim Fnθ (x) =1 n∑ I(Y jθ < u), u ∈ R
n j=1

FNθ0 (u−) = lim FNθ0 (x) =1 N∑ I(Xk < u), u ∈ R.
N k=1

x→u−

and similarly
x→u−Therefore, (5) can be rewritten asno√√Θ̂ = θ ∈ Θ : |Fnθ (Xk ) − FNθ0 (Xk )| ≤ q1−α / N, |Fnθ (Xk −) − FNθ0 (Xk −)| ≤ q1−α / N, ∀k = 1, · · · , N . (6)Consequently, given any θ ∈ Θ, we may simulate Y1θ , · · · ,Ynθ ∼ Pθ and readily check if it is insideΘ̂. If it is, then we say this θ is eligible. However, there are potentially infinitely many θ in Θ to test.
To handle this, we may choose a finite number of values, say θ1 , · · · , θm ∈ Θ, as “representatives” anddetermine their eligibility. The resulting eligible points, collectively denotedon√√Θ̂m = θ1 , . . . , θm ∈ Θ : |Fnθi (Xk ) − FNθ0 (Xk )| ≤ q1−α / N, |Fnθi (Xk −) − FNθ0 (Xk −)| ≤ q1−α / N, ∀k = 1, · · · , N ,provide a discretized approximation to Θ̂. Here, θ1 , · · · , θm can either be deterministic grid points in Θ or
randomly sampled. For instance, when Θ is bounded, we may sample θ1 , · · · , θm ∼ Uni f (Θ). Otherwise,
we may sample θi ’s from some proper (heavy-tailed) distribution on Θ. In the high-dimensional case,θ1 , . . . , θm can be strategically chosen using stochastic gradient-based search with randomly sampled initialpoints.
Putting everything together, we propose Algorithm 1 to compute the eligibility set which can be utilized
to calibrate the unknown parameter θ .

Bai and Lam
Algorithm 1 Constructing eligibility set of θ .
Require: The output data X1 , · · · , XN . The number of candidate θ ’s m. The simulation replication size n.
The confidence level 1 − α.
Ensure: An eligibility set Θ̂m .
1: Generate θ1 , · · · , θm ∈ Θ;
2: for i = 1, · · · , m do
3:
Generate a random sample Y1θi , · · · ,Ynθi ∼ Pθi ;4:for k = 1, · · · , N do
5:Compute ai,k = n1 ∑nj=1 I Y jθi ≤ Xk − N1 ∑Nl=1 I(Xl ≤ Xk ) ;6:Compute bi,k = n1 ∑nj=1 I Y jθi < Xk − N1 ∑Nl=1 I(Xl < Xk ) ;7:
end for
8: end for√√
9: Construct Θ̂m = {θi : ai,k ≤ q1−α / N and bi,k ≤ q1−α / N, ∀k = 1, · · · , N};
10: return The eligibility set Θ̂m .
5

THEORETICAL STATISTICAL GUARANTEESAs discussed in Section 4, our choice of (d, η) satisfies the asymptotic version of (3), that is,√lim P sup |F θ0 (x) − FNθ0 (x)| ≤ q1−α / N ≥ 1 − α.
N→∞

x∈R

However, in Algorithm 1, we substitute F θ0 with Fnθ0 , and thus additional randomness is introduced.
To justify our approach, we need to show that a similar asymptotic confidence property holds with thissubstitution. Theorem 2 demonstrates a sufficient condition for this.
Theorem 2 Suppose that X1 , · · · , XN is an i.i.d. random sample from Pθ0 and Y1θ0 , · · · ,Ynθ0 is another i.i.d.
random sample from Pθ0 . FNθ0 and Fnθ0 are respectively the empirical distribution functions of the two
random samples. F θ0 denotes the cumulative distribution function of Pθ0 . If n = ω(N) as N → ∞, then


√
θ0
θ0
lim P sup |Fn (x) − FN (x)| ≤ q1−α / N ≥ 1 − α.
n,N→∞Proof.

x∈R

Since
sup |Fnθ0 (x) − FNθ0 (x)| ≤ sup
x∈R
x∈R



θ0
θ0
θ0
θ0
|Fn (x) − F (x)| + |F (x) − FN (x)|≤ sup |Fnθ0 (x) − F θ0 (x)| + sup |FNθ0 (x) − F θ0 (x)|,
x∈R

x∈R

we get that


√
θ0
θ0
P sup |Fn (x) − FN (x)| > q1−α / N
x∈R


√
θ0
θ0
θ0
θ0
≤P sup |Fn (x) − F (x)| + sup |FN (x) − F (x)| > q1−α / N
x∈R
x∈R




√
√
θ0
θ0
θ0
θ0
≤P sup |Fn (x) − F (x)| > λ q1−α / N + P sup |FN (x) − F (x)| > (1 − λ )q1−α / N
x∈R

x∈R

Bai and Lam√
for any λ ∈ (0, 1). It is known that n supx∈R |Fnθ0 (x) − F θ0 (x)| ⇒ supt∈R |BB(F θ0 (t))| and similarly√
N supx∈R |FNθ0 (x) − F θ0 (x)| ⇒ supt∈R |BB(F θ0 (t))| as n, N → ∞ where ⇒ stands for convergence in distribution. n = ω(N) as N → ∞ implies that n/N → ∞ as N → ∞, and thus




p
√
√
θ0
θ0
θ0
θ0
n sup |Fn (x) − F (x)| > λ q1−α n/N → 0
P sup |Fn (x) − F (x)| > λ q1−α / N = P
x∈R

x∈R

as N → ∞ for any λ ∈ (0, 1). Hence,




√
√
θ0
θ0
θ0
θ0
lim P sup |Fn (x) − FN (x)| > q1−α / N ≤ lim P sup |FN (x) − F (x)| > (1 − λ )q1−α / N
N→∞
N→∞
x∈R
x∈R


= P sup |BB(F θ0 (t))| > (1 − λ )q1−α
t∈R
!
≤P

sup |BB(t)| > (1 − λ )q1−α.

t∈[0,1

By the arbitrariness of λ and the definition of q1−α , we finally get that√
lim P sup |Fnθ0 (x) − FNθ0 (x)| > q1−α / N ≤ P

N→∞

x∈R

!
sup |BB(t)| ≥ q1−α= α.

t∈[0,1]

Note that the probability of interest in Theorem 2 relates closely to the type I error in statistical inference.
According to the theorem, conditional on being chosen as a “representative”, θ0 is correctly selected to
put inside the eligibility set with probability at least 1 − α asymptotically. In this sense, we may regard thediscretized eligibility set Θ̂m as (upon a suitable interpolation) an approximate confidence region for θ0 .
On the other hand, we are also interested in the type II error, that is, the probability that some
representative does not belong to the true (discretized) eligibility set yet is accepted. The type II error
here governs the conservativeness level of Θ̂m as an approximate confidence region. Note that, in the caseF θ 6= F θ0 and as n and N go to infinity, the difference between Fnθ and FNθ0 gradually converges to thedifference between F θ and F θ0 , which is a positive constant. Theorem 3 captures the rate under which
this occurs and provides an upper bound for the probability that F θ 6= F θ0 yet θ is eligible.
Theorem 3 Suppose that X1 , · · · , XN is an i.i.d. random sample from Pθ0 and Y1θ , · · · ,Ynθ is an i.i.d.
random sample from Pθ . FNθ0 and Fnθ are respectively the empirical distribution functions of the two
random samples. F θ0 and F θ denote the cumulative distribution functions of Pθ0 and Pθ . Suppose thatsupx∈R |F θ (x) − F θ0 (x)| > 0. For any ε1 , ε2 > 0 such that ε1 + ε2 < supx∈R |F θ (x) − F θ0 (x)|, if

N>
then

q1−α
θ
supx∈R |F (x) − F θ0 (x)| − ε1 − ε2

2
,





√
2
2
θ0
θ
P sup |Fn (x) − FN (x)| ≤ q1−α / N ≤ 2 e−2nε1 + e−2Nε2 .
x∈R

Proof.

We know that
|Fnθ (x) − FNθ0 (x)| ≥ |F θ (x) − F θ0 (x)| − |Fnθ (x) − F θ (x)| − |FNθ0 (x) − F θ0 (x)|,Bai and Lam
and then we get thatsup |Fnθ (x) − FNθ0 (x)| ≥ sup |F θ (x) − F θ0 (x)| − sup |Fnθ (x) − F θ (x)| − sup |FNθ0 (x) − F θ0 (x)|.
x∈R

x∈R

x∈R

x∈R

Therefore,


√
θ0
θ
P sup |Fn (x) − FN (x)| ≤ q1−α / N
x∈R


√
≤P sup |F θ (x) − F θ0 (x)| − sup |Fnθ (x) − F θ (x)| − sup |FNθ0 (x) − F θ0 (x)| ≤ q1−α / N
x∈R
x∈R
x∈R




θ0
θ
θ
θ0
≤P sup |Fn (x) − F (x)| > ε1 + P sup |FN (x) − F (x)| > ε2 .
x∈R

x∈R

√
The last inequality is obtained using the fact that q1−α / N < supx∈R |F θ (x) − F θ0 (x)| − ε1 − ε2 under theconditions in the theorem. By the refined Dvoretzky–Kiefer–Wolfowitz (DKW) inequality (Kosorok 2008),
we get that




2
θ0
θ
θ
−2nε12
θ0
P sup |Fn (x) − F (x)| > ε1 ≤ 2e
, P sup |FN (x) − F (x)| > ε2 ≤ 2e−2Nε2 ,
x∈R

x∈R

which concludes the proof.
With Theorem 3, we may easily develop a type II error guarantee for Algorithm 1. More specifically,
we can bound the probability that Θ̂m contains some θ ∈ Θ such that supx∈R |F θ (x) − F θ0 (x)| is large.
Theorem 4 shows the details.
Theorem 4 We follow Algorithm 1 to obtain Θ̂m . For any ε, ε1 , ε2 > 0 such that ε1 + ε2 < ε, if

2
q1−α
N>
,
ε − ε1 − ε2
then





2
2
P ∃i = 1, · · · , m s.t. sup |F θi (x) − F θ0 (x)| > ε, θi ∈ Θ̂m ≤ 2m e−2nε1 + e−2Nε2 .
x∈R

Proof.

By applying Theorem 3, we have that


θi
θ0
P ∃i = 1, · · · , m s.t. sup |F (x) − F (x)| > ε, θi ∈ Θ̂m
x∈R


m
θi
θ0
≤ ∑ P sup |F (x) − F (x)| > ε, θi ∈ Θ̂m
i=1
m

x∈R

i=1

x∈R



√
θ0
θi
θ0
θi
≤ ∑ P sup |Fn (x) − FN (x)| ≤ q1−α / N sup |F (x) − F (x)| > ε
x∈R



2
2
≤2m e−2nε1 + e−2Nε2 .

Theorems 2 and 4 together justify that the eligibility set Θ̂m given by Algorithm 1 can help us calibrate theunknown parameter θ , especially when m, n, N are sufficiently large and chosen properly. More specifically,
the following corollaries provide instructions on how to choose the sizes.

Bai and Lam√
Figure 1: Apply Algorithm 1 to calibrate µ ( N supx∈R |Fnθi (x) − FNθ0 (x)| against θi ).
Corollary 5 We follow Algorithm 1 to obtain Θ̂m . If log m = o(N) and n = Ω(N) as N → ∞, then for any
ε > 0,

lim P ∃i = 1, · · · , m s.t. sup |F θi (x) − F θ0 (x)| > ε, θi ∈ Θ̂mm,n,N→∞= 0.

x∈R

Corollary 6 We follow Algorithm 1 to obtain Θ̂m . If m = o(N) and n = Ω(N) as N → ∞, then


p
θi
θ0
lim P ∃i = 1, · · · , m s.t. sup |F (x) − F (x)| > log m/m, θi ∈ Θ̂m = 0.
m,n,N→∞

6x∈R

NUMERICAL RESULTS

We conduct some numerical experiments to validate our approach. We consider a G/G/1 queuing system
where the inter-arrival time is known to follow a Gamma(λ , 1) distribution while the service time follows a
Weibull(µ, 1) distribution. Then the expected inter-arrival time and the expected service time are respectively
λ and µ. We assume that we know λ , µ ∈ (0, 2) and that the true values are λ0 = 1/2 and µ0 = 1. Supposethat we are only given a random sample of the average sojourn time of the first c = 10 customers in thesystem. Clearly the distribution of this output quantity is uniquely determined by λ and µ, and conversely,
given the values of λ , µ, the average sojourn time is easy to simulate. Thus this parameter calibration
problem falls into our settings.
We first consider the case where λ is known and we want to calibrate µ (i.e., θ := µ). Suppose
we have N = 100 output observations generated from the unknown µ0 = 1, which we mimic by runningsimulation. We follow Algorithm 1 by setting m = n = 1000 and sampling θ1 , · · · , θm ∼ Uni f (0, 2). For
each θi , we simulate n = 1000 versions of average sojourn time with λ = 1/2 and√ µ = θi , and then computesupx∈R |Fnθi − FNθ0 |. We choose 1 − α = 0.95. By comparing the realized value, N supx∈R |Fnθi − FNθ0 |, with√
q1−α , we determine whether θi is eligible. Figure 1 plots N supx∈R |Fnθi (x) − FNθ0 (x)| against θi , and thered horizontal line corresponds to q1−α ≈ 1.36. Thus the dots below the red line constitute the eligibility
set Θ̂m . From the graph, we find that all the eligible θi ’s are relatively close to the true value θ0 = 1. Morespecifically, the smallest and the largest values in the eligibility set are respectively 0.90 and 1.15, which
has greatly shrunk the original a priori known range (0, 2).
Additionally, in the above setting, we may construct an approximate confidence interval with thesmallest and the largest values in Θ̂m . After repeating the above process for L = 100 times, we find thatamong the confidence intervals that we obtain, 99 cover the true value 1. Given that the target confidence
level is 95%, the coverage probability is satisfactory, though a little conservative.
Now we consider the case that both λ and µ are unknown and we calibrate them simultaneously (i.e.,θ := (λ , µ)). This problem is non-identifiable. It is easy to verify that, as λ and µ grow simultaneously
in a certain way, the distribution of the average sojourn time of the first c = 10 customers does notchange. Similar to the previous setup, given an output data set of size N = 100, we run Algorithm 1 withm = n = 1000, 1 − α = 0.95, but now we generate θ1 , · · · , θm ∼ Uni f (0, 2)2 . In Figure 2, we show theBai and Lam(a)√
N supx∈R |Fnθi (x) − FNθ0 (x)| against θi .

(b) The eligibility set Θ̂m .

Figure 2: Apply Algorithm 1 to calibrate (λ , µ).
√three-dimensional scatter plot of N|Fnθi (x) − FNθ0 (x)| against θi = (λi , µi ) as well as the two-dimensional
visualization of the eligibility set Θ̂m . Particularly, in Figure 2(b), the red cross mark denotes the true value(1/2, 1), which lies approximately in the band. We therefore see that, when faced with non-identifiability,
our approach can still result in a region that contains the truth, where the size of the region is related tothe degree of non-identifiability (i.e., the range of parameters that exactly match in terms of the output
distribution). This demonstrates how our approach can bypasss non-identifiability and provides empiricalsupport to our proposed framework.
7

CONCLUSION

In this paper, we have studied a framework to calibrate unknown parameter only using the output-level
data, which only makes very mild assumptions on the model. Indeed, theoretically our framework may
deal with the situation where the relationship between the input parameter and the output distribution iscomplicated or even a black box.
In particular, we consider a finite number of values of the parameter as “representatives” and then foreach of them, we generate a simulated sample, which may be compared to the real data using the KS
statistic. The eligibility set consists of the “representatives” whose simulated sample is “close enough” tothe true data sample, and then we employ this set to estimate the parameter of interest.
We have also developed statistical guarantees for this procedure. To be more specific, we analyze thetype I error as well as the type II error of the eligibility set in terms of the data size and the simulation size.
These theoretical results may provide some guidance on how to choose the simulation size according to thedata size. Moreover, the numerical results support that our framework works well even under complicatedsettings.
ACKNOWLEDGMENTS
We gratefully acknowledge support from the National Science Foundation under grants CAREER CMMI1834710 and IIS-1849280.
REFERENCES
Avellaneda, M., R. Buff, C. Friedman, N. Grandechamp, L. Kruk, and J. Newman. 2001. “Weighted Monte Carlo: a New
Technique for Calibrating Asset-Pricing Models”. International Journal of Theoretical and Applied Finance 4(01):91–119.
Balci, O., and R. G. Sargent. 1982. “Some Examples of Simulation Model Validation Using Hypothesis Testing”. In Proceedingsof the 14th Conference on Winter Simulation - Volume 2, WSC ’82, 621–629: Winter Simulation Conference.
Barton, R. R. 2012. “Tutorial: Input Uncertainty in Output Analysis”. In Proceedings of the 2012 Winter Simulation Conference(WSC), 1–12. IEEE.
Ben-Tal, A., L. El Ghaoui, and A. Nemirovski. 2009. Robust Optimization, Volume 28. Princeton University Press.

Bai and Lam
Bertsimas, D., D. B. Brown, and C. Caramanis. 2011. “Theory and Applications of Robust Optimization”. SIAM review 53(3):464–
501.
Chick, S. E. 2006. “Bayesian Ideas and Discrete Event Simulation: Why, What and How”. In Proceedings of the 2006 wintersimulation conference, 96–106. IEEE.
Currin, C., T. Mitchell, M. Morris, and D. Ylvisaker. 1991. “Bayesian Prediction of Deterministic Functions, with Applications
to the Design and Analysis of Computer Experiments”. Journal of the American Statistical Association 86(416):953–963.
Delage, E., and Y. Ye. 2010. “Distributionally Robust Optimization under Moment Uncertainty with Application to Data-Driven
Problems”. Operations research 58(3):595–612.
Glasserman, P., and B. Yu. 2005. “Large Sample Properties of Weighted Monte Carlo Estimators”. Operations Research 53(2):298–
312.
Goeva, A., H. Lam, H. Qian, and B. Zhang. 2019, aug. “Optimization-Based Calibration of Simulation Input Models”. Operations
Research 67(5):1362–1382.
Goeva, A., H. Lam, and B. Zhang. 2014. “Reconstructing Input Models via Simulation Optimization”. In Proceedings of theWinter Simulation Conference 2014, 698–709. IEEE.
Henderson 2003. “Input Model Uncertainty: Why Do We Care and What Should We Do about It?”. In Proceedings of the2003 Winter Simulation Conference, 2003., Volume 1, 90–100 Vol.1.
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.
Kosorok, M. R. 2008. Introduction to Empirical Processes and Semiparametric Inference. New York: Springer.
Kraan, B., and T. Bedford. 2005. “Probabilistic Inversion of Expert Judgments in the Quantification of Model Uncertainty”.
Management science 51(6):995–1006.
Lam, H. 2016. “Advanced Tutorial: Input Uncertainty and Robust Analysis in Stochastic Simulation”. In 2016 Winter Simulation
Conference (WSC), 178–192. IEEE.
Lam, H., X. Zhang, and M. Plumlee. 2017. “Improving Prediction from Stochastic Simulation via Model Discrepancy Learning”.
In 2017 Winter Simulation Conference (WSC), 1808–1819. IEEE.
Nelson, B. L. 2016, feb. “’Some Tactical Problems in Digital Simulation’ for the Next 10 Years”. Journal of Simulation 10(1):2–11.
Plumlee, M., and H. Lam. 2017. “An Uncertainty Quantification Method for Inexact Simulation Models”. arXiv preprint
arXiv:1707.06544.
Sargent, R. G. 2010. “Verification and Validation of Simulation Models”. In Proceedings of the 2010 winter simulation
conference, 166–183. IEEE.
Schruben, L. W. 1980, mar. “Establishing the Credibility of Simulations”. Simulation 34(3):101–105.
Song, E., B. L. Nelson, and C. D. Pegden. 2014. “Advanced Tutorial: Input Uncertainty Quantification”. In Proceedings of theWinter Simulation Conference 2014, 162–176. IEEE.
Tarantola, A. 2005. Inverse Problem Theory and Methods for Model Parameter Estimation, Volume 89. siam.
Wiesemann, W., D. Kuhn, and M. Sim. 2014. “Distributionally Robust Convex Optimization”. Operations Research 62(6):1358–
1376.
Zhang, X., and L. Zou. 2016. “Simulation Metamodeling in the Presence of Model Inadequacy”. In 2016 Winter Simulation
Conference (WSC), 566–577. IEEE.

AUTHOR BIOGRAPHIESYUANLU BAI is a PhD student in the Department of Industrial Engineering and Operations Research at Columbia University.
She received an M.S. in Operations Research from Columbia University, and a B.S. in Statistics as well as a B. Ec in Economics
from Peking University. Her research interest lies in stochastic simulation, such as uncertainty quantification and extreme risk
analysis. Her email address is yb2436@columbia.edu.
HENRY LAM is an Associate Professor in the Department of Industrial Engineering and Operations Research at Columbia
University. He received his Ph.D. degree in statistics from Harvard University in 2011, and was on the faculty of Boston
University and the University of Michigan before joining Columbia in 2017. His research focuses on Monte Carlo simulation,
uncertainty quantification, risk analysis, and stochastic and robust optimization. His email address is henry.lam@columbia.edu.