Proceedings of the 2022 Winter Simulation Conference
B. Feng, G. Pedrielli, Y. Peng, S. Shashaani, E. Song, C.G. Corlu, L.H. Lee, E.P. Chew, T. Roeder, andP. Lendermann, eds.
BATCHING ON BIASED ESTIMATORS
Shengyi He
Henry Lam
Department of Industrial Engineering and Operations Research
Columbia University
500 West 120th Street
New York, NY 10027, USA

ABSTRACT
Existing batching methods are designed to cancel the variability parameter but not the bias of estimators,
and thus are applied typically in the setting of unbiased estimation. We provide a batching scheme thatcancel out the bias and variability parameters of estimators simultaneously, yielding asymptotically exact
confidence intervals for biased estimation problems. We apply our batching method to finite difference
estimators. We extend our method to the multivariate case in constructing confidence regions. We validate
our theory and analyze the effect of the number of batches through numerical examples.
1

INTRODUCTION

Batching methods are widely used in simulation analysis. The basic idea of batching methods is to divide thedata into batches. To construct a confidence interval (CI), these methods utilize pivotal statistics that cancelout the variability parameter by judiciously combining the batch estimates. They thus construct CIs without
knowing the values of these variability parameters. These schemes are especially useful for problems wherethe variance of the output is hard to compute, such as quantile (Nakayama 2014) whose variance estimation
involves density estimation, and in serially dependent problems and steady-state estimation (Asmussen and
Glynn 2007; Nakayama 2007).
The current literature on batching methods focus on the estimation variability, but not the bias. To
construct CIs for biased estimators, however, one could consider a batching strategy for the bias as well.
Biased estimation arise in many places in the simulation literature, a basic example being the class of finitedifference estimators, which is widely used in stochastic gradient estimation under black-box settings where
the underlying dynamics of the simulation model is inaccessible (Lam et al. 2021). For finite-difference
estimators, under the arguably optimal selection of the perturbation parameter, the asymptotic distributionof the (scaled) estimation error is not centered at 0 (Glynn 1989). Moreover, the center of this asymptotic
distribution depends on model characteristics that are unknown a priori, so an explicit correction for thebias may not be straightforward. In this case, classical batching methods would fail because they do notaccount for this bias. A bootstrap correction for the bias will not work either, since the estimator based on
the resampled data has the same amount of bias as the original data. An explicit correction for the bias
using some type of bias estimators is discussed in Section 4.4 of Hall (1992) for kernel density estimators(KDE). However, the bias estimator requires the order and differentiability of the kernels, which make it
hard to be applied to general biased estimators.
In this paper, we propose a batching method that aims to cancel out both the bias and variance parameterssimultaneously, so that it can be applied to biased estimation problems. The way we handle the bias is totake batch estimates with different batch sizes and cancel out the bias by scaling and subtracting these batch
estimates. The use of different batch sizes has some resemblance with subsampling techniques (Politis

He and Lam
et al. (1999); Bickel et al. (1997)) where the resample size is chosen smaller than the original size in abootstrap, but the motivation in our case is different as we aim to cancel out the bias, instead of removing
smoothness conditions in a full-size bootstrap or handling serial dependence. The variability parameter,on the other hand, can be handled in a similar way as the usual batching methods. Combining them,
we construct a pivotal statistic with a limiting t-distribution, which yields a CI with asymptotically exact
coverage. We also extend our method to the multivariate case where we construct ellipsoidal confidence
regions. We apply our method to finite-difference estimators both in the univariate and multivariate cases,
where the latter is based on simultaneous perturbation.
The rest of the paper is organized as follows. Section 2 introduces biased problems. Section 3 discusses
the challenges of usual batching and bootstrap methods. Section 4 proposes our batching method for biased
problems. Section 5 extends our method to the multivariate case. Section 6 gives numerical examples.
2

BIASED ESTIMATION PROBLEMS

Consider an unknown target quantity ψ. Given n i.i.d. data X1 , . . . , Xn , we use an estimator ψ̂n that depends
on the data and a tuning parameter δn to estimate ψ. Here, we consider biased problems in which wecannot obtain unbiased samples easily.
A basic example of biased estimators is finite-difference estimators in stochastic gradient estimation.
Suppose f (·) is a black-box function that can only be noisily evaluated via unbiased samples. To estimate
the derivative of f (·) at a point x0 , we would use the finite-difference estimator given by
Forward Finite Difference (FFD): ψ̂n,FFD =orCentral Finite Difference (CFD): ψ̂n,CFD =(i)(i)
1 n fˆ1 (x0 + δn ) − fˆ2 (x0 )
∑
n i=1
δn

(i)
(i)
1 n fˆ1 (x0 + δn ) − fˆ2 (x0 − δn )
∑
n i=1
2δn

(i)
(i)
where fˆ1 (·) and fˆ2 (·), i = 1, 2, . . . , n are two independent noisy runs of f , and δn is a perturbation
parameter.
For biased problems such as the finite-difference estimators depicted above, we typically have a central
limit theorem (CLT)nβ (ψ̂n − ψ) ⇒ N µ, σ 2(1)

for some β , when δn is appropriately chosen in terms of n (Fox and Glynn 1989). Here, unlike standard
CLT, because of the bias we have a non-zero-mean normal limit and β could be smaller than the canonicalrate of 1/2. For example, for FFD and CFD, under enough smoothness conditions, we have that (Glynn
1989)



δ f ′′ (x0 ) 2Var fˆ2 (x0 )1/4′,, if n1/4 δn → δ(2)nψ̂n,FFD − f (x0 ) ⇒ N −2
δ2
by choosing δn of order n−1/4 andδ 2 f ′′′ (x0 ) 2Var fˆ2 (x0 )n1/3 ψ̂n,CFD − f ′ (x0 ) ⇒ N −,, if n1/6 δn → δ
3
δ2

(3)

by choosing δn of order n−1/6 . The convergence rates in (2) and (3) are optimal for the respective classesof estimators. That is, if we choose δn at different rates than depicted, then the convergence rates would
not improve. On a high level, this is because the chosen orders of δn balance the bias and variance in theestimation, and any distortions would increase either bias or variance that leads to a larger magnitude ofoverall error.

He and Lam
Our goal is to construct a CI for ψ by using (1). In many problems such as finite differences as shown
above, it is difficult to estimate σ 2 and µ and thus a direct use of the CLT in CI construction. Moreover,as we will show, classical batching methods and bootstrap methods also face challenges.
3

CHALLENGES IN STANDARD BATCHING AND BOOTSTRAP METHODS

There are several variants of batching methods depending on the construction of the point estimator andthe scheme in aggregating different batches. The ideas of these methods are similar in the sense that they
all aim to cancel out the variability of the point estimator via dividing by a quantity that represents thevariability among the batch estimates. We discuss one of the variant in detail here and similar discussion
will also hold for other variants of batching methods.
One common batching method constructs a level 1 − α CI as
Sbatch(4)ψ̄ ± tm−1,1−α/2 √
m
where
ψ̄ =
2
Sbatch=1 m (i)∑ ψ̂n ,m i=1

2
1 m  (i)
ψ̂
−
ψ̄
,
n
∑
m − 1 i=1(i)ψ̂n are independent estimators of ψ, and tm−1,1−α/2 is the 1 − α quantile of t distribution with m − 1
degrees of freedom. This CI is based on the observation that when µ = 0 in (1), we have that√m (ψ̄ − ψ)r(5) ⇒ tm−1
m1
m−1 ∑i=1(i)ψ̂n − ψ̄2

as n → ∞. In (5), the asymptotic distribution does not depend on σ since it is canceled out by taking theratio. However, if µ deviates from 0, the asymptotic distribution of the above statistic would change sinceµ is not canceled out. Therefore, the CI given in (4) is not valid when µ ̸= 0.
Bootstrap methods also have issues due to the bias. Suppose we want to use bootstrap to approximate
the distribution of nβ (ψ̂n − ψ) where ψ̂n is calculated based on i.i.d. data X1 , . . . , Xn . The general principleof bootstrap suggests to use the distribution ofnβ (ψn∗ − ψ̂n )conditional on X1 , . . . , Xn , where ψn∗ is the estimate based on resampled data {X1∗ , . . . , Xn∗ }, sampled withreplacement from {X1 , . . . , Xn }. In the finite difference setting (say FFD), we have thatψ̂n,FFD = EP̂ [X]
where
P̂ :=
and Xi is the i-th estimator:
Xi =

1 n
∑ δXi
n i=1

(i)
(i)
fˆ1 (x0 + δn ) − fˆ2 (x0 ).
δnHe and Lam
Here δXi is point mass at Xi . So the estimator based on resampled data is given by∗
ψn,FFD= EP∗ [X]

where
P∗ :=

1 n
∑ δXi∗ .
n i=1∗
as a function of P∗ , it is clear that
From the linearity of ψn,FFDhi
∗E nβ ψn,FFD− ψ̂n,FFD |{X1 , . . . , Xn } = 0.

But on the other hand, from (2) we know that (under some conditions, e.g., uniform integrability)hiδ f ′′ (x0 )E nβ (ψ̂n,FFD − ψ) → −
̸= 0.
2


∗
Comparing the above two displayed relations, we conclude that the distirbution of nβ ψn,FFD− ψ̂n,FFD
conditional on the data cannot approximate the true distribution of nβ (ψ̂n,FFD − ψ).
Essentially, the reason that the usual bootstrap cannot work is that any estimator with the same choiceof δn will share the same amount of bias. Therefore, comparisons between the original estimator and theestimator based on resampled data will not give any information about the bias. In Section 4.4.2 of Hall
(1992), a similar issue is discussed for KDE, which is a related type of biased problems whose amountof bias is determined by the bandwith. Hall (1992) proposes an explicit correction for the bias, but the
proposed formula for the correction is involved and depends on properties of the problem, such as theorder of the kernel and the differentiability of the kernel function, which makes it hard to implement.
The above challenges faced by the standard batching and bootstrap suggests that we look for ways tocancel out the bias and variance simultaneously. An idea is to construct batch estimators with different
amount of bias so that information on the bias can be obtained. This motivates the batching method that
we will describe next.
4

A BATCHING METHOD FOR BIASED ESTIMATION

We consider a modified batching approach that handles the bias and variability of biased estimation
simultaneously. The proposed method is given in Algorithm 1. In Algorithm 1, we first construct a new( j)quantity νn,s = nβ ψ̂n − sβ ψ̂s which is governed by a (scaled) unbiased CLT. To explain, notice that replacing
n with s in (1), we also have thatsβ (ψ̂s − ψ) ⇒ N µ, σ 2(6)
as s → ∞. By subtracting (1) by (6), we have that for independent ψ̂n and ψ̂s ,

 


nβ ψ̂n − sβ ψ̂s − nβ − sβ ψ ⇒ N 0, 2σ 2 .
( j)

(7)

as n → ∞, s → ∞. Therefore, up to a scaling of nβ − sβ , νn,s is governed by a CLT with mean zero in the( j)limit. The rest of the steps in Algorithm 1 handle the variance using the variability among {νn,s } j=1,2,...,min a similar way as usual batching. We have the following theorem regarding the asymptotic exactness ofthe output of Algorithm 1.
Theorem 1 Suppose that (1) holds. Then in Algorithm 1, we have that for any α ∈ (0, 1), and I =Iα , Iα,lower or Iα,upper ,P (ψ ∈ I ) → 1 − αas n, s → ∞.

He and Lam
Algorithm 1 Batching for biased estimation
Require: nominal level 1 − α, batch sizes n ̸= s, number of batches m, the value of β in (1)for j = 1, 2, . . . , m (independently) do( j)ψ̂n ← an estimator based on n i.i.d. samples and tuning parameter δn( j)ψ̂s ← an estimator based on another s i.i.d. samples and tuning parameter δs
( j)
( j)
( j)
νn,s ← nβ ψ̂n − sβ ψ̂s
end for( j)ν̄ ← m1 ∑mj=1 νn,s

2
( j)
1
S2 ← m−1∑mj=1 νn,s − ν̄


h
 
i
Iα ← ν̄ − tm−1,1−α/2 √Sm / nβ − sβ , ν̄ + tm−1,1−α/2 √Sm / nβ − sβiIα,lower ← −∞, ν̄ + tm−1,1−α √Sm / nβ − sβh Iα,upper ← ν̄ − tm−1,1−α √Sm / nβ − sβ , ∞return Two sided CI Iα,lower and one-sided CIs Iα,lower , Iα,upper .
Proof.

( j)

From the construction of νn,s and the observation (7), we have that


 
( j)
νn,s − nβ − sβ ψ⇒ (Z j ) j=1,...,m(8)

j=1,...,m
i.i.dwhere Z j ∼ N 0, 2σ 2 . Therefore, we have thatν̄ − nβ − sβ ψZ̄d√
√ = tm−1⇒
S/ m
SZ / m(9)2

where Z̄ = (1/m) ∑mj=1 Z j and SZ2 = (1/(m − 1)) ∑mj=1 (Z j − Z̄) , and tm−1 is a t-distribution with degree of
freedom m − 1. The convergence in distribution follows from (8) and the continuous mapping theorem,dand the equality in distribution = follows from elementary properties of independent normal variables.
By inverting (9), we derive the desired result.
The use of the batch size s that is smaller than n has some resemblance with subsampling techniques(Politis et al. (1999); Bickel et al. (1997)) where the resample size in a bootstrap is chosen smaller than thefull data size. However, the motivation here is different from that literature. While they use subsampling
to remove smoothness conditions required by full-size bootstraps or to handle serial dependence, we use
different batch sizes in order to combine these batches in a way to eliminate the bias. Suppose further thats = ρn, i.e., ρ is the subsample factor of s with respect to n. In this case, the two-sided CI in Algorithm 1
can be written as1
SYȲ ± tm−1,1−α/2 √1 − ρα
m
( j)

( j)

2

where Y j = ψn − ρ α ψs , and Ȳ = (1/m) ∑mj=1 Y j and SY2 = (1/(m − 1)) ∑mj=1 (Y j − Ȳ ) .
5

GENERALIZATION TO THE MULTIVARIATE CASE

The analysis in Section 2 can be extended to the multivariate case, where we want a confidence region for
multiple biased estimators. Suppose that we are interested in a multidimensional unknown target quantity
Ψ ∈ Rd and have biased estimators given by Ψ̂n which depend on a tuning parameter δn . In this case, wewill start from a multivariate analog of (1):nβ (Ψ̂n − Ψ) ⇒ N(µ, Σ)(10)He and Lam
where µ ∈ Rd , Σ ∈ Rd×d .
For example, for the gradient estimation of a performance function w.r.t. multiple parameters, simultaneous perturbation finite difference (SPFD) can be employed (Spall 1992). Suppose f (·) : Rd → R can
only be noisily evalutaed and we want to estimate the gradient of f at x0 ∈ Rd . Then we would take arandom vector h ∈ Rd and then for each of i = 1, 2, . . . , d, an estimator for ∇ f (x0 )i is given byfˆ1 (x0 + δn h) − fˆ2 (x0 − δn h)
2δn hi(11)

Therefore, with n simulation runs, the estimator is given by
1 n
Ψ̂n,SPFD = ∑
n j=1
( j)

( j)
( j)
fˆ1 (x0 + δn h( j) ) − fˆ2 (x0 − δn h( j) )

!
(12)

( j)

2δn hii=1,2,...,d( j)where fˆ1 (·) and fˆ2 (·), j = 1, 2, . . . , n are two independent noisy runs of f , h( j) , j = 1, 2, . . . , n are
independent samples with the same distribution as h, and δn is a (deterministic) perturbation parameter. As
discussed in Lam et al. (2019), if we assume that h has mean zero, the components of h are independent,
each component of h has finite inverse second moment and smoothness conditions hold for f , then thebias and variance of the estimator given in (11) are of orders δn2 and δn−2 respectively, which is the sameas CFD. Therefore, when δn is of order n−1/6 , following the general result provided in Theorem 2 of Fox
and Glynn (1989) and using Cramér–Wold device, we have that under regularity conditions,n1/3 (Ψ̂n,SPFD − ∇ f (x0 )) ⇒ N(µ, Σ)for some µ ∈ Rd and Σ ∈ Rd×d .
For the multivariate case, an analog of Algorithm 1 is provided in Algorithm 2. In Algorithm 2, the way
we handle the bias and variance is similar to the univariate case, which yields a statistic with asymptotic
Hotelling T 2 distribution.
Algorithm 2 Batching for biased estimation (the multivariate version)
Require: nominal level 1 − α, batch sizes n ̸= s, number of batches m, the value of β in (10), dimension dfor j = 1, 2, . . . , m (independently) do( j)Ψ̂n ← an estimator based on n i.i.d. samples and tuning parameter δn( j)Ψ̂s ← an estimator based on another s i.i.d. samples and tuning parameter δs
( j)
( j)
( j)
Vn,s ← nβ Ψ̂n − sβ Ψ̂s
end for( j)V̄ ← m1 ∑mj=1 Vn,s


⊤
( j)
j
1
S ← m−1∑mj=1 Vn,s − V̄ Vn,s − V̄
o
n
⊤

2
2
(Here, Td,m−1,1−αthe 1 − αIα,d ← y ∈ Rd : m V̄ − (nβ − sβ )y S−1 V̄ − (nβ − sβ )y ≤ Td,m−1,1−αquantile of Hotelling T 2 with dimension d and m − 1 degrees of freedom.)
return Confidence region Iα,dFor Algorithm 2, we have that
Theorem 2 Suppose that (10) holds. Moreover, suppose that m ≥ d + 1. Then in Algorithm 2,P(Ψ ∈ Iα,d ) → 1 − αas n, s → ∞.

He and Lam
Proof.

Similar to the one-dimensional case, from (1) andnβ (Ψ̂n − Ψ) ⇒ N(µ, Σ),

we have that




( j)
Vn,s − nβ − sβ Ψ ⇒ (Z j ) j=1,2,...,mwhere Z j ∼ N(0, 2Σ) are independent for each j = 1, 2, . . . , m. Therefore,⊤m V̄ − (nβ − sβ )y S−1 V̄ − (nβ − sβ )y ⇒ mZ̄⊤ S−1
Z Z̄
where Z̄ = 1n ∑mj=1 Z j and
SZ =


⊤
1 m
Z j − Z̄ Z j − Z̄ .
∑m − 1 j=1d2
Note that mZ̄⊤ S−1
Z Z̄ = Td,m−1 . This gives the desired result.

6

NUMERICAL EXAMPLES

6.1 A Toy Example
To illustrate the effectiveness of the method, we study the problem of constructing a 95% CI for the gradientoff (θ ) := θ + θ 2 + sin θ .
at θ = 0. Suppose we can only generate estimators fˆ(θ ) = f (θ ) + ε where the noise ε is a standard
normal. We consider both FFD and CFD. For FFD, we set δn = n−1/4 . For CFD, we set δn = n−1/6 .
We define the total simulation budget as N = m(s + n). We let ρ = 1/2. We test m = 2, 5, 10, 15, 20 and
N = 900, 9000, 90000, 900000. Note that once N and m are specified, s and n are determined accordingly
by the constraint for the simulation budget and the choice of ρ. We replicate the experiments 1000 times
to estimate the coverage probabilities and the half length of the CIs.
The results for FFD and CFD are shown in Table 1 and Table 2 respectively. We can see that for both
FFD and CFD, the empirical coverages are close to the nominal level. Indeed, all of the empirical coverages
reported in Tables 1 and 2 are within the range 95%±2%. The difference in their performances appears
in the half length. For a fixed simulation budget, when m increases, the average half length decreases.
For example, in Table 1, when the simulation budget is fixed as 900 and m increases from 2 to 30, theaverage half length reduces from 22.6 to 2.2. A similar observation also holds for CFD. This can also betheoretically√justified. Indeed, from the form of the CI given in Algorithm 1, we can see that the half length
is of order ( m(nβ − sβ ))−1 . Since β < 1/2, the effect of increasing m supersedes the effect of decreasingn, s when m(n + s) and s/n are fixed. The half length also decreases as the simulation budget increases,
but at a slow rate. For example, for FFD, when the simulation budget increases from 900 to 900000 (103times larger), the average half length (when m = 10) only reduces by 81% (1-0.6/3.2). The half length for
CFD has a faster decrease rate than FFD (which reduces by 90% when the simulation budget increase from
900 to 900000), but the amount of decrease in half length is still small compared to the increase in thesimulation budget. Essentially, this is due to the slow rate in the CLT result (1). From the above analysis,
and also suggested by the order of the half length, we can see that increasing m or increasing m and thesimulation budget simultaneously is more effective than just increasing the simulation budget, if we want
to shorten the length of the CI.√Indeed, when N and m increases proportionally such that n and s are fixed,
the order of length given by ( m(nβ − sβ ))−1 will decrease at rate N −1/2 . On the other hand, if we fix mand increase N, the rate is only N −β . Lastly, comparing CFD and FD, we find that the half length for CFD
is significantly smaller than the half length of FFD. For example, when the simulation budget is 9000 and
m = 10, the lengths for FFD and CFD are 1.8 and 0.4 respectively. This is due to the fact that CFD has afaster rate of convergence than FFD in the CLTs (2) and (3).

He and Lam
simulation budget

m empirical coverage average half length2
95.9% (0.6%)
22.6 (16.7)
595.3% (0.7%)
4.4 (1.6)
900
1095.9% (0.6%)
3.2 (0.8)
20
95.8% (0.6%)
2.5 (0.4)30
94.6% (0.7%)
2.2 (0.3)
2
95.0% (0.7%)
12.7 (9.1)
595.6% (0.6%)
2.5 (0.9)9000
1094.8% (0.7%)
1.8 (0.4)
20
95.0% (0.7%)
1.4 (0.2)
30
96.4% (0.6%)
1.2 (0.2)2
93.8% (0.8%)
7.0 (5.4)
594.0% (0.8%)
1.4 (0.5)90000
1095.7% (0.6%)
1.0 (0.2)
20
93.6% (0.8%)
0.8 (0.1)
30
94.8% (0.7%)
0.7 (0.1)
2
94.7% (0.7%)
3.9 (3.0)
595.3% (0.7%)
0.8 (0.3)
90000010
95.2% (0.7%)
0.6 (0.1)
20
94.3% (0.7%)
0.4 (0.1)
30
94.8% (0.7%)
0.4 (0.1)
Table 1: Empirical coverages and average half lengths for FFD in the toy example. The standard errors
are shown in parentheses.
6.2 Sensitivity Analysis for an M/M/1 Queue
We consider the example in Lam et al. (2019). For an M/M/1 queue, we define the performance function
as the expected value of the average system time of the first 10 customers. The parameters are the arrival
rate λ and the service rate µ. Our goal is to construct 95% CIs for the gradients of the performance w.r.t.
λ and µ.
First, we construct 95% CIs for the two gradients separately. We focus on CFD with δn = n−1/6 .
As in the toy example, the total simulation budget is defined as N = m(s + n) and ρ = 1/2. We test
m = 2, 5, 10, 20, 30 and N = 900, 1800, 3600. We use the true values for the gradients given in Lam et al.
(2019). The estimated empirical coverage probabilities and half length are reported in Table 3. From Table
3, we have similar observations as in the toy example: all of the empirical coverages are close to thenominal level, the half length of the CI decreases as simulation budget increases or as m increases, but fora fixed m, the rate that the half length decreases as the simulation budget increases is slow.
Next, we study the 95% confidence region for the whole gradient vector w.r.t. λ and µ. We use SPFD
and choose the distribution of h as (h1 , h2 ) where h1 and h2 are independent taking 1 and −1 with equal
probability. We choose δn = n−1/6 as in CFD. We test m = 5, 10, 20, 30 and N = 900, 1800, 3600 (note that
we can not take m = 2 as in the one-dimensional case since we need m ≥ d + 1 in Theorem 2).
The empirical coverages and average areas of the confidence regions are reported in Table 4. Again,
we observe that the empirical coverages are close to the nominal level. The average area decreases fasteras the simulation budget increases or as m increases than the one-dimensional case. For example, whenm = 30 and the simulation budget increases from 900 to 3600, the average area becomes 33% (0.015/0.045)
of the original average area. On the other hand, in Table 3, the half length for the gradient w.r.t. µ becomes
61% (0.14/0.23) of the original for the same change of simulation budget. This is because the order of thearea of the confidence region is the square of the order of the half lengths.

He and Lam
simulation budget

m empirical coverage average half length2
95.2% (0.7%)
5.18 (3.72)
594.1% (0.7%)
1.15 (0.42)
900
1094.7% (0.7%)
0.86 (0.21)
20
95.4% (0.7%)
0.72 (0.12)30
94.4% (0.7%)
0.66 (0.09)
2
95.1% (0.7%)
2.50 (1.86)
594.8% (0.7%)
0.53 (0.20)
9000
1095.3% (0.7%)
0.40 (0.09)
20
93.4% (0.8%)
0.33 (0.05)
30
95.3% (0.7%)
0.31 (0.04)
2
94.8% (0.7%)
1.11 (0.85)
594.3% (0.7%)
0.25 (0.09)
90000
1094.6% (0.7%)
0.19 (0.04)
20
95.6% (0.6%)
0.15 (0.03)
30
94.9% (0.7%)
0.14 (0.02)
2
94.9% (0.7%)
0.52 (0.38)
595.6% (0.6%)
0.12 (0.04)
90000010
94.9% (0.7%)
0.09 (0.02)
20
95.2% (0.7%)
0.07 (0.01)
30
94.7% (0.7%)
0.07 (0.01)
Table 2: Empirical coverages and average half lengths for CFD in the toy example. The standard errors
are shown in parentheses.
simulationbudgetmempirical
average half
empirical
average halfcoverage (λ )length (λ )coverage (µ)
length (µ)
2 94.9% (0.7%) 1.67 (1.28) 94.8% (0.7%) 1.83 (1.33)
5 95.4% (0.7%) 0.38 (0.13) 94.2% (0.7%) 0.39 (0.15)
900
10 95.0% (0.7%) 0.28 (0.07) 93.2% (0.8%) 0.29 (0.07)
20 94.0% (0.8%) 0.23 (0.04) 95.5% (0.7%) 0.25 (0.04)
30 94.6% (0.7%) 0.22 (0.03) 94.8% (0.7%) 0.23 (0.03)
2 94.2% (0.7%) 1.34 (1.05) 94.7% (0.7%) 1.40 (1.06)
5 95.1% (0.7%) 0.30 (0.11) 96.1% (0.6%) 0.31 (0.11)
1800
10 95.0% (0.7%) 0.22 (0.05) 94.8% (0.7%) 0.24 (0.06)
20 94.5% (0.7%) 0.19 (0.03) 93.6% (0.8%) 0.20 (0.03)
30 94.9% (0.7%) 0.17 (0.02) 96.3% (0.6%) 0.18 (0.02)
2 95.1% (0.7%) 1.07 (0.81) 95.3% (0.7%) 1.13 (0.86)
5 94.4% (0.7%) 0.23 (0.08) 95.1% (0.7%) 0.25 (0.09)3600
10 94.6% (0.7%) 0.18 (0.04) 96.2% (0.6%) 0.18 (0.04)
20 95.1% (0.7%) 0.15 (0.02) 96.2% (0.6%) 0.15 (0.03)
30 93.9% (0.8%) 0.14 (0.02) 95.1% (0.7%) 0.14 (0.02)
Table 3: Empirical coverages and average half lengths for CFD in the M/M/1 queueing example. Columns
with (λ ), (µ) correspond to gradients w.r.t. λ and µ, respectively. The standard errors are shown inparentheses.

He and Lam
simulation budget

m empirical coverage average area
5
94.7% (0.7%)
0.191 (0.252)
900
1095.4% (0.7%)
0.076 (0.059)
20
95.1% (0.7%)
0.050 (0.024)
30
94.7% (0.7%)
0.045 (0.019)
595.0% (0.7%)
0.115 (0.136)
1800
1095.1% (0.7%)
0.042 (0.029)
20
94.6% (0.7%)
0.030 (0.015)
30
94.2% (0.7%)
0.025 (0.010)
594.7% (0.7%)
0.063 (0.070)
3600
1096.2% (0.6%)
0.026 (0.018)
20
94.4% (0.7%)
0.017 (0.008)
30
93.9% (0.8%)
0.015 (0.006)
Table 4: Empirical coverages and average areas for SPFD in the M/M/1 queueing example. The standard
errors are shown in parentheses.
ACKNOWLEDGMENTS
We gratefully acknowledge support from the National Science Foundation under grants CAREER CMMI1834710 and IIS-1849280.
REFERENCES
Asmussen, S., and P. W. Glynn. 2007. Steady-State Simulation, 96–125. New York, NY: Springer New York.
Bickel, P. J., F. Götze, and W. R. van Zwet. 1997. “Resampling Fewer Than n Observations: Gains, Losses, and Remedies
for Losses”. Statistica Sinica 7(1):1–31.
Fox, B. L., and P. W. Glynn. 1989. “Replication Schemes For Limiting Expectations”. Probability in the Engineering andInformational Sciences 3(3):299–318.
Glynn, P. 1989. “Optimization Of Stochastic Systems Via Simulation”. In 1989 Winter Simulation Conference Proceedings,edited by E. A. MacNair, K. J. Musselman, and P. Heidelberger, 90–105. Piscataway, New Jersey: Institute of Electricaland Electronics Engineers, Inc.
Hall, P. 1992. The Bootstrap and Edgeworth Expansion. Springer.
Lam, H., H. Li, and X. Zhang. 2021. “Minimax Efficient Finite-Difference Stochastic Gradient Estimators Using Black-Box
Function Evaluations”. Operations Research Letters 49(1):40–47.
Lam, H., X. Zhang, and X. Zhang. 2019. “Enhanced Balancing of Bias-Variance Tradeoff in Stochastic Estimation: A Minimax
Perspective”. arXiv:1902.04673. To appear in Operations Research.
Nakayama, M. K. 2007. “Fixed-width Multiple-comparison Procedures Using Common Random Numbers for Steady-state
Simulations”. European Journal of Operational Research 182(3):1330–1349.
Nakayama, M. K. 2014, November. “Confidence Intervals for Quantiles Using Sectioning When Applying Variance-Reduction
Techniques”. ACM Trans. Model. Comput. Simul. 24(4):1–21.
Politis, D., D. Wolf, J. Romano, M. Wolf, P. Bickel, P. Diggle, and S. Fienberg. 1999. Subsampling. Springer Series in Statistics.
Springer New York.
Spall, J. 1992. “Multivariate Stochastic Approximation Using a Simultaneous Perturbation Gradient Approximation”. IEEE
Transactions on Automatic Control 37(3):332–341.

AUTHOR BIOGRAPHIESSHENGYI HE is a PhD student in the Department of Industrial Engineering and Operations Research at Columbia University.
He received his B.S. degree in statistics from Peking University in 2019. His research interests include variance reduction anduncertainty quantification via stochastic and robust optimization. His email address is sh3972@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. His research interests include efficientHe and Lam
methodologies and statistical uncertainty quantification for Monte Carlo computation, predictive modeling and data-driven
optimization. His email address is khl2114@columbia.edu. His website is http://www.columbia.edu/ khl2114/.