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

UNCERTAINTY QUANTIFICATION ON SIMULATION ANALYSIS DRIVEN BY RANDOMFORESTS
Amirhossein MeisamiMark P. Van Oyen

Henry Lam

Department of Industrial and Operations
Engineering
University of Michigan
1205 Beal Ave
Ann Arbor, MI 48109, USA

Department of Industrial Engineering andOperations Research
Columbia University
500 W. 120th Street
New York, NY 1002, USA

ABSTRACT
We consider simulation-based estimation where the inputs are calibrated from predictive models generated
by random forests (RFs). RF is a common technique to produce ensemble predictions by aggregating many
individual decision trees. This problem arises in data-driven applications such as individualized surgery
operations scheduling and supply chain management. We investigate the estimation of the output variance
contributed from the noises in the input prediction, which can be used to construct output confidence
intervals. In particular, we study the integration of simulation runs with a recently proposed infinitesimal
jackknife estimator that can reduce the computational burden from double layers of bootstrapping. We
demonstrate our scheme with an elementary numerical example.
1

INTRODUCTION

Modern simulation analysis has been increasingly linked with learning-based predictive modeling; see,e.g., the recent literature in simulation analytics (Jiang, Hong, and Nelson 2016). This paper presents a
study along this line, and considers simulation analysis driven by input distributions that are generated frompredictive models. More precisely, consider the estimation of ψ(F1 , . . . , Fm ) where ψ(·) is an output summaryof interest that can be simulated from the input distributions F1 , . . . , Fm . Here, Fi ’s are distributions given
individual sets of covariates that are calibrated by learning from historical data. Our interest is uncertainty
quantification schemes for the simulation output ψ(F1 , . . . , Fm ) under this setting, in particular the estimationof the variance stemming from the noises in building the predictive model, and subsequently the constructionof an output confidence interval.
This particular problem arises from the authors’ own encounter in individualized surgery operations
scheduling (Meisami et al. 2017). In this context, ψ(F1 , . . . , Fm ) is the expected daily total waiting time ina surgery room. F1 , . . . , Fm denote the distributions of the surgery times for m scheduled patients. These
distributions come from a predictive model given each patient’s characteristics, such as age, gender, past
medical history etc. Clearly, an informed decision-making based on the performance measure requires
proper quantification of the variability due to the errors coming from the input prediction.
Another natural instance of this problem, in the area of supply chain and inventory management, isthe common newsvendor problem. For instance, in an m-product multi-stage newsvendor problem, thedecision maker has to choose the order quantity of each product at different times based on the predicted
demand such that cumulative back-order and inventory cost is minimized. The decision maker often has
access to information such as customer demographics, seasonality, economic conditions, etc. and use suchdata to predict the product demand distributions F1 , . . . , Fm . The quantity ψ(F1 , . . . , Fm ) here refers to theMeisami, Lam, and Van Oyen
expected cumulative cost under a given ordering policy. Much like the surgery scheduling scenario, it isoperationally important to capture the variation and uncertainty in deciding a desirable policy.
One common method in building predictive models for Fi ’s with covariates, which we focus on in thispaper, is the quantile regression forest (QRF) (Meinshausen 2006). This is a generalization of random
forest (RF) that produces ensemble prediction by aggregating many decision trees together, where each tree
is grown by random splitting and a bootstrap replication of the data set (known as bootstrap aggregating,or bagging in short; Breiman 1996, Breiman 2001). RF and bagging have been shown, both theoreticallyand empirically, to have desirable properties (e.g., variance reduction, stability; Büchlmann and Yu 2002,
Efron 2014), and have become one of the most popular off-the-shelf learning methods. QRF differs fromthe standard RF in that the output of each tree is a probability distribution on the response variable, sothat their aggregation give rise to an overall predicted distribution. Along a similar vein as RF, QRF has
been shown to demonstrate superior predictive power than standard tools like quantile regression in certain
scenarios (Meinshausen 2006).
Our current work on uncertainty quantification, through estimating the variance contributed from input
noises and its subsequent use in confidence interval construction, follows the stream of literature on input
uncertainty (e.g., Barton (2012), Song, Nelson, and Pegden (2014), Lam (2016)). In particular, we utilize
the general technique of the bootstrap and the delta method. Cheng and Holland (1997) considers theestimation of input-induced variance via both the delta method and the bootstrap, the latter used together
with an analysis-of-variance adjustment to reduce bias. Cheng and Holland (1998) and Cheng and Holland(2004) suggest a more efficient refinement of the delta method called the two-point method. Song and
Nelson (2015) studies the bootstrap on a mean-variance model in capturing input uncertainty.
To estimate the sampling variance from QRF (or RF) for a prediction problem, standard bootstrapping
would require running two layers of resampling, one on the original data and one on the resampled data
used to construct the QRFs. In simulation driven by these QRFs, there is an additional effort of simulation
runs, thus posing enormous computational burden. As the main focus of this paper, we investigate theintegration of simulation with an infinitesimal jackknife (IJ) estimator, proposed by the recent work Efron(2014) and Wager, Hastie, and Efron (2014) that avoids double-layer bootstrapping for RFs in prediction
problems, in estimating the output variance contributed from the input noises of the QRF. Throughout thispaper we will focus on the essential ideas and omit the rigorous mathematical developments.
In the rest of this paper, Section 2 will first describe the basics of QRF and how to use it in simulation
analysis. Section 3 then introduces the IJ estimator and discusses how to blend it efficiently in output
uncertainty quantification. Section 4 shows a numerical example.
2

SIMULATION ANALYSIS DRIVEN BY QUANTILE REGRESSION FORESTS

We first describe some basics of QRF, followed by a discussion on how it is used for input modeling insimulation analysis.
2.1 Basics of Quantile Regression Forest
We follow the discussion in Meinshausen (2006). Given a collection of i.i.d. data (X j ,Y j ), j = 1, . . . , n,where X j ∈ B for some appropriate space B is the covariate or feature vector, and Y j ∈ R is the response,
the goal is to estimate F(y|X = x) = P(Y ≤ y|X = x) for all y ∈ R, the conditional probability distributionof the response Y given the feature X = x.
To start, we first grow a tree by splitting the feature space B into non-overlapping rectangular subspaces,
say Rl , l = 1, . . . , L, each representing a leaf. This tree is generated randomly, by using a bootstrap resample{(Xi∗ ,Yi∗ )}i=1,...,n , i.e., a sample of size n drawn from {(X j ,Y j )} j=1,...,n under sampling-with-replacement,
and also randomly splitting the feature space according to some independent random source denoted by
θ . For each feature value x, we denote Rl(x,θ ) as the leaf that x lies on. For convenience, we denote(X ∗ ,Y ∗ ) = ((X j∗ ,Y j∗ ) : j = 1, . . . , n) as the bootstrapped data set. The output of a tree, given x, is theMeisami, Lam, and Van Oyenempirical distribution of Y ∗ inside the same leaf, i.e.,t(y; x, (X ∗ ,Y ∗ ), θ ) =∑ j:X j∗ ∈Rl(x,θ ) I(Y j∗ ≤ y)#{ j : X j∗ ∈ Rl(x,θ ) }We grow these random trees many times, say B times, where each tree contains the bootstrap sample(X b ,Y b ) = ((X jb ,Y jb ) : j = 1, . . . , n) and random realization θ b . Then the QRF outputs the average of alltrees, i.e., given x, it outputs
1 B(1)F̂(y|X = x) = ∑ t(y; x, (X b ,Y b ), θ b )B b=1
Besides its superior empirical performance (e.g., compared to quantile regression, especially for extremal
quantiles), QRF has been shown to be consistent, nonparametrically, when the leave sizes of the trees(measured by the proportions of data) are grown at a suitable rate relative to the sample size and thesplitting mechanism are properly taken (Meinshausen 2006). However, convergence rate results appear
only available in the case of subsampled trees (i.e., the bootstrap resample size grows at a smaller rate thann); see, e.g., Mentch and Hooker (2016), Wager and Athey (2017).
2.2 Input ModelingConsider the performance measure described in the introduction ψ(F1 , . . . , Fm ), where given Fi ’s, ψ(·) canbe estimated by running simulation. An example would be the expectation-type measureψ(F1 , . . . , Fm ) = EF1 ,...,Fm [h(Y1 , . . . , Ym )](2)

where Yi denotes a sequence of i.i.d. variates generated from Fi , and h is some cost function that can beevaluated given the Yi ’s.
We consider the situation where each input distribution Fi is a conditional distribution given a feature
value xi , i.e., Fi (y) = P(Y ≤ y|X = xi ). As mentioned in the introduction, these Fi ’s could represent the
distributions of the random surgery times given patients’ observed characteristics. Using the QRF, a natural
estimate for ψ(F1 , . . . , Fm ) will be ψ̂(F̂1 , . . . , F̂m ), where F̂i (y) = F̂(y|X = xi ) is given by (1), and ψ̂ is asimulation approximation of ψ. In the context of (2), this meansψ̂(F̂1 , . . . , F̂m ) =1 R∑ h(Yr1 , . . . , Yrm )R r=1

where each Yri consists of independent realizations of variates generated from input distribution F̂i . Under
the assumptions that guarantee nonparametric consistency of QRF and mild regularity conditions on ψ(·),ψ̂(F̂1 , . . . , F̂m ) is consistent for estimating ψ(F1 , . . . , Fm ) when the bootstrap size B in the QRF and thesimulation size R grow to ∞, and hence ψ̂(F̂1 , . . . , F̂m ) is a justified estimate of the performance measure.
For reliable decision-making, we need a measure of the variability on the point estimate of ψ(F1 , . . . , Fm ).
This variability comes from both the noises of the simulation replications (called stochastic or simulation
uncertainty) and the input data (called input uncertainty); see, e.g., Barton (2012), Song, Nelson, and Pegden(2014), Lam (2016). This variability can be captured by a confidence interval (CI), with a confidence levelsay 1 − α, given by
"
#
r
r
2
2
τ
τ
ψ̂(F̂1 , . . . , F̂m ) − z1−α/2 σ 2 + , ψ̂(F̂1 , . . . , F̂m ) + z1−α/2 σ 2 +(3)RRwhere z1−α/2 is the 1 − α/2 quantile of the standard normal distribution. The variance component σ 2 is the
input-induced variance, whereas τ 2 /R is the simulation-induced variance. Interval (3) has assumed that theMeisami, Lam, and Van Oyenbootstrap size B is large enough to make the resampling variance negligible, and a Gaussian approximation
is asymptotically valid. This assumption is not formally justified here. In Section 3, we will consider amodified version of (3) and provide further discussion.
The quantity τ 2 is relatively easy to estimate, for instance by taking the sample variance of all thesimulation replications. The quantity σ 2 , in the nonparametric setting, typically involves estimating thevariance of the so-called influence function of ψ(·), which can be viewed as a functional derivative withrespect to the input distributions (Hampel et al. 2011, Huber 2011). It could be difficult to efficiently
estimate σ 2 by directly using the form of the influence function (e.g., in a long-horizon problem). Therefore,bootstrapping has been used as a major tool for estimating σ 2 (e.g., Cheng and Holland 1997, Song and Nelson2015). In the current context, this is done by repeatedly drawing bootstrap resamples from {(Xi ,Yi )}i=1,...,n .
For each resample, fit the QRF pretending the resample is the original data, and calculate ψ̂(F̂1 , . . . , F̂m ).
The variance of these outcomes gives the bootstrap variance that approximates σ 2 .
Note that the above scheme for estimating σ 2 requires two layers of bootstrapping, the first corresponding
to the generation of resamples on which QRFs are trained on, and the second to the generation of trees inconstructing the QRFs and subsequently the drive of simulation runs. The overall computational effort ison the order of the product of the outer bootstrap size and the sum of the tree and the simulation replication
sizes.
3

SIMULATION-BASED INFINITESIMAL JACKKNIFE ESTIMATOR

Given the difficulties in approximating σ 2 in (3), we consider an alternate approach that requires looking at
a slightly different point estimator. For convenience, let F = (F1 , . . . , Fm ) and F∗ (y) = (t(y; xi , (X ∗ ,Y ∗ ), θ ∗ ) :i = 1, . . . , m) as the vector of input distributions predicted by each tree. Our point estimator depicted
in Section 2.2 can be expressed as ψ̂(Ê[F∗ ]) where Ê[·] denotes the sample average of B independent
realizations of trees. As the numbers of trees and simulation replications grow to ∞, this approximatesψ(E∗ [F∗ ]) where E∗ [·] denotes the expectation generated by the bootstrapped trees. CI (3) is derived fromthe conjectured central limit theorem ψ(E∗ [F∗ ]) ∼ N(ψ(F), σ 2 ).
Our considered alternate point estimator takes the formÊ[ψ̂(F∗ )](4)that interchanges Ê and ψ̂. As the numbers of trees and simulation replications grow to ∞, the estimator(4) approximatesE∗ [ψ(F∗ )](5)
where E∗ [·] denotes the expectation with respect to the bootstrap. We note that (5) is equivalent toE∗ [ψ ∗ (F∗ )],(6)

where ψ ∗ (F∗ ) denotes an unbiased simulation replication in estimating ψ(F∗ ) given F∗ (e.g., h(Y1 , . . . , Ym )with Yi being i.i.d. sequence generated from Fi∗ ), and with slight abuse of notation, E∗ [·] in (6) denotes
the expectation generated from the combination of a bootstrapped tree and the simulation replication. Notethat (6) can be viewed as the limiting output of a bagging scheme, where each “base learner” is now arealization of ψ ∗ (F∗ ). From this view, an implementable bagging scheme to approximate (6) isÊ[ψ ∗ (F∗ )] =1 B
∑ h(Yb1 , . . . , Ybm )
B b=1

(7)

where (Ybi )i=1,...,m is generated from Fb = (t(y; xi , (X b ,Y b ), θ b ) : i = 1, . . . , m). This is similar to using
the established QRF, but instead of outputting t(y; x, (X b ,Y b ), θ b ) for each tree, we output a copy of thesimulation output generated by the input models Fb .

Meisami, Lam, and Van Oyen
Our subsequent analysis assumes that when n is large and under smoothness conditions on ψ(·),
F∗ ≈ E∗ [F∗ ] and E∗ [ψ ∗ (F∗ )] = E∗ [ψ(F∗ )] ≈ ψ(E∗ [F∗ ]) well enough so that a central limit theorem holds
approx.
for E∗ [ψ ∗ (F∗ )] and implies E∗ [ψ ∗ (F∗ )] ∼ N(ψ(F), σ 2 ). Thus, together with the simulation variability,a (1 − α)-level CI for ψ(F) is
#
"
r
r
τ2
τ2
∗ ∗
∗ ∗
2
2
(8)
Ê[ψ (F )] − z1−α/2 σ + , Ê[ψ (F )] + z1−α/2 σ +RRThe main purpose of using Ê[ψ ∗ (F∗ )] instead of ψ̂(F̂) in Section 2.2 is the availability of an estimate forσ 2 without resorting to a double-layer bootstrapping, given by the infinitesimal jackknife (IJ) estimator
proposed by Efron (2014) and Wager, Hastie, and Efron (2014), in the formn

V̂IJ = ∑ Cov∗ (Ni∗ , ψ ∗ (F∗ ))2

(9)

i=1

where Ni∗ is the number of counts of the observation (Xi ,Yi ) appearing in the bootstrap resample for building
F ∗ , and Cov∗ (·, ·) denotes the covariance generated by the boostrapping. Under some conditions usingpsubsampled trees, it is shown that V̂IJ /σ 2 → 1 (Wager and Athey 2017). However, as far as we know, there
have been no such convergence results available for QRFs with full-size trees, and we may need to resort
to the use of subsampled trees for a rigorous analysis leading to (8).
A bias-reduced estimator for (9) using a finite tree size B (Wager, Hastie, and Efron 2014) is given byBnd 2i − n ∑ (ψ b (Fb ) − ψ ∗ (F∗ ))2
V̂IJB = ∑ Cov
B2 b=1
i=1
where

(10)

B

d i = 1 ∑ (Nib − 1)(ψ b (Fb ) − ψ ∗ (F∗ ))
Cov
B b=1
Nib and ψ b (Fb ) denote the realizations of Ni∗ and ψ ∗ (F∗ ) in the b-th tree, and ψ ∗ (F∗ ) = (1/B) ∑Bb=1 ψ b (Fb ).
As explained in Wager, Hastie, and Efron (2014), the first term in (10) is the plug-in estimator of (9), andthe second term corrects the bias of the plug-in estimator from order σ 2 /(nB) to order nσ 2 /B. We notethat Wager, Hastie, and Efron (2014) focuses on the mean prediction, but their analysis clearly applies toany type of function on the response, including ψ ∗ (F∗ ).
We have used one intra-tree simulation replication in the illustration above. Note that (6) remains
unchanged if one increases the intra-tree simulation replication size from one to any number ((5) is anextreme that uses an infinite number of simulation replications). The conjectured asymptotic analysis
remains the same except that the intra-tree variance for ψ ∗ (F∗ ) decreases in the usual Monte Carlo rate.
Thus, one can use V̂IJB with ψ b (·) redefined accordingly. Our numerical experiment in the next section
indicates that under a fixed budget, a single intra-tree replication does appear to perform superior to otherchoices.
4

NUMERICAL EXAMPLE

In this section we provide numerical results on the method discussed in Section 3 using the well-known“Auto MPG” data set that was analyzed in Wager, Hastie, and Efron (2014). We use 7 features (numberof cylinders, displacement, horsepower, weight, acceleration, year, and origin) over a training set of 312
observations to build our QRF. For this proof-of-concept example, we fix four feature values in a separate
test set, and consider a performance measureψ(F1 , . . . , F4 ) = EF1 ,...,F4 [h(Y1 , . . . ,Y4 )](11)

Meisami, Lam, and Van Oyen
where h is a simple sum functionh(y1 , . . . , y4 ) = y1 + y2 + y3 + y4
and Yi ’s are four independent miles-per-gallon (MPG) responses generated from the QRF given the fourfeature values. We use simulation to estimate ψ(F1 , . . . , F4 ) (even though it can be computed exactly).
We should mention that, although we choose to use a very simple h, the same procedure essentiallyapplies when h is the total waiting and overtime cost of a surgery schedule and Yi ’s denote the predicted
surgery durations. Similarly, h can be the overall back-order and holding cost in a multi-product newsvendor
problem and Yi ’s denote the predicted product demands.
We first show two results (Figures 1 and 2) following the analysis in Wager, Hastie, and Efron (2014).
Figure 1 shows the point estimate of (11), given by (7), and the standard deviation estimate, given by (10),on 20 instances of (11) with different feature values. We use one intra-tree simulation replication in our
experiment. The behaviors of these values show similar patterns as Figure 1 in Wager, Hastie, and Efron(2014), except that the variation of the standard deviation estimate appears more uniform regardless of the
distance to the actual values. This may be explained by the summation operation h that gives an averaging
effect.

Figure 1: Point and standard deviation estimates on (11) with the “Auto MPG” data set, following therepresentation of Figure 1 in Wager, Hastie, and Efron 2014. We depict the results on 20 different instancesof (11). The dot and the half-length of each interval represent the point and standard deviation estimates
respectively.
Next, Figure 2 shows a comparison between the IJ variance estimate and the estimate generated bydirect bootstrapping. The latter requires a double-layer bootstrap, first on the generation of the resample
used to build QRF, and second on the generation of the trees within each QRF. We use a resample sizeof 1, 000 in each layer and one intra-tree simulation replication. We use B = 2000 and also one intra-tree
simulation replication in the IJ estimate. Figure 2 shows that the IJ variance estimate is generally close to
the direct bootstrap, suggesting the validity of the former.
We now consider varying the intra-tree simulation replication size, denoted R, as discussed at the endof Section 3. We consider three different scenarios under a total simulation budget of 2000 replications,first with 200 trees and 10 intra-tree replications, second with 400 trees and 5 intra-tree replications, andthird with 2000 trees and only one intra-tree replication. Table 1 shows that the variance estimates areMeisami, Lam, and Van Oyen

Figure 2: Variance estimates using IJ and direct boostrapping, on 20 different instances of (11). Solid line
shows IJ, and dashed line shows direct bootstrapping. The half-height of the gray band is the standard
deviation of the bootstrap distribution in using the direct bootstrapping.
most uniform in the B = 2000 case, which are also significantly higher than the other two cases. Since
Figure 2 demonstrates that the estimates using B = 2000 matches closely with direct bootstrapping, they
can serve as a reasonable benchmark. This implies that increasing R under a fixed budget in this example
may degrade the estimation reliability.
5

CONCLUSION

We have studied the estimation of the variance contributed from input noise, as a way to quantify input
uncertainty and construct a confidence interval, in stochastic simulation driven by input models calibrated
from QRFs. Standard bootstrapping on this problem requires double layers of resampling and simulation
runs. We study the integration of an IJ scheme recently proposed in the context of prediction problems,
which avoids the use of double-bootstrapping, into the simulation-based estimation of the input-induced
variance. We present a strategy and demonstrate it through an elementary numerical example. Future work
includes the mathematical analysis and more extensive testing of our approach.

ACKNOWLEDGMENTS
We gratefully acknowledge support from the National Science Foundation under grants CMMI-1542020,
CMMI-1523453, CMMI-1548201 and CAREER CMMI-1653339.
REFERENCES
Barton, R. R. 2012. “Tutorial: Input Uncertainty in Output Analysis”. In Proceedings of the 2012 Winter
Simulation Conference, edited by C. Laroque, J. Himmelspach, O. R. R. Pasupathy, and A. Uhrmacher,
1–12. Piscataway, New Jersey: IEEE.
Breiman, L. 1996. “Bagging Predictors”. Machine Learning 24 (2): 123–140.
Breiman, L. 2001. “Random Forests”. Machine Learning 45 (1): 5–32.
Büchlmann, P., and B. Yu. 2002. “Analyzing Bagging”. Annals of Statistics 1:927–961.

Meisami, Lam, and Van Oyen
Table 1: Variance estimates on (11) with different combinations of tree and intra-tree replication sizes
under fixed simulation budget.
Variance Estimates
Mean Sum

B = 200 & R = 10

B = 400 & R = 5

B = 2000 & R = 1

106.5696.97
112.54121.85106.73108.49
99.24
97.87
87.07
69.25
76.76
69.55
75.04
80.46
76.76
79.70
86.58
111.49113.14108.55

0.46
1.85
11.99
0.66
0.63
0.50
0.30
0.71
0.44
1.90
3.52
0.94
0.44
0.54
0.27
0.52
0.25
0.12
0.24
0.33

0.95
1.07
0.36
0.48
0.59
0.34
0.50
0.43
0.76
0.48
0.66
0.30
0.64
0.66
1.02
0.43
0.95
0.47
0.37
0.23

2.06
2.30
2.28
2.10
2.08
2.17
2.07
2.09
2.04
2.15
2.16
2.16
2.29
2.12
2.10
2.07
2.01
2.24
2.10
2.15

Cheng, R. C., and W. Holland. 1997. “Sensitivity of Computer Simulation Experiments to Errors in Input
Data”. Journal of Statistical Computation and Simulation 57 (1-4): 219–241.
Cheng, R. C., and W. Holland. 1998. “Two-point Methods for Assessing Variability in Simulation Output”.
Journal of Statistical Computation Simulation 60 (3): 183–205.
Cheng, R. C., and W. Holland. 2004. “Calculation of Confidence Intervals for Simulation Output”. ACM
Transactions on Modeling and Computer Simulation (TOMACS) 14 (4): 344–362.
Efron, B. 2014. “Estimation and Accuracy After Model Selection”. Journal of the American Statistical
Association 109 (507): 991–1007.
Hampel, F. R., E. M. Ronchetti, P. J. Rousseeuw, and W. A. Stahel. 2011. Robust Statistics: The Approach
Based on Influence Functions, Volume 114. New York: John Wiley & Sons.
Huber, P. J. 2011. Robust Statistics. Berlin Heidelberg: Springer.
Jiang, G., L. J. Hong, and B. L. Nelson. 2016. “A Simulation Analytics Approach to Dynamic Risk
Monitoring”. In Proceedings of the 2016 Winter Simulation Conference, edited by T. M. K. Roeder,
P. I. Frazier, R. Szechtman, E. Zhou, T. Huschka, and S. E. Chick, WSC ’16, 437–447. Piscataway,
NJ, USA: Institute of Electrical and Electronics Engineers, Inc.
Lam, H. 2016. “Advanced Tutorial: Input Uncertainty and Robust Analysis in Stochastic Simulation”. In
Winter Simulation Conference (WSC), 2016, edited by T. M. K. Roeder, P. I. Frazier, R. Szechtman,
E. Zhou, T. Huschka, and S. E. Chick, 178–192. Piscataway, NJ, USA: IEEE.
Meinshausen, N. 2006. “QuantileRegression Forests”. Journal of Machine Learning Research 7 (Jun):
983–999.

Meisami, Lam, and Van Oyen
Meisami, A., H. Lam, M. Van Oyen, C. Stromblad, and N. Kastango. 2017. “An Individualized Learning
Methodology for the Surgery Scheduling Problem”. Working paper, University of Michigan Industrial
and Operations Engineering Dept..
Mentch, L., and G. Hooker. 2016. “Quantifying Uncertainty in Random Forests via Confidence Intervalsand Hypothesis Tests”. Journal of Machine Learning Research 17 (26): 1–41.
Song, E., and B. L. Nelson. 2015. “Quickly Assessing Contributions to Input Uncertainty”. IIE Transactions 47 (9): 893–909.
Song, E., B. L. Nelson, and C. D. Pegden. 2014. “Advanced Tutorial: Input Uncertainty Quantification”.
In Proceedings of the 2014 Winter Simulation Conference, edited by A. Tolk, S. Diallo, I. Ryzhov,
L. Yilmaz, S. Buckley, and J. Miller, 162–176. Piscataway, New Jersey: IEEE.
Wager, S., and S. Athey. 2017. “Estimation and Inference of Heterogeneous Treatment Effects Using
Random Forests”. Journal of the American Statistical Association (just-accepted).
Wager, S., T. Hastie, and B. Efron. 2014. “Confidence Intervals for Random Forests: The Jackknife andthe Infinitesimal Jackknife”. Journal of Machine Learning Research 15 (1): 1625–1651.
AUTHOR BIOGRAPHIESAMIRHOSSEIN MEISAMI is a Ph.D. candidate in the Department of Industrial and Operations Engineering at the University of Michigan, Ann Arbor. His research interests include data-driven decision
making and integration of machine learning and optimization. His e-mail address is meisami@umich.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 at Harvard University, and wasa faculty member in the Department of Mathematics and Statistics at Boston University and the Department
of Industrial and Operations Engineering at the University of Michigan before joining Columbia. His
research focuses on stochastic simulation, risk analysis, and simulation optimization. His email address is
khl2114@columbia.edu.
MARK VAN OYEN Mark Van Oyen is a Professor of Industrial and Operations Engr. (IOE) at the University of Michigan. His interests include the analysis, design, and control of stochastic systems (models
and applications). His current research emphasizes healthcare operations and medical decision making.
He co-authored papers that won the 2016 Manufacturing and Service Operations Management (MSOM)Best Paper Award, 2016 MSOM Service Science SIG Best Paper Award, and the 2010 Pierskalla Award.
He has served as Associate Editor for Operations Research, Naval Research Logistics, IIE Transactions,and IIE Transactions on Healthcare Syst. Engr. and Senior Editor for Flexible Services & Manufacturing.
He was a faculty member of the Northwestern Univ. Sch. of Engr. (1993-2005) and Loyola Univ. of
Chicago’s Sch. of Bus. Admin. (1999-2005). His email address is vanoyen@umich.edu.