Hostname: page-component-74656f9884-8r8nv Total loading time: 0 Render date: 2026-11-08T06:44:24.152Z Has data issue: false hasContentIssue false
  • English
  • Français

NEW ESTIMATORS FOR EFFICIENT GI/G/1 SIMULATION

Published online by Cambridge University Press:  23 March 2005

Chia-Li Wang  and
Ronald W. Wolff
Chia-Li Wang
Affiliation:
Department of Applied Mathematics, National Dong Hwa University, Hualien, Taiwan, Republic of China, E-mail: cwang@mail.ndhu.edu.tw
Ronald W. Wolff
Affiliation:
Department of Industrial Engineering and Operations Research, University of California, Berkeley, California 94720, E-mail: wolff@ieor.berkeley.edu
Rights & Permissions [Opens in a new window]

Abstract

For simulating GI/G/1 queues, we investigate estimators of stationary delay-in-queue moments that were suggested but not investigated in our recent article and we develop new ones that are even more efficient. Among them are direct spread estimators that are functions of a generated sequence of spread idle periods and are combinations of estimators. We also develop corresponding conditional estimators of equilibrium idle-period moments and delay moments. We show that conditional estimators are the most efficient; in fact, for Poisson arrivals, they are exact. In simulation runs with both Erlang and hyperexponential arrivals, conditional estimators of mean delay are more efficient than a published method that estimates idle-period moments by factors well over 100 and by factors of over 800 to several thousand for estimating stationary delay variance.

Information

Type
Research Article
Information
Probability in the Engineering and Informational Sciences , Volume 19 , Issue 2 , April 2005 , pp. 219 - 239
Copyright
© 2005 Cambridge University Press

1. INTRODUCTION

For simulating GI/G/1 queues, we compare the efficiency of estimators of moments of stationary delay in queue D, including standard regenerative estimators that we call direct delay (DD), which are notoriously inefficient in heavy traffic. A method in Minh and Sorli [7], which we call direct idle (DI), estimates average delay d from a generated sequence of idle periods. A method in Wang and Wolff [8] (abbreviated here as WW), called direct equilibrium (DE), estimates delay moments from a generated sequence of equilibrium idle periods. In this article, we continue our investigation of delay moment estimators that are functions of (any kind of) idle periods and we call any estimator of this type an idle-period estimator (IPE).

An alternative approach to simulating GI/G/1 queues, based on conjugate random walks, is briefly described in WW. For a fuller description, see Asmussen [1]; for research results, see Asmussen [2,3] and Ensor and Glynn [5].

We will compare the relative (statistical) efficiency of some estimator a with an alternative estimator b of the same quantity. The statistical efficiency of a (over b) is the factor by which the asymptotic variance of a is smaller than the asymptotic variance of b, for runs of the same length (number of customers or cycles). We use the regenerative method to estimate efficiency factors.

As shown in Section 4.1 of WW, the rate of convergence of IP and DD estimators depends on the relative size of moments of either the interarrival or service distribution. We call distributions with a large coefficient of variation irregular. Distributions irregular in this sense usually have large higher moments as well. It is a rough predictor of the relative efficiency of the estimators that we investigate.

Direct idle estimators are more efficient than DD in heavy traffic. Although not mentioned in Minh and Sorli [7], the efficiency factor of DI over DD is greatly increased when the service distribution is irregular. WW show that the DE estimator of d is twice as efficient as DI for Poisson arrivals (independent and identically distributed (i.i.d.) exponential interarrival times). For a particular irregular interarrival distribution, they obtained an efficiency factor of about 20 for DE over DI. See Section 2.1 for a description of DE estimators.

Of course, computational efficiency also is important. We use the inverse transform method for the distributions considered here. One measure of computational efficiency is to count the relative number of uniform random variables used by different methods. Sometimes, however, this may not reflect actual computation times. In Section 9, we present a table of computational efficiency as ratios of actual computation times for different estimators.

As background, we present well-known results for the GI/G/1 queue. Let Si ∼ G be the service time of arrival i, Ti ∼ A be the time between arrivals i and i + 1, and Xi = Si − Ti, i = 1,2,…, where {Si} and {Ti} are independent sequences of i.i.d. random variables. For generic S, T, and X, let λ = 1/E(T), μ = 1/E(S), and ρ = λ/μ. Assume ρ < 1, E(X2) < ∞, and first-in-first-out (FIFO) order of service. Heavy traffic means that ρ is near 1.

Let Di be the delay in queue of arrival i, the time between that customer's arrival and when service on that customer begins. The Di satisfy the recursion

where a+ ≡ max(a,0). When ρ < 1, the Di in (1) converge in distribution to a proper random variable D, and the process {Di} regenerates whenever an arrival finds the system empty. When we also have E(S2) < ∞, the average delay

and when

, w.p.1, as m → ∞.

Typical performance measures are d and the variance V(D) = E(D2) − d2. From generated Si and Ti values and initial D1, D1 = 0 in this article, suppose we generate (say) m Di values from (1). The DD estimator of E(Dr) is

Suppose that arrivals 1, K1 + 1,K1 + K2 + 1,… find the system empty. Cycles {D1,…,DK1}, {DK1+1,…,DK1+K2},… are i.i.d., as are cycle lengths K1,K2,…. Related regenerative processes may be defined in continuous time, where cycles are called busy cycles. The first busy cycle has length C1 = T1 + ··· + TK1. The first busy period is B1 = S1 + ··· + SK1, and the first idle period is I1 = C1 − B1, where Cj, Bj, and Ij,j ≥ 1, are each i.i.d. sequences. Both E(K1) and E(C1) are finite.

For a fixed number of (say) n cycles, the DD estimator of E(Dr), (3), is of the same form, where now m = Mn is a random variable defined as

An alternative in Minh and Sorli [7] exploits the well-known equation

where I is an idle period. Under their method, Ij, j ≥ 1, are generated using (1). After n cycles, the DI estimators of E(I) and E(I2) are the arithmetic averages of the Ij and the Ij2. To estimate d, substitute these estimators into (4).

As in WW, our IPEs are defined in terms of two distributions.

For Z ∼ F, Z ≥ 0 and E(Z) < ∞, and Fc = 1 − F, let Ze ∼ Fe and Zs ∼ Fs be corresponding equilibrium and spread random variables and distributions, where

, with moments

Because delay moments may be expressed in terms of equilibrium idle-period (Ie) moments, we will compare the efficiency of estimators of the latter quantities.

In this article, we introduce antithetic DE (ADE) and independent DE (IDE) estimators in Section 2, alone and in combination. (DE in WW is denoted ODE here.) We investigate a corresponding family of estimators of Ie moments, called ODS, ADS, and IDS, based on generated values of the spread idle period, Is. For Poisson arrivals, we derive the asymptotic efficiency of DE and DS (direct spread) estimators of d, alone and in combination, in Section 3. In Section 4, we develop IPEs of delay variance V(D), and when arrivals are Poisson, we derive their asymptotic efficiency.

Theoretical and empirical results for Poisson arrivals, where delay moments are known, help validate our simulation runs, and similar results, although not exact, are shown empirically to be good approximations for other interarrival distributions.

In Section 5, we present conditional estimators CE, CI, and CS, corresponding to DE, DI, and DS, respectively, and show that for estimating any moment of Ie, the CE and CI estimators are identical. For E(Ie), the CS estimator also is identical. For Poisson arrivals, conditional estimators are exact. In Section 6, we show that a conditional estimator of E(Ier), r > 0, is more efficient than the corresponding IPE and more efficient than any combination of estimators considered in Section 2.

For Section 6 results, we apply the well-known conditional variance formula

which holds for any random variable Y and class of events

. Variance is reduced by estimating E(Y) from generated samples of

rather than Y.

In Section 7, we report on the relative efficiency of DD, IP, and conditional estimators of d for simulation runs, using combinations of Erlang, exponential, and hyperexponential distributions as interarrival and service distributions. In Section 8, we report on the relative efficiency of V(D) estimators for similar runs.

Combinations of either DE and or DS estimators are more efficient than these estimators used alone. Although expected, we are able to calculate the improvement when arrivals are Poisson. Among the many asymptotic results derived for this case, we consider this to be rather surprising: Each of the estimators ODE, ADE, and IDE (or ODS, ADS, and IDS) of either d or V(D) have the same efficiency when used alone. DS estimators are more efficient than corresponding DE estimators.

The CE (CI) estimators are the most efficient and most useful of those we consider. For both 2-Erlang and hyperexponential arrivals, CE estimators of d are more efficient than corresponding DI estimators by factors well over 100. For similar estimators of V(D), factors ranging from over 800 to several thousand are achieved. Comparing CE with DD, these factors multiply. Thus, CE estimators of delay moments appear to be far more efficient than known alternatives for regular as well as irregular arrivals, and for moderate as well as heavy traffic. In Section 9, we present a method to apply CE estimators to any interarrival distribution.

2. DE AND DS ESTIMATORS OF THE AVERAGE DELAY AND MOMENTS OF THE EQUILIBRIUM IDLE PERIOD

The method developed in WW for generating DE estimators of d and E(Ier),r > 0, also applies to DS estimators. We take up V(D) estimators in Section 4.

Let W = D + S be the stationary waiting time in the system, and let T and Ts be, respectively, an interarrival time and a spread interarrival time, each independent of W. Define U ∼ uniform on (0,1), independent of all else. Idle period I has the representation I = (T − W|T > W).

As discussed in WW and Wolff [9, p.484], Ie has two representations:

where UTs and (1 − U)Ts are ∼Te. The same approach gives a representation of Is, not mentioned in these references, but essential for constructing DS estimators:

2.1. DE Estimators

Because UTs ∼ Te, we write (Te − W|Te > W) ∼ Ie in place of the first representation in (7), where Te and W are independent.

Now, simulate sequences of realized values of Ti and Si over n cycles. From (1), this generates sequences of Di and Wi = Di + Si. For i.i.d. Ui, uniform on (0,1), let Ti = A−1(Ui). In parallel, using the same Ui, generate an i.i.d. sequence of Te values by Tei = Ae−1(Ui). The realized Ie values are the Tei − Wi values on the subsequence where Tei > Wi. Let δi be the indicator of the event Tei > Wi.

This is the setup in WW, where the DE estimator of E(Ier), r > 0, is

Let Rj be the sum of the δi over cycle j, the number of generated Ie values in cycle j. Now, call (9) the original DE (ODE) estimator of E(Ier). For r = 1 and from (4), we get the ODE estimator of d.

Except for Poisson arrivals, the generated Ie values are neither independent nor identically distributed. Nevertheless, the following results were shown in WW.

For every u and v, the limiting fraction of i for which Tei ≤ u and Wi ≤ v is, w.p.1, P(Te ≤ u)P(W ≤ v), where W is a stationary waiting time. In the long run, the generated Ie values are ∼ Ie, in the sense that the limiting fraction of them that take on values ≤t is P(Ie ≤ t) for every t, and as n → ∞,

It is easy to create new DE estimators by generating Ie values in different ways. For antithetic DE (ADE), Ti = A−1(Ui) and Tei = Ae−1(1 − Ui), for the same Ui. For independent DE (IDE), the sequence of Ui used to generate the Ti is independent of the sequence used to generate the Tei. For any IDE sequence, there is also a corresponding antithetic version, AIDE. Each Tei sequence generates a sequence of Ie values in the same way and produces estimators of the same form, the right-hand side of (9). By the argument in WW, (10) holds for them.

At first, we expected ADE and IDE estimators to be less efficient than ODE because (as can be shown) Rj is more variable for them [E(Rj) must be the same for all three]. For Poisson arrivals, however, we find as Theorem 2 in Section 3 that the three estimators have the same asymptotic efficiency.

This suggests using combinations of estimators. For example, generate ODE and ADE quantities on the same run, with ratio estimators of form (9): AODE/BODE and AADE/BADE, respectively. The combination of these estimators, denoted by BDE, is (AODE + AADE)/(BODE + BADE). Let BIDE be the estimator that combines IDE and AIDE, and let CDE be the estimator that combines all four.

Empirically, BDE turns out to be twice as efficient as ODE, more in light traffic. In some cases, CDE is four times as efficient as ODE.

2.2. DS Estimators

Although DS estimators were not mentioned in WW, they are easy to generate, based on (8). We generate i.i.d. sequences Ti = A−1(Ui) and Tsi = As−1(Ui), i ≥ 1. Let Ui′, i ≥ 1, be another i.i.d. sequence, ∼ uniform on (0,1), and independent of all else. The realized Is values are Tsi − Wi, on the subsequence where Ui′Tsi > Wi. Let δi′ be the indicator of the event Ui′Tsi > Wi.

Although new, the estimator denoted as ODS of the rth moment of Is, r > 0,

is of the same form as the ODE estimator, and by the same argument,

As with DE estimators, we will consider a whole family of DS estimators, with the same naming conventions: ADS, IDS, BDS, BIDS, and CDS.

There is one difference, however. It takes two uniform variates, Ui and Ui′, for each arrival i to generate a sequence of Is values. There are three possible antithetic versions: 1 − Ui, Ui′; Ui, 1 − Ui′; and 1 − Ui,1 − Ui′. We use the third version. This matters when original and antithetic versions are combined. The same issue arises even for DE estimators when it is not practical to generate random variables as inverses (e.g., Ti = A−1(Ui)), and more than one uniform variable is used to generate a single Ti or Tei. In our antithetic version, we use 1 − Ui for every Ui.

Because from (5), E(Ier) = E(Isr)/(r + 1), for r > 0, we readily obtain an estimator of E(Ier) from (11). For r = 1, we obtain a DS estimator of d by replacing a DE estimator of E(Ie) by a DS estimator of E(Is), divided by 2. Because dividing by 2 cuts variance by a factor 4, this appears to give DS estimators an advantage. In fact, we show as Theorem 4 in Section 3 that for Poisson arrivals, the ODS estimator of d is twice as efficient as ODE.

3. PROPERTIES OF ESTIMATORS; RESULTS FOR POISSON ARRIVALS

Consider an estimator of some quantity θ,

where {(Nj,Rj)} is an i.i.d. sequence of vectors and θ = E(Nj)/E(Rj). These are the standard assumptions for a ratio estimator. Regenerative estimators are of this form, and we will think of Nj and Rj as generated within some cycle j.

As n → ∞ and when E(Nj2) and E(Rj2) are finite, we have the Regenerative Central Limit Theorem [9, p.122]: For the generic pair N and R,

where

denotes convergence in distribution, and asymptotic variance

is estimated in our simulations by

Suppose generic N and R are of form

where the Yk are realized values of some random quantity Y. Thus,

is simply the arithmetic average of the Yk over n cycles. In most of our simulation runs, the Yk are neither independent nor identically distributed, but have the statistical properties of Y in the long run [as defined for generated Ie values between (9) and (10)], their arithmetic average converges to E(Y); that is, θ = E(Y).

Usually, a regenerative simulation is for a fixed number of full cycles. When the first cycle has a different distribution, usually because the simulation does not begin with a regeneration point at time 0, it is readily shown that (13) still holds.

An alternative is to fix the number of generated Yk to m values, where the estimator of θ is Ym = (Y1 + ··· + Ym)/m. It is well known (e.g., see Asmussen [1, p.136]), that the asymptotic efficiency of this estimator is the same as (13), when normalized for equivalent run length:

In an important special case below, we have Condition

is an i.i.d. sequence, Yk ∼ Y, and R1, R1 + R2,… are stopping times for this sequence, where Y and R (and hence N) have finite second moments.

As the Yk are i.i.d. under

, the standard central limit theorem applies:

Equating asymptotic variances in (17) and (18), we have the following theorem.

Theorem 1: Under Condition

,

This is a remarkable result because the variance of R plays no role. Similarly, aside from the requirement that R be a stopping time, the detailed nature of the dependence between R and {Yk} is irrelevant. Nevertheless, (19) is an elementary consequence of well-known properties of regenerative estimators.

We now derive efficiency results for estimators of d when arrivals are Poisson.

We first compare the efficiency of the ODE, ADE, and IDE estimators of E(Ie) and d, where the generated Yk are realized values of equilibrium idle periods. For ODE,

. For either ADE or IDE, E(R) = 1. For either case, let Tei be the generated equilibrium interarrival time for arrival i. This arrival generates the Yk value Tei − Wi if Tei > Wi. Now, Tei ∼ exp(λ) and is independent of Wi. Owing to the memoryless property of the exponential, the generated Yk value is ∼exp(λ) and is independent of Wi and past history. Thus, the generated Yk values form an i.i.d. sequence. For IDE, R is independent of the entire sequence, so is a stopping time. For ADE, {R = k} and Yk are dependent; nevertheless, R is a stopping time. For all three estimators, (19) holds, and we have the following theorem.

Theorem 2: Under Poisson arrivals, the asymptotic efficiencies of the ODE, ADE, and IDE estimators of d are the same, with asymptotic variance

.

Now, combine IDE with either ODE or ADE. As under Poisson arrivals, the Yk generated by IDE are independent of all else, and we now have E(R) = 2; the asymptotic variance of these runs is

. Thus, the (statistical) efficiency is improved by a factor of 2, compared with either ODE or ADE alone. Of course, using IDE requires generating additional random numbers. This reduces computational efficiency, compared with ODE alone, but by less than a factor of 2.

Combining ODE and ADE into BDE, we know empirically (see Table 1a) that efficiency is improved by a factor of about 2, more in light traffic. (For BDE, no additional random numbers are required.) The improvement in light traffic is easy to understand. In very light traffic, most cycles have only one customer. ODE and ADE each generate one Yk value, where these values are negatively correlated. In fact, as traffic intensity ρ → 0, the correlation is easily calculated:

which gives 11.3 as the efficiency factor by which BDE is better than DI. Apparently, this effect improves the efficiency of BDE in moderate traffic as well. In very light traffic, however, it is likely that DD is more efficient than any of these estimators.

Intriguingly, the empirical efficiencies of BDE and BIDE for these runs are almost identical. Are the true values equal? We believe that they are (slightly) different.

Combining ODE, ADE, IDE, and AIDE into CDE, we get (empirically) an efficiency improvement of 4 or more in Table 1a, compared with one of these estimators alone. In fact, the estimated efficiency factor of CDE is almost exactly the sum of the estimated factors for BDE and BIDE (results for IDE and BIDE are not reported in Table 1a). We now derive a result that partially explains why this is so.

For a generic cycle in (14), let subscript a denote BDE quantities and subscript b denote BIDE quantities. First, we show the following result:

Let

be the σ-field of events generated by all the random variables that determine the queuing system and the generated BDE quantities. It is easy to see that

because for any generated busy period, waiting time values, and BDE values, each generated IDE and AIDE value is ∼ exp(λ), independent of the past, and

is a stopping time for the Yk sequence corresponding to IDE and AIDE. (Actually, Rb is a sum of stopping times, one for the IDE sequence and the other for AIDE.) Thus, (Na − θRa) and (Nb − θRb) are uncorrelated, and (21) holds. Note that these quantities are not independent. As E(Ra) = E(Rb), we get the asymptotic variance of CDE by dividing (21) by 4E(Ra2). We have the following theorem:

Theorem 3: Under Poisson arrivals, the asymptotic variance of the CDE estimator of d is the sum the asymptotic variances of the BDE and BIDE estimators, divided by 4.

Thus, if BDE and BIDE were equally efficient, CDE would be twice as efficient.

The same argument applies to generating additional independent pairs of type IDE and AIDE. Each new pair increases (21) by Var[Nb − θRb]. Of course, these exact results hold only for Poisson arrivals.

We now investigate DS estimators of d and E(Ie) when arrivals are Poisson. In this case, both Ts and Is are ∼2-Erlang, with mean 2/λ and variance 2/λ2.

For ODS (ADS and IDS are similar), we generate sequences {Tsi} and {Ui′Tsi}. For every i where Ui′Tsi > Wi, we generate an Is value:

Now, (1 − Ui′)Tsi and Ui′Tsi are each ∼exp(λ) and they are independent. As (Ui′Tsi − Wi|Ui′Tsi > Wi) is ∼exp(λ), (22) is ∼Is and is independent of Wi and the past. Thus, this procedure generates a sequence of i.i.d. Is values for any of these runs used alone, with stopping times that have the same mean: E(R) = 1. Unlike ODE, however, ODS is not guaranteed to generate exactly one Is value per cycle. As far as efficiency is concerned, this does not matter.

Because we are estimating E(Ie), not E(Is) = 2E(Ie), our estimator is the arithmetic average of the generated Is values, divided by 2, with asymptotic variance V(Ts)/4 = ½λ2. From these results and Theorem 1, we have the following result.

Theorem 4: Under Poisson arrivals, the asymptotic efficiencies of the ODS, ADS, and IDS estimators of d are the same, with asymptotic variance ½λ2. Used alone, they are twice as efficient as the estimators in Theorem 2.

By the same argument used for Theorem 3, we have the following theorem.

Theorem 5: Under Poisson arrivals, the asymptotic variance of the CDS estimator of d is the sum the asymptotic variances of the BDS and BIDS estimators, divided by 4.

4. ESTIMATION OF DELAY VARIANCE

In (10) in WW, V(D) is the sum of known terms, plus variance V(Ie):

We now investigate properties of estimators of V(D) via (23).

Under any of ODE, ADE, or IDE, let Iek be the kth generated value of Ie, k = 1,…,m, with arithmetic average Iem, where m is the number of generated Ie values after n cycles. The DE estimator of V(Ie) is

Before going on, we state a preliminary result. Let {Pm}, {Qm}, and {Zm} be sequences of random variables, F be a distribution function, and

.

Lemma 1:

Although not quite Lemma 1, similar results are proven in Chung [4, p.92] by methods that cover this case.

Let

. Apply Lemma 1 to (24). When convergence occurs,

converge to the same distribution.

In terms of the idle period itself, V(Ie) may be written as in Marshall [6]:

The DI estimator of V(Ie) is

Apply Lemma 1 to (28) with

. When convergence occurs,

converge to the same distribution.

Suppose that arrivals are Poisson. Both the Iek and Ik are i.i.d. ∼ exp(λ), and under ODE, m = n, as Rj = 1. Under either ADE or IDE, E(Rj) = 1. Because the summed quantities on the right in (26) are i.i.d., we have by Theorem 1 that the asymptotic efficiencies under ODE, ADE, and IDE are the same. Now, apply (14) to the right-hand expressions in (26) and (29). We have, after some algebra, the following result.

Theorem 6: Under Poisson arrivals and any of ODE, ADE, or IDE,

Thus, for estimating V(D) under Poisson arrivals, any of the DE estimators, used alone, are

times more efficient than DI.

Under renewal arrivals, we must estimate the asymptotic variance of each variance estimator [e.g., the variances of the Normal distribution in (30) and (31)].

Because the left-hand variance estimators in (26) and (29) are not conventional ratio estimators, we estimate the variances of the right-hand quantities, which are. However, they contain the unknown quantity E(Ie). This situation is similar to estimating (14) by (15). We replace the unknown quantity by an estimator for it,

, that converges to the true value w.p.1.

In both (26) and (29), m is the generated number of either Ie or I values after n cycles. The estimated asymptotic variance of these quantities is given by (15), but to use it, we must break up the run into cycles. For (26), Nj is the sum of terms

over cycle j, and Rj is the number of these terms. Breaking up (29) is similar, where Rj is the sum of generated I values in cycle j. Note that

in (15) is now an estimator of V(Ie), [e.g., (24), (28), or the estimator on the right-hand side of (26) or (29), with

substituted for E(Ie)]. Under DE, let

and estimate the asymptotic variance from

We now show how to estimate V(Ie) from generated Is values, where

From generated values Is1,…,Ism, with mean Ism, our estimator is

Apply Lemma 1 to (32) with Zm = Ism /2E(Ie). When convergence occurs,

converge to the same distribution.

For Poisson arrivals, the Isk values in (33) are i.i.d. ∼ exp(λ) × exp(λ). The asymptotic variance is the variance of (2Isk2 − 3E(Ie)Isk)/6. After some algebra, we get the following result.

Theorem 7: Under Poisson arrivals and any of ODS, ADS, or IDS,

Comparing this result with Theorem 6, any of the DS estimators of V(D) are approximately 2.2 times as efficient as DE and 6.5 times as efficient as DI.

As with DE and DS estimators of d, these estimators of V(D) may be combined. For the next result, we use the argument for Theorem 3, except that now we are dealing with functions of generated quantities. To illustrate for IDE with stopping time R, {Iek,Ie(k+1),…} is independent of {R ≥ k} and, hence, so is the sequence defined by f (Iek) = Iek2 − E(Ie)Iek; see (26). We have the following theorem.

Theorem 8: Under Poisson arrivals, the asymptotic variance of the CDE estimator of V(D) is the sum of the asymptotic variances of the BDE and BIDE estimators, divided by 4. The asymptotic variance of the CDS estimator of V(D) is the sum of the asymptotic variances of the BDS and BIDS estimators, divided by 4.

Under renewal arrivals, estimating the asymptotic variance of DS variance estimators may be done in the same manner described earlier for DE and DI estimators.

5. CONDITIONAL ESTIMATORS OF EQUILIBRIUM-IDLE-PERIOD MOMENTS

For an IPE, there is an opportunity to generate an I, Ie, or Is value for every i. If a value is generated at i, the numerator and denominator of an estimator increase; otherwise, they do not. For a corresponding conditional estimator, the numerator and denominator of this estimator increase for every i by a corresponding expected value, given Wi. We do not generate the actual values, except, of course, for the idle periods at the end of busy cycles. As earlier, our estimators are for n busy cycles.

Corresponding to (9), the DE estimator of E(Ier), we define two functions:

where gE(W)pE(W) = E [(Te − W)+r|W]. The corresponding CE estimator is

For Poisson arrivals, gE(W) = E(Tr), independent of W, and

, a constant. The asymptotic (and actual) variance of the CE estimator is zero! This is a consequence of the memoryless property of the exponential and is true only for this case. For renewal arrivals, we show in Section 6 that CE is better than DE.

Because Te has density ae(x) = λAc(x), we rewrite (36) as

The DI estimator of E(Ier) is constructed from estimated moments of I. The same is true of CI. The CI estimator of E(Ir) is

where from (138) in Wolff [9, p.37]

As E(Ier) = E(Ir+1)/(r + 1)E(I), we have

which is identical to (37). We have the following result.

Lemma 2: The CE and CI estimators of E(Ier), r > 0, are identical.

For the DS estimator (11) of E(Isr), we define

The corresponding CS estimator of E(Ier) is

We now compare the CS and CE estimators. First, observe that as UTs ∼ Te, pS(W) = pE(W). Now, write

Further conditioning this expression on Ts, Ts > W, it is readily seen that [(1 − U)Ts|UTs > W,W,Ts] and [(UTs − W)|UTs > W,W,Ts] are both uniformly distributed on (0,Ts − W). Integrating over the conditional distribution of Ts, given (UTs > W,W), we have

where each is ∼[(Te − W)|Te > W,W]. Combining these results with (36) and (41), we have a result that holds for the first moment.

Lemma 3: The CE, CI, and CS estimators of E(Ie) are identical.

The CS estimators of higher moments may, in general, be different; when different, we expect them to be very difficult to compute. For Poisson arrivals, CS estimators of higher moments are the same. The following is an example.

For Poisson arrivals, the quantities in (42) are ∼exp(λ), independent of W, and they are independent. It is easy to work out

For r = 2, (41) is 2/λ2 = E(T2), a constant, in agreement with (36).

6. LOWER BOUND

We now show that the asymptotic variance of each conditional estimator of E(Ier) is a lower bound on the asymptotic variance of the corresponding IPE and also on all possible combinations of the same estimator, as explained below. We show this for CE estimators, which covers CI. For CS, the form of the proof is identical.

Over cycle j in (9), let

Then rewriting (9),

where (Nj,Rj), j ≥ 1, are i.i.d. vectors. Thus, (43) is a standard ratio estimator. As n → ∞ and when E(Nj2) and E(Rj2) are finite,

with asymptotic variance

For each cycle j with the corresponding conditional estimator, we condition on

, which is Kj and the collection of Wi within that cycle. The asymptotic variance of this estimator is

From (6), we have

Each of the estimators ODE, ADE, and IDE is of form (9); they differ only in the way the Tei are generated. The conditional estimator is identical for each of these estimators; hence, (46) holds for each of them.

When these estimators are combined, two or more Tei values are generated for each i. When, for example, ODE and ADE are combined, the estimator is the ratio of sums of the numerators and denominators of the ODE and ADE estimators. Within each cycle, the corresponding CE estimator sums identical terms, so the CE estimator corresponding to the combined estimator is the same. Thus, (46) still holds, where the right-hand variance is for the combined estimator. Now, IDE may be replicated, and this result holds for any number and combination of estimators.

The same argument applies to DS estimators. We have the following result.

Theorem 9: The asymptotic variance of a conditional estimator of E(Ier), whether DE, DS, or DI, is a lower bound on the asymptotic variance of the corresponding IPE used alone, or in any combination with estimators of the same type.

Applying the same approach to the right-hand side of (26), we also have the following.

Theorem 10: The asymptotic variance of the CE (CI) estimator of V(Ie) and V(D) is a lower bound on the corresponding asymptotic variances of the corresponding DE and DI estimators.

7. SIMULATION RESULTS FOR ESTIMATING AVERAGE DELAY

We now report simulation results for estimating expected delay d. The interarrival and service distributions are one of the following: exponential, 2-Erlang, and hyperexponential with balanced means, denoted by M, E2, and BH, respectively. Hyperexponential H is a mixture of exponentials,

and BH is H with β/μ1 = (1 − β)/μ2. The squared coefficient of variation of each of these distributions is 1, ½, and selected by us to be 10, respectively. Simulation run lengths for these comparisons are 106 busy cycles. This is excessive for our primary purpose, which is to compare the CE, DE, and DS estimators with each other and with DI. We made long runs to stabilize the ratios in the second column.

In Tables 1a–1e, the second column reports the (estimated) factor by which DI is more efficient than DD. Each subsequent column reports the factor by which the designated estimator is more efficient than DI. The product of these factors is the factor by which that estimator is more efficient than DD.

Efficiency Improvement for Estimating d of Various Queues

To save space, we omit results for ADE, IDE, BIDE, ADS, IDS, and BIDS runs, here and in Tables 2a–2e, except for their contributions to combinations that are reported. In most cases, the omitted results are nearly identical to corresponding reported results. When there is a measurable difference between ODE and ADE (or ODS and ADS), the corresponding IDE (or IDS) results are in between. For BH arrivals, ODS is somewhat (up to about 10%) more efficient than ADS. For E2 arrivals, ADS is slightly better than ODS, but they are very close. For Poisson arrivals, of course, the efficiencies are the same.

Efficiency Improvement for Estimating V(D) of Various Queues

Table 1a is consistent with Theorems 2–5 and, of course, the fact that the CE estimator has zero variance. BDE is better than ODE and ADE by a factor of 2 or more, with the factor increasing in light traffic. Irregular service reduces the efficiency of DD, but has little effect on IPEs.

In Tables 1d and 1e, CDE (CDS) falls well short of twice the efficiency of BDE (BDS). Ultimately, the CE bound must cause the efficiency increase from generating additional Tei or Tsi values to be subadditive; that is, the combined improvement is less that the sum of the individual improvements. Subadditivity also occurs in Tables 1b and 1c, but the effect is smaller and the bound is further away.

8. SIMULATION RESULTS FOR ESTIMATING DELAY VARIANCE

These results are for the same set of runs and table layout as in Section 7. Tables 2a–2e are consistent with Theorems 6–8. These tables exhibit patterns observed earlier.

9. CONCLUSION

As observed in WW, the advantage DE has over DI for estimating d arises from the fact that variance V(Ie) is a function of only the first three moments of I, whereas the variance of the DI estimator, which requires separate estimation of E(I2), depends on the fourth moment of I. As moments of I are closely related to moments of T, irregular arrivals strongly favor DE over DI. The explanation for DS is similar, and our simulation runs are consistent with these observations. For estimating V(D), DE and DS do even better, relative to DI, for similar reasons. In both cases, DS is more efficient than DE, as our analysis for Poisson arrivals would lead us to expect. In general, however, using DS requires generating more uniform random variables and is less efficient computationally.

For DE or DS estimators used alone, this is the primary source of increased efficiency over DI. It does not matter much how we generate these estimators (ODE, ADE, and so on), and results are only weakly dependent on ρ and the service distribution. Initially, combinations such as BDE and BDS produce substantial and worthwhile efficiency gains. Eventually, however, the CE bound takes over, and the efficiency gains from adding more combinations become subadditive, except in the unique case where arrivals are Poisson.

The CE (or what turned out to be the equivalent CI) estimators not only provide bounds, but they also give achievable efficiency gains that are compelling.

The CI estimators could have been considered much earlier, but no one thought of them. They are certainly much easier to explain than CE, as there is no need to discuss the generation of equilibrium idle periods.

For E2 and BH arrivals, the CE integrations are easy to perform analytically. When analytical integrations are not possible, one can prepare precalculated numerical values of the integrals as functions of waiting time W and then use a hashing function for efficient hash-table lookup.

In Table 3, we report computational efficiency as the ratio of actual running times of our simulation runs. The computation of conditional probability and expectation involves exponential functions, which, apparently, consumes more computer time than the generation of extra uniform random variables. Nevertheless, the computational efficiency of DI over the others is not more than a factor of 2 in Table 3. We believe computational efficiency will be similar for a broader class of phase-type distributions but have not investigated this matter beyond what is presented here.

Computational Efficiency for Estimating d

The interarrival distributions in our runs have infinite tails. We also experimented on one distribution with finite tail, the U/BH/1 queue, where interarrival times are uniformly distributed on (0,2/λ). CE estimators of d are more efficient than DI by factors from 70 to 7, as ρ increases from 0.1 to 0.9. Because this distribution is very regular, ODE does not perform well here. Also note that for deterministic interarrival times, the CI, and thus CE, estimator is exactly the same as DI.

The CE estimators of delay moments appear to be far more efficient statistically than known alternatives for regular as well as irregular arrivals and for moderate as well as heavy traffic. Compared with DD, efficiency factors sometimes are astronomical. Using hash-table lookup, CE applies to any interarrival distribution.

References

REFERENCES

Asmussen, S. (1987). Applied probability and queues. New York: Wiley.
Asmussen, S. (1990). Exponential families and regression in the Monte Carlo study of queues and random walks. Annuals of Statistics 18: 1851–1867.Google Scholar
Asmussen, S. (1992). Queueing simulation in heavy traffic. Mathematics of Operations Research 17: 84–111.Google Scholar
Chung, K.L. (1974). A course in probability theory, 2nd ed. New York: Academic Press.
Ensor, E.A. & Glynn, P.W. (2000). Simulating the maximum of a random walk. Journal of Statistical Planning and Inference 85: 127–135.Google Scholar
Marshall, K.T. (1968). Some relationships between the distributions of waiting time, idle time, and interoutput time in the GI/G/1 queue. SIAM Journal on Applied Mathematics 16: 324–327CrossRefGoogle Scholar
Minh, D.L. & Sorli, R.M. (1983). Simulating the GI/G/1 queue in heavy traffic. Operations Research 31: 966–971.Google Scholar
Wang, C.-L. & Wolff, R.W. (2003). Efficient simulation of queues in heavy traffic. ACM TOMACS 13: 62–81.Google Scholar
Wolff, R.W. (1989). Stochastic modeling and the theory of queues. Englewood Cliffs, NJ: Prentice-Hall.
Figure 0

Efficiency Improvement for Estimating d of Various Queues

Figure 1

Efficiency Improvement for Estimating V(D) of Various Queues

Figure 2

Computational Efficiency for Estimating d