偏誤修正模擬法在均態模擬上的效能比較

偏誤修正模擬法在均態模擬上的效能比較

計畫類別： 個別型計畫

計畫編號： NSC93-2213-E-004-011-

執行期間： 93 年 08 月 01 日至 94 年 07 月 31 日

執行單位： 國立政治大學資訊管理學系

計畫主持人： 謝明華

計畫參與人員： 黃瓊玉, 林立人

報告類型： 精簡報告

處理方式： 本計畫可公開查詢

中 華 民 國 94 年 10 月 31 日

行政院國家科學委員會補助專題研究計畫

成果

報告

偏誤修正模擬法在均態模擬上的效能比較

計畫類別：□ 個別型計畫 □ 整合型計畫

計畫編號：NSC 93－2213－E－004－011－

執行期間： 93 年 8 月 1 日至 94 年 7 月 31 日

計畫主持人：謝明華

計畫參與人員：黃瓊玉, 林立人

成果報告類型(依經費核定清單規定繳交)：精簡報告

執行單位： 國立政治大學資訊管理學系

中 華 民 國 94 年 10 月 31 日

中文摘要

當模擬技師模擬一個隨機系統時，他通常是對均衡狀態的效能指標感興趣。對這樣的

指標，古典的點估計方法只簡單地將從所模擬過程的適當的函數做時間平均。因為模擬無

法起始於（未知）的均態分配，所以古典的點估計通常是有偏誤的。

在具回復性 (regenerative) 的均態模擬中。有數種偏誤修正的模擬方法被提出。雖

然這些方法的漸進特性 (asymptotic property) 相當類似。但它們實證效能的比較卻十分

缺乏。在本計畫中，我們將應用這些方法在各種 queueing 和庫存管理的模型上。我們將

研究這些方法實際應用在這些模型上的實證結果。我們研究的問題將包括：

一、點估計偏誤特性。

二、信賴區間的覆蓋準確度。

本計畫的研究成果可成其他模擬技師的有用指南。

關鍵字：偏誤修正估計，回復，模擬，均態估計

I

英文摘要

When simulating a stochastic system, simulationists often are interested in estimating

various steady-state performance measures. The classical point estimator for such a measure

involves simply taking the time average of an appropriate function of the process being simulated.

Since the simulation can not be initiated with the (unknown) steady-state distribution, the

classical point estimate is generally biased. In the context of regenerative steady-state simulation,

a variety of other point estimates have been developed in an attempt to minimize the bias. In this

project we will investigate the results of applying these point estimates to various queueing and

inventory management models. The bias of these point estimates of the performance measure and

the coverage probability of associated confidence intervals will be investigated for these models.

Conclusions will be drawn from this experimental work as to which methods are most effective

in reducing bias.

報告內容

本 計 劃 已 有 期 刊 論 文 發 表 (ACM Transactions on Modeling and Computer Simulation

(TOMACS). 2004. 14: 325-343. ) 故以論文內容作為成果報告內容

。

Empirical Performance of Bias-Reducing

Estimators for Regenerative Steady-State

Simulations

MING-HUA HSIEH

National Chengchi University, Taiwan and

DONALD L. IGLEHART and PETER W. GLYNN Stanford University

When simulating a stochastic system, simulationists often are interested in estimating various steady-state performance measures. The classical point estimator for such a measure involves simply taking the time average of an appropriate function of the process being simulated. Since the simulation can not be initiated with the (unknown) steady-state distribution, the classical point estimator is generally biased. In the context of regenerative steady-state simulation, a variety of other point estimators have been developed in an attempt to minimize the bias. In this paper, we provide an empirical comparison of these estimators in the context of four different continuous-time Markov chain models. The bias of the point estimators and the coverage probabilities of the associated confidence intervals are reported for the four models. Conclusions are drawn from this experimental work as to which methods are most effective in reducing bias.

Categories and Subject Descriptors: G.3 [Probability and Statistics]: probabilistic algorithms; I.6.1 [Simulation and Modeling]: Simulation Theory

General Terms: Algorithms, Performance, Theory

Additional Key Words and Phrases: Bias-reducing estimators, regeneration, simulation, steady-state estimation

1. INTRODUCTION

Let Y = (Y (t) : t ≥ 0) be a real-valued stochastic process representing the output of a simulation. Suppose that Y satisfies a law of large numbers (LLN)

This research was supported by the Army Research Office under Contract no. DAAG55-97-1-0377 and the National Science Council, Taiwan, under Grant no. NSC 93-2213-E-004-011.

Authors’ addresses: M.-H. Hsieh, Management Information Systems Dept., National Chengchi University, No. 64, Zhi-nan Road, Section 2, Wen-shan District, Taipei 116, Taiwan ROC; email: mhsieh@mis.nccu.edu.tw; D. L. Iglehart and P. W. Glynn, Department of Management Science and Engineering, Stanford University, Stanford, CA.

Permission to make digital or hard copies of part or all of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or direct commercial advantage and that copies show this notice on the first page or initial screen of a display along with the full citation. Copyrights for components of this work owned by others than ACM must be honored. Abstracting with credit is permitted. To copy otherwise, to republish, to post on servers, to redistribute to lists, or to use any component of this work in other works requires prior specific permission and/or a fee. Permissions may be requested from Publications Dept., ACM, Inc., 1515 Broadway, New York, NY 10036 USA, fax:+1 (212) 869-0481, or permissions@acm.org.

C

2004 ACM 1049-3301/04/1000-0325 $5.00

of the form ¯ Y (t)= 1 t
t 0 Y (s)ds⇒ α (1)

as t → ∞, for some (deterministic) constant α, where ⇒ denotes convergence

in distribution. The steady-state simulation problem is concerned with the

effi-cient estimation of the steady-state meanα, and the construction of associated

confidence intervals. Because steady-state performance measures are funda-mental to the analysis of an enormous variety of stochastic systems arising in production, inventory, telecommunications, computing, and service industries, the steady-state simulation problem is of great practical importance.

One of the major difficulties associated with steady-state simulation is that the system is typically initialized in a state that is atypical of steady-state behavior. For example, a factory flow simulation is often initialized with no work-in-process on the factory floor. Such an initial condition gives rise to an “initial transient” during which the values of Y observed are not characteristic

of steady-state. As a consequence, the obvious estimator ofα suggested by (1),

namely the time-average ¯Y (t), is biased as an estimator ofα. In other words,

E ¯Y (t)= α.

As might be expected, this bias problem is greatly magnified when paral-lel processors are used to accelerate the steady-state calculation. In particular, if independent replications of Y are simulated on each of the processors, the initial transient will be replicated on each processor, so that the presence of potentially massive parallelism then has no mitigating impact on the initial transient problem; see Glynn and Heidelberger [1991a, 1991b, 1992a] for de-tails. In this article, we focus on the traditional computing environment, in which only a single processor is available.

In such a context, the bias in ¯Y (t) as an estimator ofα (also known as the

“initial bias”) can be mitigated in two different ways. One approach is to delete that initial segment of the simulation that is “contaminated” by initial bias. Such an initial bias deletion approach has been studied by many authors; see, for example, Cash et al. [1992], Glynn [1995], Goldsman et al. [1994], Nelson [1992], Schruben [1982], Schruben et al. [1983], White [1997], and White et al. [2000]. An alternative is to consider an estimator, based on simulating Y over

[0, t], that attempts to compensate for the bias present in ¯Y (t). We refer to such

estimators as “bias reducing” estimators.

In this article, we provide an empirical investigation of six bias-reducing estimators that have been proposed in the context of processes Y that are re-generative. A great variety of stochastic systems exhibit regenerative structure, including discrete and continuous time Markov chains, as well as a broad class of discrete-event simulations; see Glynn [1989] for a discussion of the regen-erative structure of discrete event simulations. Theoretically all “well-posed” steady-state simulation problem are regenerative, see Glynn [1994]. However, identifying the regeneration times can be difficult for a general discrete-event simulation, see Henderson and Glynn [2001]. Therefore, one should notice the limitation that the bias-reducing estimators we studied can only apply to sim-ulations with identifiable regeneration times. This article may be viewed as

Empirical Performance of Bias-Reducing Estimators • 327

an update of an empirical study of Iglehart [1975] for regenerative processes, in which three of the six estimators considered here were investigated on a different time scale.

2. THEORETICAL BACKGROUND

Suppose that Y is the simulation output that is derived from the simulation of a stochastic system that can be modeled in terms of a Markov process X . In particular, suppose that Y (t)= f (X (t)), where X = (X (t) : t ≥ 0) is a Markov process living on a state space S, and f : S→  is a real-valued performance measure. Assuming that X exhibits positive recurrent behavior, it can be shown in substantial generality that

E ¯Y (t)= α +b

t + o(exp(−βt)) (2)

as t → ∞, for some constants b and β (where β > 0), where o(a(t)) denotes a functionψ(t) such that ψ(t)/a(t) → 0 as t → ∞. See, for example, Glynn [1984], for such a result in the setting of finite-state continuous-time Markov chains (CTMC’s).

Here, the time parameter t corresponds to simulated time. In other words, ¯

Y (t) is obtained by averaging the values of Y (·) associated with simulating X to

the deterministic time horizon t. It is important to recognize that different bias expansions hold when the time horizon is specified differently. For example, when X is simulated to the completion of nth regenerative cycle at (simulated) time T (n), one obtains a bias expansion of the form

E ¯Y (T (n))= α + ˜b1 n + ˜b2 n2 + o  1 n2  (3) as n→ ∞, where ˜b2is generally nonzero. See Iglehart [1975] for a discussion

of (3). Note that because ˜b2 = 0, the bias expansion (3) is qualitatively (and

quantitatively) different from the expansion (2). The key point here is that bias expansions are heavily impacted by the time scale which one chooses to use for specifying the termination time for the simulation. In this article, we will focus exclusively on estimators obtained from a given simulated time horizon [0, t], so that the relevant bias expansion is (2).

A corresponding expansion to (2), focused on the variance of ¯Y (t), also exists

in substantial generality: Var ¯Y (t)= σ 2 t + c t2 + o  1 t2  (4) as t→ ∞; see Glynn [1984] for details in the context of continuous-time Markov chains. It turns out that the time-average variance constantσ2is unaffected by the choice of the initial condition for X (although c is impacted by the ini-tialization). Thus, the mean square error (MSE) of ¯Y (t) can be easily computed

via the formula

MSE( ¯Y (t))= Var ¯Y (t) + (E ¯Y (t) − α)2=σ

2 t + b2_{+ c} t2 + o  1 t2  (5)

as t → ∞. Note that the first-order term in the MSE expression, namely σ2_{/t, is}

unaffected by the choice of the initial condition. Thus, the initial transient has only a second-order impact on the quality of the sample mean estimator ¯Y (t).

This suggests that use of a bias-reducing estimator, in preference to ¯Y (t), will

offer only second-order improvements to the quality of the simulation output procedure.

This brings us to the question of how to measure the quality of a particular steady-state estimator. Given an estimator ˆα(t) for the steady-state mean α, this paper will consider the two following measures:

(1) mean absolute percentage error,

100· |E ˆα(t) − α|/|α|; (2) median absolute percentage error,

100· |(median of the distribution of ˆα(t)) − α|/|α|;

We will evaluate both of these quality measures via simulation. In particular, given a simulated time horizon t and an estimator ˆα(t), we perform m indepen-dent iindepen-dentically distributed (i.i.d.) replications of the random variable (r.v.) ˆα(t), thereby producing ˆα1(t),. . . , ˆαm(t), and estimate (1) via

100· | 1 m m i=1αˆi(t)− α| |α| and estimate (2) via

100· |(sample median of { ˆα1(t),. . . , ˆαm(t)}) − α|/|α|

In most carefully designed simulations, the simulator will wish to produce a con-fidence interval forα, in addition to an estimator for α. Approximately valid con-fidence intervals can easily be derived from the central limit theorem (CLT) that holds for all the bias reducing estimators ˆα(t) considered in this paper, namely

√

t( ˆα(t) − α) ⇒ σ N(0, 1) (6)

as t → ∞. The constant σ appearing here is precisely the time-average variance constant of expressions (4) and (5). As noted earlier, σ2 _{does not}

depend on the choice of the initial condition. Consequently, the quality of the confidence interval is affected by the initialization only in a second-order sense. To gain an appreciation for the magnitude of this second-order effect, it is appropriate to consider the coverage error of the confidence interval. Specifi-cally, an approximate 100(1− δ)% confidence interval for α is characterized via a choice of a point estimator ˆα(t) for α and a time-average variance constant estimator ˆv(t) forσ2_{. Given ˆα(t) satisfying (6) and ˆv(t) satisfying ˆv(t) ⇒ σ}2_as

t → ∞, a corresponding 100(1 − δ)% confidence interval for α is given by

 ˆ α(t) − z ·  ˆv(t) t , ˆα(t) + z ·  ˆv(t) t  ,

where z is chosen so that P (−z ≤ N(0, 1) ≤ z) = 1−δ. Denote the above random interval by [L(t), R(t)], so that L(t) and R(t) are, respectively, the left and right end points of the confidence interval. The coverage error associated with the

Empirical Performance of Bias-Reducing Estimators • 329

interval for simulated time horizon t is then given by |P(α ∈ [L(t), R(t)]) − (1 − δ)|.

The one-side coverage error (to the left) is|P(α < L(t)) − δ/2|, whereas the one-sided error (to the right) is|P(α > R(t)) − δ/2|.

In this article, we will numerically investigate the coverage associated with a particular choice of ˆα(t) and ˆv(t), in addition to studying the quality of the point estimator ˆα(t) itself.

3. DESCRIPTION OF THE NUMERICAL EXPERIMENT

In this section, we will provide the details of our numerical experiment. We start by describing three continuous-time Markov chain (CTMC) queueing models and one CTMC inventory model that were used as test vehicles for our simula-tions. We choose these four models both because they are somewhat represen-tative of real applications and because they are simple enough to permit us to compute exactly the theoretical values of certain parameters associated with the simulation.

Model 1. The first model considered is a single-server queueing model with a Poisson arrival process having an arrival rate of 0.95 customers per unit time, and exponential processing times with a mean of 1. Any customer that arrives to find at least 199 customers in the system balks and is permanently lost to the system. Thus, the number-in-system process X = (X (t) : t ≥ 0) is a birth–death process on S = {0, 1, . . . , 199} (often denoted as M/M/1/199). In this model (as in all our others), X (t)⇒ X (∞) as t → ∞. Our goal here is to computeα = EX(∞) = 18.993.

Model 2. The second model is a non-Markovian server analog of our first model. Here the arrival process is Poisson with arrival rate 0.95, and processing times are gamma distributed with parameters (2, 2). Any customer that arrives to find at least 199 customers in the system balks and is permanently lost to the system. The number-in-system process X = (X (t) : t ≥ 0) is certainly not a CTMC. However, a gamma(2,2) random variable is equivalent to the sum of two independent exponential random variables with mean 0.5. If we use

Z = (Z (t) : t ≥ 0) to keep track of the current phase of the service time in

process, i.e. Z (t) = 1 if the server is consuming the first exponential random variable at time t and Z (t) = 2 if the server is consuming the second one at time t. Let X (t)= (X (t), Z (t)). Then, X = ( X (t) : t ≥ 0) is a CTMC. Our goal

here is to computeα = Ef( X (∞)) = 14.487, where f ( X (t)) = X (t).

Model 3. This model is the queue-length process for two M/M/1/19 queues in tandem. Here the arrival process is Poisson with arrival rate 1. The arrived job is served at server 1 first. The service discipline is FCFS at all servers and all servers generate independent identically distributed (i.i.d.) exponentially distributed service times with mean 1/µi. When a job completes service at

server 1, it goes to server 2. If it finds at least 19 customers in the second queue, it is permanently lost to the system. When a job completes service at server 2, it goes to server 1 with probability p and leaves the system with probability

1− p. The parameter values we have selected are µ1 = 2.3, µ2 = 2.15, and

p= 0.5. Let Xi(t), i= 1, 2, denote the number of jobs in server i at time t, and

X (t) = (X1(t), X2(t)). Our performance measure isα = Ef(X (∞)) = 11.523,

where f (X (t))= X1(t)+ X2(t).

Model 4. This model is the number-in-stock process of a single item inven-tory management system. The system is operated under a stationary (s, S) pol-icy. When the number-in-stock decreases to a fixed number s, an order amount of S− s is placed. The delivery times are i.i.d. exponentially distributed with a mean of 4/3. The demands for the item occur according to a Poisson process with arrival rate 1. Demands that cannot be satisfied immediately are lost. We select S= 15 and s = 2. Let X (t) denote the number-in-stock at time t. Our per-formance measure is the steady-state operating costα = Ef(X (∞)) = 401.696, where

f (X (t))=



50X (t), if X (t)> 0; 300, if X (t)= 0.

We turn next to a discussion of the specific estimators forα that we shall con-sider in this article. The generic steady-state simulation problem that we are considering in this article involves the computation of the steady-state meanα of a real-valued process Y = (Y (t) : t ≥ 0). Given Y and a simulated time hori-zon t, our first steady-state estimator is the “classical” (time-average) sample mean, namely α1(t)= 1 t
t 0 Y (s) ds.

As mentioned in the Introduction, all our bias reducing estimators take advan-tage of the presence of regenerative structure. In particular, the times at which the underlying continuous-time Markov chain returns to any fixed state z∈ S is a sequence of regeneration times for Y . Let (T (n) : n≥ 0) be the sequence of times at which Y regenerates, and set

i =

T (i)

T (i−1)Y (s) ds,

τi = T(i) − T(i − 1),

N (t) = max{n ≥ 0 : T(n) ≤ t},

for i≥ 1 and t ≥ 0. A fundamental theoretical assumption for the validity of the bias reducing estimators that we shall study is that the regenerative process

Y be non-delayed, so that a regeneration occurs at t = 0 (and so that T(0) =

0). Thus, once the initial state for the simulation is chosen, this necessarily determines the regenerative return state to be used.

Our first bias-reducing estimator is due to Meketon and Heidelberger [1982]. It is defined as α2(t)= N (t)+1 i=1 i N (t)+1 i=1 τi .

Empirical Performance of Bias-Reducing Estimators • 331

They established theoretically thatα2(t) does indeed reduce bias, as evidenced

by the bias expansion

Eα2(t)= α + o(1/t) (7)

as t → ∞. (Compare with (2).) This bias expansion can be derived by using Taylor series argument, see Glynn and Heidelberger [1990] for such an argu-ment. All the bias-reducing estimators that we shall discuss here have the same theoretical characteristic, namely, that their bias is of order o(1/t) as t → ∞.

The second bias-reducing estimator is obtained by using Taylor expansions for the bias of a general function of sample means based on N (t) observations. It was introduced by Glynn and Heidelberger [1990] and is given by

α3(t)=  α1(t), N (t)= 0; α1(t)+_t12 N (t) i=1 iτi− α1(t) N (t) i=1 τi2 , N (t)≥ 1.

See Glynn and Heidelberger [1990] for a discussion of its bias-reducing properties.

Our third bias-reducing estimator is based on jack-knifing and is given as

α4(t)=    α1(t), N (t)≤ 1; N (t) _{N (t)} i=1i N (t) i=1τi  − N (t)−1 N (t) N (t) i=1  j=ij  j=iτj  , N (t)≥ 2.

The bias-reducing theoretical behavior ofα4(t) was established by Glynn and

Heidelberger [1992b].

Our final bias-reducing estimator is based on a renewal-theoretical asymp-totic expansion for the bias of the sample mean up to time t of a regenerative process. It was proposed in Glynn [1994] and is given by

α5(t)=  α1(t), N (t)= 0; α1(t)+_t12 N (t) i=1  Ai− α1(t)τi2/2  , N (t)≥ 1, where Ai =
_τi 0 sY(T (i− 1) + s) ds.

We also have an interest in studying the degree to which initial bias affects the quality of an estimator. To this end, consider the estimator

α6(t)= α1(t)+ (α − Eα1(t)).

Note thatα6(t) is unbiased as an estimator for α. However, α6(t) can not be

employed practically, because it requires the ability to pre-compute both the steady-state meanα and the “transient” expectation Eα1(t). Of course, for the

class of continuous-time Markov chain models we have selected, both α and

Eα1(t) can be computed numerically. (The quantityα can be found from the

lin-ear equations that characterize the steady-state of a continuous-time Markov chain, whereas Eα1(t) can be computed, as we shall see next, by solving

numer-ically a system of linear differential equations.) By empirnumer-ically studyingα6(t),

we can see the degree of improvement that is theoretically possible when the initialization bias is entirely eliminated.

Suppose that Y takes the form Y (t)= f (X (t)), where X = (X (t), t ≥ 0) is a

S-valued continuous-time Markov chain with finite state space and f : S→ 

is a “performance measure” that we choose to encode as a column vector (i.e.,

f = ( f (x) : x ∈ S) is a column vector.) Let w(t, x) = t Exα1(t)= Ex

t

0 f (X (s)) ds,

where Ex(·) denotes the expectation operator under which X (0) = x. For h > 0,

the Markov property implies that

w(t+ h, x) = Ex

h

0

f (X (s))ds+ Exw(t, X (h))

= f (x)h + o(h) + Exw(t, X (h)).

But Exw(t, X (h))= w(t, x) + (Aw(t))(x)h + o(h), where A = (A(x, y) : x, y ∈ S)

is the generator (or “rate matrix”) of X and w(t)= (w(t, x) : x ∈ S) is encoded as a column vector. Hence,

h−1(w(t+ h, x) − w(t, x)) = f (x) + (Aw(t))(x) + o(1) (8) as h ↓ 0. So, letting h ↓ 0 in (8), we conclude that w(t) exists and equals

f + Aw(t). So, Exα1(t) can be computed by solving the initial value problem

w(t)= f + Aw(t) subject to w(0)= 0.

This can be solved numerically by a finite difference method. An alternative is available when X is irreducible. Then, X possesses a unique stationary distri-butionπ = (π(x) : x ∈ S), which we encode as a row vector. If is the matrix in which each row equalsπ, it is known that (A − )−1exists and

w(t)= (A − )−1[(exp( At)− I)( f − f )] + f t; See Glynn [1994].

As indicated in Section 2, the quality of the above estimators needs to be as-sessed. Based on Section 2’s discussion, we will compute, for each combination of estimator and model, the mean absolute percentage error and the median absolute percentage error. This calculation will be based on 5000 independent replications of the model. In addition, our experiments are concerned with cov-erage error, and the degree to which covcov-erage error is affected by initialization issues. Thus, we will also numerically investigate the coverage errors described in Section 2.

4. NUMERICAL RESULTS

We have run simulations for a variety of cases by varying the time interval. Since we are interested in evaluating the performance of bias-reducing steady-state estimators, we selected, for each model, three combinations of initial state and simulated time that have non-trivial initial bias (1% to 25%), see Figures 1–4 for details. For each combination of simulated time and model, we run 5000 simulations. For each simulation run, we compute the estimatorsαi,

i = 1, . . . , 6, described in Section 3. In addition, we also run another set of

Empirical Performance of Bias-Reducing Estimators • 333

Fig. 1. Mean APE and median APE for various estimators in Model 1 withλ = .95 and µ = 1 and initialized at X (0)= 0. Here, t is the simulated time horizon (chosen to be 64Eτ1, 256Eτ1, and

1024Eτ1) for estimatorsα1(t) and estimatorsα3(t) toα7(t). Forα2(t), the simulated time t2is the

average (over 5000 replications) of the time required to complete the cycle in progress at t.

Recall thatα1(t) is the classical point estimator which having “high-bias” (bias

of order 1/t, see Eq. (2)); α2(t) to α5(t) are “low-bias” estimators (i.e., bias of

order o(1/t), see Eq. (7)); and α6(t) andα7(t) are unbiased estimators. The

esti-matorsα2(t) toα6(t) can be regarded as “bias-reducing” estimators andα6(t) is

the “best” one. Initial transient deletion methods aim to collect sample paths that are approximately in steady-state. Viewed in this light,α7(t) and other

initial transient deletion methods can be regarded as “approximated steady-state” estimators andα7(t) is the “best” one. Therefore, we may consider that

we have done comparisons between bias-reducing estimators and a “perfect” initial transient deletion estimator (i.e.,α7(t)).

As mentioned in Section 2, we are interested in mean absolute percentage error (APE) for these estimators. In Figures 1–4 we displayed the simulation results for mean APE. The results are as expected: compared to the classical estimatorα1(t), the low-bias estimators generally reduced the bias. The only

exception isα4(t) which has more bias thanα1(t) for some cases in Model 1 to

3. This exception theoretically can relate to the fact that Eq. (7) (theoretical property of the low-bias estimators we considered) implies that the low-bias estimators must reduce bias when the chosen t is sufficiently large. Empirically, this exception can be explained as follows:α4(t) only uses simulated information ACM Transactions on Modeling and Computer Simulation, Vol. 14, No. 4, October 2004.

Fig. 2. Mean APE and median APE for various estimators in Model 2 withλ = .95 and gamma(2,2) distributed processing times and initialized at X (0)= 0. Here, t is the simulated time horizon (chosen to be 64Eτ1, 256Eτ1, and 1024Eτ1) for estimatorsα1(t) and estimatorsα3(t) toα7(t). For α2(t), the simulated time t2is the average (over 5000 replications) of the time required to complete

the cycle in progress at t.

up to T (N (t)) which is significant shorter than t for the exceptional cases, for example, (ET(N (t)), t) ≈ (990, 1347) for the exceptional case in Model 1 (see Figure 1). In short, the empirical results indicate that the low-bias estimators are generally effective in reducing mean APE.

We are also interested in the bias-reducing effect on median APE of these estimators. The experimental results are reported in Figures 1–4. From our empirical study, we found that the bias-reducing effect on median APE of these estimators is model dependent. For Models 1 and 2, the low-bias estimators have limited use in reducing median APE. This is because the distribution of

α1(t) of these two models is highly skewed, and thus the magnitude of median

APE of α1(t) is much larger than that of mean APE of α1(t). Note that even

the unbiased estimators have the same limitation for these two models. For Models 3 and 4, the low-bias estimators are effective in reducing median APE, since the distributions of α1(t) of these two models are more symmetric, and

the magnitude of median APE of α1(t) and mean APE ofα1(t) are about the

same. In brief, the empirical results indicate that the low-bias estimators can effectively reduce median APE if the distribution ofα1(t) is near symmetric and

Empirical Performance of Bias-Reducing Estimators • 335

Fig. 3. Mean APE and Median APE for Various Estimators in Model 3 withλ = 1, µ1 = 2.3, µ2 = 2.15, and p = 0.5 and initialized at X (0) = (0, 0). Here, t is the simulated time horizon

(chosen to be 8Eτ1, 16Eτ1, and 32Eτ1) for estimatorsα1(t) and estimatorsα3(t) toα7(t). Forα2(t),

the simulated time t2is the average (over 5000 replications) of the time required to complete the

cycle in progress at t.

We turn next to the coverage issues of the confidence interval forα. Before proceeding to the experimental results, let us introduce the time-average vari-ance constant (TAVC) estimators we used.

The first TAVC estimatorσ2 _{is its theoretical value. We use it as a}

perfor-mance benchmark.

The second TAVC estimator t· Var(αi(t)) is based on Eq. (4). Here, Var(αi(t))

is the sample variance ofαi(t). That is,

 Var(αi(t))= 1 4999 5000 k=1 α(k) i (t), for i= 1, . . . , 7,

whereα(k)_i (t) is the kth replication ofαi(t).

The third TAVC estimator vi(t) is defined as

vi(t)=    0, N (t)= 0, 1; N (t) j=1(j−αi(t)τj)2 N (t) j=1τj , N (t)≥ 2.

for i = 1, . . . , 6. It is the classical TAVC estimator for the regenerative method.

Fig. 4. Mean APE and median APE for various estimators in Model 4 with s = 2 and S = 15 and initialized in State S= 15. Here, t is the simulated time horizon (chosen to be 2Eτ1, 4Eτ1,

and 8Eτ1) for estimatorsα1(t) and estimatorsα3(t) toα7(t). Forα2(t), the simulated time t2is the

average (over 5000 replications) of the time required to complete the cycle in progress at t.

The last TAVC estimator t· vJ(t)/N(t) is based on jackknife theory and vJ(t)

is defined as vJ(t)=  0, N (t)= 0, 1; 1 N (t)−1 N (t) j=1( ˜αj(t)− α4(t))2, N (t)≥ 2, where ˜ αj(t)= N(t) N (t) i=1 i N (t) i=1 τi  − (N(t) − 1)  i= ji  i= jτi  .

It is noteworthy that only the third and fourth TAVC estimators can be obtained from a single replication. Using these four TAVC estimators, we con-structed four different 100(1− δ)% confidence intervals of α. Namely,

 αi(t)− z σ √ t,αi(t)+ z σ √ t  , (9)  αi(t)− z   Var(αi(t)),αi(t)+ z   Var(αi(t))  , (10)

Empirical Performance of Bias-Reducing Estimators • 337  αi(t)− z  vi(t) t ,αi(t)+ z  vi(t) t  , (11) and  αi(t)− z  vJ(t) N (t),αi(t)+ z  vJ(t) N (t)   , (12)

where z is chosen so that P (−z ≤ N(0, 1) ≤ z) = 1−δ. Let [L(t), R(t)] denote any above random interval. To access the coverage quality of confidence intervals, we computed the following measures:

m= ˆP(α ∈ [L(t), R(t)]) = 1 5000 5000 k=1 Iα ∈ L(k)(t), R(k)(t), l = ˆP(α < L(t)) = 1 5000 5000 k=1 Iα < L(k)_(t)_, and r= ˆP(α > R(t)) = 1 5000 5000 k=1 Iα > R(k)(t),

where I (·) is an indicator function and L(k)_{(t) and R}(k)_{(t) are independent}

copies of L(t) and R(t). We chooseδ = 0.1, thus the ideal value of (l, m, r) is (5%, 90%, 5%). Note that the coverage errors described in Section 2 are the ab-solute difference between (l , m, r) and (5%, 90%, 5%). The experimental results of these measures for various confidence intervals are reported in Tables I–IV. A few interesting things were observed as follows:

(1) The coverage quality of confidence intervals is dominated by the quality of TAVC estimator. For example, in Model 1 and 2, confidence intervals formed by usingα1(t) and vi(t) have very bad coverage quality. Replacing

α1(t) by unbiased estimators improves the cover quality only a bit. However,

replacing vi(t) byσ2has dramatic improvement. This is as expected, since

bias has only a second-order effect on the quality of confidence intervals, see Eq. (5).

(2) The distribution ofα1(t) plays an important role in coverage quality. This

can be seen from the following observations. Confidence intervals con-structed from unbiased estimators and σ2 _{are expected to have (l , m, r)}

close to (5%, 90%, 5%). This is observed in Models 3 and 4, but not for Mod-els 1 and 2. The problem in ModMod-els 1 and 2 is the highly skewed distribution ofα1(t). The skewed distribution ofα1(t) also leads to a high median APE.

(3) Confidence intervals constructed from Var(αi) have comparable coverage

quality of those formed by usingσ2_.

Table I. l , m, and r for Various Confidence Intervals in Model 1 (Displayed in Cells Below TAVC Estimator.) Estimator ¯t σ2 _{t· Var(α} i) α1(t) 1347.3 0.00% 99.50%± 0.16% 0.50% 1.44% 94.62%± 0.52% 3.94% α2(t) 1699.4 0.00% 98.98%± 0.23% 1.02% 0.00% 94.28%± 0.54% 5.72% α3(t) 1347.3 0.00% 99.52%± 0.16% 0.48% 0.98% 95.08%± 0.50% 3.94% α4(t) 1347.3 0.00% 99.54%± 0.16% 0.46% 4.72% 91.28%± 0.66% 4.00% α5(t) 1347.3 0.00% 99.52%± 0.16% 0.48% 1.26% 94.86%± 0.51% 3.88% α6(t) 1347.3 0.00% 99.12%± 0.22% 0.88% 0.00% 93.38%± 0.58% 6.62% α7(t) 1347.3 0.00% 96.34%± 0.44% 3.66% 0.00% 93.76%± 0.56% 6.24% α1(t) 5389.3 0.00% 96.94%± 0.40% 3.06% 0.02% 94.98%± 0.51% 5.00% α2(t) 5784.1 0.00% 95.64%± 0.48% 4.36% 0.00% 94.16%± 0.55% 5.84% α3(t) 5389.3 0.00% 96.00%± 0.46% 4.00% 0.00% 94.16%± 0.55% 5.84% α4(t) 5389.3 0.00% 95.72%± 0.47% 4.28% 0.02% 95.44%± 0.49% 4.54% α5(t) 5389.3 0.00% 96.48%± 0.43% 3.52% 0.00% 94.60%± 0.53% 5.40% α6(t) 5389.3 0.00% 96.50%± 0.43% 3.50% 0.00% 93.90%± 0.56% 6.10% α7(t) 5389.3 0.00% 94.88%± 0.51% 5.12% 0.00% 93.84%± 0.56% 6.16% α1(t) 21557.1 0.16% 94.52%± 0.53% 5.32% 0.36% 93.76%± 0.56% 5.88% α2(t) 21917.3 0.12% 93.92%± 0.56% 5.96% 0.20% 93.56%± 0.57% 6.24% α3(t) 21557.1 0.10% 93.00%± 0.59% 6.90% 0.02% 93.68%± 0.57% 6.30% α4(t) 21557.1 0.16% 93.22%± 0.58% 6.62% 0.02% 94.04%± 0.55% 5.94% α5(t) 21557.1 0.10% 93.96%± 0.55% 5.94% 0.16% 93.80%± 0.56% 6.04% α6(t) 21557.1 0.02% 94.20%± 0.54% 5.78% 0.24% 93.50%± 0.57% 6.26% α7(t) 21557.1 0.10% 94.20%± 0.54% 5.70% 0.22% 93.68%± 0.57% 6.10% Estimator ¯t vi(t) t· vJ(t)/N(t) α1(t) 1347.3 65.16% 25.46%± 1.01% 9.38% 54.34% 36.88%± 1.12% 8.78% α2(t) 1699.4 57.80% 27.84%± 1.04% 14.36% 47.40% 38.38% ± 1.13% 14.22% α3(t) 1347.3 64.44% 24.14%± 1.00% 11.42% 52.86% 38.60% ± 1.13% 8.54% α4(t) 1347.3 76.68% 14.46%± 0.82% 8.86% 65.38% 34.04%± 1.10% 0.58% α5(t) 1347.3 64.86% 24.98%± 1.01% 10.16% 53.64% 37.70% ± 1.13% 8.66% α6(t) 1347.3 47.98% 33.68%± 1.10% 18.34% 35.60% 51.94% ± 1.16% 12.46% α1(t) 5389.3 43.48% 47.54%± 1.16% 8.98% 37.46% 58.54%± 1.15% 4.00% α2(t) 5784.1 39.76% 47.06%± 1.16% 13.18% 33.88% 58.88% ± 1.14% 7.24% α3(t) 5389.3 41.54% 45.50%± 1.16% 12.96% 35.82% 60.16% ± 1.14% 4.02% α4(t) 5389.3 47.40% 41.46%± 1.15% 11.14% 41.40% 57.70% ± 1.15% 0.90% α5(t) 5389.3 42.54% 46.62%± 1.16% 10.84% 36.82% 59.18% ± 1.14% 4.00% α6(t) 5389.3 37.92% 49.82%± 1.16% 12.26% 31.00% 64.22% ± 1.12% 4.78% α1(t) 21557.1 26.66% 70.22%± 1.06% 3.12% 24.50% 73.84%± 1.02% 1.66% α2(t) 21917.3 25.28% 70.00%± 1.07% 4.72% 22.96% 74.52%± 1.01% 2.52% α3(t) 21557.1 25.26% 69.26%± 1.07% 5.48% 23.08% 74.74%± 1.01% 2.18% α4(t) 21557.1 27.04% 68.40%± 1.08% 4.56% 24.80% 73.90%± 1.02% 1.30% α5(t) 21557.1 26.00% 69.98%± 1.07% 4.02% 23.76% 74.28%± 1.02% 1.96% α6(t) 21557.1 24.38% 71.82%± 1.05% 3.80% 22.34% 75.74%± 1.00% 1.92%

Note that m is reported in terms of a 95% confidence interval and ¯t is the average simulated time for each steady-state mean estimator.

(4) Confidence intervals constructed from vJ(t) have consistently better

cover-age quality than those constructed from vi(t).

(5) Confidence intervals constructed from low-bias estimators have similar cov-erage quality to those constructed from unbiased estimators, and have bet-ter coverage quality to those constructed fromα1(t).

Empirical Performance of Bias-Reducing Estimators • 339

Table II. l , m, and r for Various Confidence Intervals in Model 2 (Displayed in Cells Below TAVC Estimator.) Estimator ¯t σ2 _{t· Var(α} i) α1(t) 1347.4 0.00% 99.34%± 0.19% 0.66% 0.52% 94.70%± 0.52% 4.78% α2(t) 1613.4 0.00% 98.78%± 0.26% 1.22% 0.00% 94.08%± 0.55% 5.92% α3(t) 1347.4 0.00% 99.30%± 0.19% 0.70% 0.18% 94.70%± 0.52% 5.12% α4(t) 1347.4 0.00% 98.28%± 0.30% 1.72% 0.86% 94.80%± 0.52% 4.34% α5(t) 1347.4 0.00% 99.30%± 0.19% 0.70% 0.28% 94.70%± 0.52% 5.02% α6(t) 1347.4 0.00% 98.96%± 0.24% 1.04% 0.00% 93.06%± 0.59% 6.94% α7(t) 1347.4 0.00% 95.96%± 0.46% 4.04% 0.00% 93.86%± 0.56% 6.14% α1(t) 5389.5 0.00% 95.94%± 0.46% 4.06% 0.00% 94.24%± 0.54% 5.76% α2(t) 5676.1 0.00% 94.88%± 0.51% 5.12% 0.00% 93.72%± 0.56% 6.28% α3(t) 5389.5 0.00% 94.32%± 0.54% 5.68% 0.00% 93.30%± 0.58% 6.70% α4(t) 5389.5 0.02% 94.10%± 0.55% 5.88% 0.02% 94.24%± 0.54% 5.74% α5(t) 5389.5 0.00% 95.12%± 0.50% 4.88% 0.00% 93.72%± 0.56% 6.28% α6(t) 5389.5 0.00% 95.30%± 0.49% 4.70% 0.00% 93.32%± 0.58% 6.68% α7(t) 5389.5 0.00% 94.98%± 0.51% 5.02% 0.00% 94.74%± 0.52% 5.26% α1(t) 21557.9 0.38% 94.86%± 0.51% 4.76% 0.98% 93.56%± 0.57% 5.46% α2(t) 21845.6 0.34% 94.20%± 0.54% 5.46% 0.48% 93.70%± 0.57% 5.82% α3(t) 21557.9 0.32% 93.50%± 0.57% 6.18% 0.26% 93.70%± 0.57% 6.04% α4(t) 21557.9 0.36% 94.00%± 0.55% 5.64% 0.26% 94.34%± 0.54% 5.40% α5(t) 21557.9 0.36% 94.22%± 0.54% 5.42% 0.40% 93.90%± 0.56% 5.70% α6(t) 21557.9 0.20% 94.62%± 0.52% 5.18% 0.50% 93.60%± 0.57% 5.90% α7(t) 21557.9 0.38% 93.44%± 0.58% 6.18% 0.50% 93.08%± 0.59% 6.42% Estimator ¯t vi(t) t· vJ(t)/N(t) α1(t) 1347.4 60.70% 29.28%± 1.06% 10.02% 51.38% 41.12% ± 1.14% 7.50% α2(t) 1613.4 53.86% 30.80%± 1.07% 15.34% 45.12% 42.52% ± 1.15% 12.36% α3(t) 1347.4 59.50% 27.58%± 1.04% 12.92% 49.54% 43.18% ± 1.15% 7.28% α4(t) 1347.4 70.36% 19.28%± 0.92% 10.36% 59.86% 39.82% ± 1.14% 0.32% α5(t) 1347.4 60.24% 28.40%± 1.05% 11.36% 50.48% 42.16% ± 1.15% 7.36% α6(t) 1347.4 47.38% 35.06%± 1.11% 17.56% 34.94% 54.62% ± 1.16% 10.44% α1(t) 5389.5 37.24% 54.92%± 1.16% 7.84% 33.22% 63.64%± 1.12% 3.14% α2(t) 5676.1 34.36% 54.20%± 1.16% 11.44% 30.50% 63.92% ± 1.12% 5.58% α3(t) 5389.5 36.04% 51.80%± 1.16% 12.16% 31.28% 65.16% ± 1.11% 3.56% α4(t) 5389.5 40.50% 49.16%± 1.16% 10.34% 35.56% 63.22% ± 1.12% 1.22% α5(t) 5389.5 36.64% 53.62%± 1.16% 9.74% 32.34% 64.36%± 1.11% 3.30% α6(t) 5389.5 33.62% 56.40%± 1.15% 9.98% 28.24% 68.06%± 1.08% 3.70% α1(t) 21557.9 24.90% 72.18%± 1.04% 2.92% 23.06% 75.32%± 1.00% 1.62% α2(t) 21845.6 23.68% 71.92%± 1.05% 4.40% 21.68% 75.82%± 1.00% 2.50% α3(t) 21557.9 23.72% 71.80%± 1.05% 4.48% 21.64% 76.20%± 0.99% 2.16% α4(t) 21557.9 25.30% 71.58%± 1.05% 3.12% 23.28% 75.74%± 1.00% 0.98% α5(t) 21557.9 24.26% 72.10%± 1.04% 3.64% 22.38% 75.70%± 1.00% 1.92% α6(t) 21557.9 23.30% 73.32%± 1.03% 3.38% 20.96% 77.12%± 0.98% 1.92%

Note that m is reported in terms of a 95% confidence interval and ¯t is the average simulated time for each steady-state mean estimator.

5. CONCLUSIONS

We have studied seven estimation methods for four stochastic models. For the four low-bias estimators we have studied, we found

(1) In terms of mean APE:α2(t) andα3(t) have superior performance to others.

The performance ofα5(t) is generally better than that ofα1(t), but on average ACM Transactions on Modeling and Computer Simulation, Vol. 14, No. 4, October 2004.

Table III. l , m, and r for Various Confidence Intervals in Model 3 (Displayed in Cells Below TAVC Estimator.) Estimator ¯t σ2 _{t· Var(α} i) α1(t) 589.7 6.46% 90.80%± 0.67% 2.74% 8.22% 88.30%± 0.75% 3.48% α2(t) 814.7 3.16% 94.84%± 0.51% 2.00% 6.56% 89.60%± 0.71% 3.84% α3(t) 589.7 5.16% 91.58%± 0.65% 3.26% 6.50% 89.44%± 0.72% 4.06% α4(t) 589.7 19.90% 60.78%± 1.14% 19.32% 4.16% 88.66%± 0.74% 7.18% α5(t) 589.7 5.66% 91.38%± 0.65% 2.96% 7.40% 88.76%± 0.73% 3.84% α6(t) 589.7 3.06% 92.02%± 0.63% 4.92% 4.12% 89.94%± 0.70% 5.94% α7(t) 589.7 3.46% 91.42%± 0.65% 5.12% 4.22% 90.08%± 0.70% 5.70% α1(t) 1179.3 5.50% 91.62%± 0.64% 2.88% 7.02% 89.26%± 0.72% 3.72% α2(t) 1399.6 3.54% 93.56%± 0.57% 2.90% 5.42% 90.08%± 0.70% 4.50% α3(t) 1179.3 4.26% 91.34%± 0.65% 4.40% 4.94% 89.96%± 0.70% 5.10% α4(t) 1179.3 9.42% 75.88%± 1.00% 14.70% 1.40% 92.24%± 0.62% 6.36% α5(t) 1179.3 4.90% 91.54%± 0.65% 3.56% 5.86% 89.90%± 0.70% 4.24% α6(t) 1179.3 3.56% 92.04%± 0.63% 4.40% 4.58% 89.98%± 0.70% 5.44% α7(t) 1179.3 4.22% 90.44%± 0.68% 5.34% 4.38% 90.08%± 0.70% 5.54% α1(t) 2358.7 5.62% 90.92%± 0.67% 3.46% 6.14% 90.18%± 0.69% 3.68% α2(t) 2573.9 4.10% 92.30%± 0.62% 3.60% 5.36% 90.40%± 0.69% 4.24% α3(t) 2358.7 4.24% 90.52%± 0.68% 5.24% 4.22% 90.62%± 0.68% 5.16% α4(t) 2358.7 6.00% 85.74%± 0.81% 8.26% 2.92% 91.70%± 0.64% 5.38% α5(t) 2358.7 4.96% 90.80%± 0.67% 4.24% 5.10% 90.52%± 0.68% 4.38% α6(t) 2358.7 4.30% 91.10%± 0.66% 4.60% 4.58% 90.46%± 0.68% 4.96% α7(t) 2358.7 4.08% 90.86%± 0.67% 5.06% 4.64% 89.86%± 0.70% 5.50% Estimator ¯t vi(t) t· vJ(t)/N(t) α1(t) 589.7 29.86% 55.82%± 1.16% 14.32% 16.08% 74.58% ± 1.01% 9.34% α2(t) 814.7 21.12% 64.14%± 1.12% 14.74% 11.18% 78.46% ± 0.96% 10.36% α3(t) 589.7 29.72% 52.14%± 1.16% 18.14% 15.70% 74.96% ± 1.01% 9.34% α4(t) 589.7 39.30% 38.62%± 1.13% 22.08% 26.84% 67.02% ± 1.09% 6.14% α5(t) 589.7 29.70% 54.26%± 1.16% 16.04% 15.86% 74.82% ± 1.01% 9.32% α6(t) 589.7 23.36% 57.00%± 1.15% 19.64% 11.80% 77.08% ± 0.98% 11.12% α1(t) 1179.3 18.62% 70.72%± 1.06% 10.66% 9.98% 87.06%± 0.78% 2.96% α2(t) 1399.6 14.24% 73.70%± 1.02% 12.06% 7.22% 89.60%± 0.71% 3.18% α3(t) 1179.3 17.22% 67.90%± 1.09% 14.88% 9.00% 87.76%± 0.76% 3.24% α4(t) 1179.3 22.82% 53.78%± 1.16% 23.40% 12.90% 83.42% ± 0.87% 3.68% α5(t) 1179.3 18.04% 69.92%± 1.07% 12.04% 9.48% 87.44%± 0.77% 3.08% α6(t) 1179.3 15.24% 70.20%± 1.06% 14.56% 7.64% 88.48%± 0.74% 3.88% α1(t) 2358.7 13.58% 78.40%± 0.96% 8.02% 8.44% 88.62%± 0.74% 2.94% α2(t) 2573.9 11.18% 79.56%± 0.94% 9.26% 6.94% 89.80%± 0.70% 3.26% α3(t) 2358.7 11.80% 76.18%± 0.99% 12.02% 7.08% 89.36%± 0.72% 3.56% α4(t) 2358.7 14.50% 69.90%± 1.07% 15.60% 8.50% 87.22%± 0.78% 4.28% α5(t) 2358.7 12.66% 77.76%± 0.97% 9.58% 7.64% 89.08%± 0.73% 3.28% α6(t) 2358.7 11.48% 77.84%± 0.97% 10.68% 6.80% 89.56%± 0.71% 3.64%

Note that m is reported in terms of a 95% confidence interval and ¯t is the average simulated time for each steady-state mean estimator.

is not as good as that ofα2(t) andα3(t). The estimatorα4(t) is unstable for

a few simulation runs because it only uses sample information up to time

T (N (t)) which can be significant shorter than t. To sum up,α2(t) andα3(t)

are the best choices.

(2) In terms of median APE: the performance ranking of low-bias estimators is about the same of that in terms of mean APE. That is,α2(t) andα3(t) are

Empirical Performance of Bias-Reducing Estimators • 341

Table IV. l , m, and r for Various Confidence Intervals in Model 4 (Displayed in Cells Below TAVC Estimator.) Estimator ¯t σ2 _{t· Var(α} i) α1(t) 62.7 1.08% 83.96%± 0.85% 14.96% 1.06% 84.12%± 0.85% 14.82% α2(t) 87.9 1.84% 94.44%± 0.53% 3.72% 3.60% 90.10%± 0.69% 6.30% α3(t) 62.7 1.68% 86.00%± 0.81% 12.32% 1.56% 86.66%± 0.79% 11.78% α4(t) 62.7 3.02% 85.70%± 0.81% 11.28% 2.06% 89.32%± 0.72% 8.62% α5(t) 62.7 2.36% 90.18%± 0.69% 7.46% 2.48% 89.54%± 0.71% 7.98% α6(t) 62.7 4.88% 89.86%± 0.70% 5.26% 4.76% 90.08%± 0.70% 5.16% α7(t) 62.7 5.22% 88.94%± 0.73% 5.84% 4.56% 90.16%± 0.69% 5.28% α1(t) 125.4 1.70% 86.66%± 0.79% 11.64% 1.70% 86.66%± 0.79% 11.64% α2(t) 150.1 2.88% 92.62%± 0.61% 4.50% 4.50% 89.38%± 0.72% 6.12% α3(t) 125.4 2.62% 87.50%± 0.77% 9.88% 2.50% 88.12%± 0.75% 9.38% α4(t) 125.4 7.62% 83.60%± 0.86% 8.78% 4.16% 90.94%± 0.67% 4.90% α5(t) 125.4 3.74% 89.94%± 0.70% 6.32% 3.86% 89.44%± 0.72% 6.70% α6(t) 125.4 4.74% 89.78%± 0.70% 5.48% 4.74% 89.80%± 0.70% 5.46% α7(t) 125.4 4.86% 89.38%± 0.72% 5.76% 4.40% 90.24%± 0.69% 5.36% α1(t) 250.8 2.64% 88.60%± 0.74% 8.76% 2.64% 88.54%± 0.74% 8.82% α2(t) 276.1 4.00% 91.46%± 0.65% 4.54% 4.78% 89.80%± 0.70% 5.42% α3(t) 250.8 3.68% 89.00%± 0.73% 7.32% 3.38% 89.60%± 0.71% 7.02% α4(t) 250.8 6.74% 86.34%± 0.80% 6.92% 4.88% 90.14%± 0.69% 4.98% α5(t) 250.8 4.70% 89.78%± 0.70% 5.52% 4.70% 89.78%± 0.70% 5.52% α6(t) 250.8 4.98% 89.86%± 0.70% 5.16% 5.00% 89.82%± 0.70% 5.18% α7(t) 250.8 5.02% 89.76%± 0.71% 5.22% 4.88% 90.08%± 0.70% 5.04% Estimator ¯t vi(t) t· vJ(t)/N(t) α1(t) 62.7 19.02% 34.40%± 1.11% 46.58% 19.18% 37.00% ± 1.12% 43.82% α2(t) 87.9 30.52% 44.18%± 1.16% 25.30% 30.36% 45.60% ± 1.16% 24.04% α3(t) 62.7 22.80% 33.26%± 1.10% 43.94% 22.00% 38.86% ± 1.13% 39.14% α4(t) 62.7 26.70% 27.92%± 1.04% 45.38% 22.18% 39.24% ± 1.14% 38.58% α5(t) 62.7 26.78% 37.04%± 1.12% 36.18% 25.40% 42.30% ± 1.15% 32.30% α6(t) 62.7 33.24% 40.02%± 1.14% 26.74% 31.98% 43.54% ± 1.15% 24.48% α1(t) 125.4 9.74% 68.38%± 1.08% 21.88% 8.10% 73.68%± 1.02% 18.22% α2(t) 150.1 14.72% 73.54%± 1.03% 11.74% 11.74% 78.94% ± 0.95% 9.32% α3(t) 125.4 12.82% 66.26%± 1.10% 20.92% 9.36% 74.48%± 1.01% 16.16% α4(t) 125.4 21.74% 57.22%± 1.15% 21.04% 12.26% 73.50% ± 1.03% 14.24% α5(t) 125.4 16.00% 67.96%± 1.09% 16.04% 11.46% 76.62% ± 0.98% 11.92% α6(t) 125.4 17.50% 68.52%± 1.08% 13.98% 12.76% 77.36% ± 0.97% 9.88% α1(t) 250.8 6.86% 78.70%± 0.95% 14.44% 4.86% 83.56%± 0.86% 11.58% α2(t) 276.1 10.08% 81.32%± 0.91% 8.60% 6.80% 86.78%± 0.79% 6.42% α3(t) 250.8 8.82% 77.64%± 0.97% 13.54% 5.34% 84.34%± 0.85% 10.32% α4(t) 250.8 13.58% 74.04%± 1.02% 12.38% 8.16% 83.46%± 0.86% 8.38% α5(t) 250.8 11.12% 78.56%± 0.95% 10.32% 7.28% 85.38%± 0.82% 7.34% α6(t) 250.8 11.62% 78.84%± 0.95% 9.54% 7.72% 85.62%± 0.82% 6.66%

Note that m is reported in terms of a 95% confidence interval and ¯t is the average simulated time for each steady-state mean estimator.

also the best choices. However, for some models, the improvement can be quite limited.

(3) In term of computational effort:α2(t) has the drawback of using longer

sim-ulated time, since it require (N (t)+1)th regeneration cycle to complete. One might think the extra simulated time is negligible, but it is not true for some situations. For example, the extra simulated time T (N (t)+ 1) − t ≈ 15Eτ1 ACM Transactions on Modeling and Computer Simulation, Vol. 14, No. 4, October 2004.

in Model 1. The estimatorα5(t) has the drawback of needing extra computer

memory for storing more detailed sample path information for computing

Ai. The work for computing α3(t) andα4(t) is negligible, since the needed

path information is also needed by TAVC estimators and thus the extra work is fairly minimal.

For the TAVC estimators, the jackknife variance estimator vJ(t) generally

gives more accurate coverage than the classical variance estimator vi(t); and

both estimators use about the same amount of computer time (which is negli-gible). However, it must be noted that vJ(t) and vi(t) both underestimated the

TAVC and thus produced confidence intervals that under covered for the short simulation runs (unlike the second TAVC estimator which does quite well but requires much more computation due to replications).

From what has been said above, we may conclude that the confidence inter-val constructed fromα3(t) and vJ(t) is our recommendation for estimating the

steady-state meanα in a regenerative steady-state simulation if the simulated time is sufficiently long. In addition, simulationists may wish to use α3(t) in

place ofα1(t), since from the experimental results, α3(t) is consistantly better

thanα1(t) in terms of mean APE and median APE, and the extra computational

work in computing α3(t) is negligible. Although it may be risky to extrapolate

the recommendations from experimental results of these four models to gen-eral discrete-event simulations, we feel these experimental experiences will be useful to other simulationists.

ACKNOWLEDGMENTS

We wish to thank the editors and the referees for suggestions that improved the exposition.

REFERENCES

CASH, C. R., DIPPOLD, D., LONG, J., NELSON, B., ANDPOLLARD, W. 1992. Evaluation of tests for initial conditions bias. In Proceedings of the 1992 Winter Simulation Conference. IEEE Computer Society, Press, Los Alamitos, Calif., 577–585.

GLYNN, P. W. 1984. Some asymptotic formulas for markov chain with applications to simulation.

J. Stat. Comput. Simul. 19, 97–112.

GLYNN, P. W. 1989. A GSMP formalism for discrete event systems. Proc. IEEE 77. 14–23. GLYNN, P. W. 1994. Some topics in regenerative steady-state simulation. Acta Appl. Math. 34,

225–236.

GLYNN, P. W. 1995. Some new results on the initial transient problem. In Proceedings of the 1995 Winter Simulation Conference. IEEE Computer Society Press, Los Alamitos, Calif. (Arlington,

Va.). 165–170.

GLYNN, P. W.ANDHEIDELBERGER, P. 1990. Bias properties of budget constrained Monte Carlo

sim-ulations. Oper. Res. 38, 801–814.

GLYNN, P. W.ANDHEIDELBERGER, P. 1991a. Analysis of initial transient deletion for replicated

steady-state simulations. Oper. Res. Lett. 10, 437–443.

GLYNN, P. W.ANDHEIDELBERGER, P. 1991b. Analysis of parallel, replicated simulations under a

completion time constraint. ACM Trans. Model. Comput. Simul. 1, 3–23.

GLYNN, P. W.ANDHEIDELBERGER, P. 1992a. Analysis of initial transient deletion for parallel

steady-state simulations. SIAM J. Sci. Stat. Comput. 13, 909–922.

GLYNN, P. W.ANDHEIDELBERGER, P. 1992b. Jackknifing under a budget constraint. ORSA J.

Empirical Performance of Bias-Reducing Estimators • 343

GOLDSMAN, D., SCHRUBEN, L.,ANDSWAIN., J. 1994. Tests for transient means in simulated time

series. Naval Res. Logist. Quart. 41, 171–187.

HENDERSON, S. G.ANDGLYNN, P. W. 2001. Regenerative steady-state simulation of discrete-event

systems. ACM Trans. Model. Comput. Simul. 11, 313–345.

IGLEHART, D. L. 1975. Simulating stable stochastic systems, V: Comparison of ratio estimators.

Naval Res. Logist. Quart. 22, 553–565.

MEKETON, M.ANDHEIDELBERGER, P. 1982. A renewal theoretic approach to bias reduction in re-generative simulations. Manage. Sci. 26, 173–181.

NELSON, B. 1992. Initial-condition bias. In Handbook of Industrial Engineering, 2 ed., G. Salvendy, Ed. Wiley, New York.

SCHRUBEN, L. 1982. Detecting initialization bias in simulation output. Oper. Res. 30, 3, 151–153.

SCHRUBEN, L., SINGH, H.,ANDTIERNEY, L. 1983. Optimal tests for initialization bias in simulation

output. Oper. Res. 31, 6, 1167–1178.

WHITE, K. P., J. 1997. An effective truncation heuristic for bias reduction in simulation output. Simulation 69, 6, 323–334.

WHITE, K. P., J., COBB, M. J.,ANDSPRATT, S. C. 2000. A comparison of five steady-state

trunca-tion heuristics for simulatrunca-tion. In Proceedings of the 2000 Winter Simulatrunca-tion Conference. IEEE Computer Society Press, Loss Alamitos, Calif., 755–760.

Received January 2003; accepted January 2004