Estimating small cell-loss ratios in ATM switches via importance sampling.
The cell-loss ratio at a given node in an ATM switch, defined as the steady-state fraction of packets of information that are lost at that node due to buffer overflow, is typically a very small quantity that is hard to estimate by simulation. Cell losses are rare events, and importance sampling is sometimes the appropriate tool in this situation. However, finding the right change of measure is generally difficult. In this article, importance sampling is applied to estimate the cell-loss ratio in an ATM switch modeled as a queuing network that is fed by several sources emitting cells according to a Markov-modulated ON/OFF process, and in which all the cells from the same source have the same destination. The charge of measure is obtained via an adaptation of a heuristic proposed by Chang et al. {1994} for intree networks. The numerical experiments confirm important efficiency improvements even for large nonintree networks and a large number of sources. Experiments with different variants of the importance sampling methodology are also reported, and a number of practical issues are illustrated and discussed.
INTRODUCTION
An Asynchronous Transfer Mode (ATM) communication switch can be modeled as
a network of queues with finite bu#er sizes. Packets of information (called cells)
join the network from several sources according to stochastic processes, and some
cells may be lost due to bu#er overflow. The long-term (or steady-state) fraction
of cells that are lost at a given node is called the cell-loss ratio (CLR) at that node.
Typical CLRs are small (e.g., less than 10 -8 ) and the cell losses tend to occur in
bunches. Cell losses are thus so rare that estimating the CLR with good precision
by straightforward simulation is very time-consuming, and in some cases practically
impossible.
Importance sampling (IS) (e.g., [Heidelberger 1995]) and the splitting method
(e.g., [Glasserman et al. 1999]) are two important candidate methods for improving
the quality of the estimator in such a situation. IS changes the probability laws
governing the system so that the rare events of interest occur more frequently,
eventually to the point of being no longer rare events. The estimator is modified
accordingly so that it remains unbiased: It is multiplied by a quantity called the
likelihood ratio. The hope is that the IS estimator is more e#cient in the sense
that the product of its variance and its computing cost is smaller than for the
regular estimator. The most di#cult problem in applying IS is (in general) to
figure out how to change the probability laws so that the variance gets reduced to an
acceptable level. Theoretically, there always exists a change of measure that reduces
the variance to zero, but it is usually much too complicated and impractical, as it
frequently leads to a non-Markovian process [Glynn and Iglehart 1989; Heidelberger
1995].
For general background on e#ciency improvement (or variance reduction), we
refer the reader to Bratley et al. [1987], Fishman [1996], Glynn [1994], and L'Ecuyer
[1994]. IS is well explained in Glynn and Iglehart [1989], Shahabuddin [1994],
Heidelberger [1995], Sadowsky [1996], and the several other references given there.
Application of IS to the simulation of communication systems is studied, e.g., by
Parekh and Walrand [1989], Chang et al. [1994], Chang et al. [1995], Bonneau
[1996], Falkner et al. [1999], and Heegaard [1998].
Chang et al. [1994] derived an asymptotically optimal change of measure, based
on the theories of e#ective bandwidth and large deviations , for estimating the probability
p that a queue length exceeds a given level x before returning to empty, given
that the queue is started from empty, for a single queue with multiple independent
arrival sources that satisfy a large deviation principle. Roughly, asymptotically
optimal means that the standard error of the IS estimator converges to zero exponentially
fast with the same decay rate (exponent) as the quantity to be estimated,
as a function of the level x. A precise definition can be found in Chang et al. [1994]
and Heidelberger [1995]. An asymptotically optimal change of measure does not
minimize the variance, but it can reduce it by several orders of magnitude. Chang
et al. [1994] extended their method to intree networks of queues, which are acyclic
tree networks where customers flow only towards the root of the tree. For intree
networks, assuming infinite bu#ers at all nodes, they obtained a lower bound on
the exponent in the (exponential) convergence rate of the variance of the IS esti-
mator, but they could not prove that this estimator is asymptotically optimal. The
Small Cell Loss Ratios in ATM via IS - 3
missing step for proving the latter is to show that the square mean converges exponentially
to 0 at the same rate (i.e., with the same exponent, and not faster) than
the variance. In a related paper, Chang [1995] obtained the exact exponent in the
convergence rate of P [q # > x] towards 0 as x #, where the random variable q #
is the steady-state queue length in an intree queueing network with infinite bu#ers.
This exponent is linked with exactly the same change of measure as in Chang et
al. [1994]. The probabilities p and P [q # > x] just described are closely related
to the CLR when x equals the bu#er size (it measures similar events), so it seems
reasonable to use the change of measure proposed by Chang et al. [1994] to estimate
the CLR as well. In fact, these authors applied their IS strategy to estimate the
steady-state average number of cell losses per unit of time, another quantity closely
related to the CLR, and observed large empirical variance reductions for examples
of queueing models with a single node and two nodes in series .
In this paper, we pursue these experiments by considering much larger networks,
some of them being non intree. We consider queueing networks having a large
number of nodes, fed by a large number of Markov-modulated ON/OFF sources
(up to 300 in our numerical examples). The nodes are organized in successive layers.
Each cell (or customer) goes through exactly one node of each layer, following a
path uniquely determined by its source. This type of queueing network models the
tra#c in an ATM switch [Dabrowski et al. 1999; Falkner et al. 1999; Giroux and
Ganti 1999]. We apply IS to estimate the CLR at any prespecified node of the
network, using a direct adaptation of the IS methodology of Chang et al. [1994].
The idea is to increase the tra#c to the target node by increasing the average
ON/OFF ratio for all the sources (and only those) feeding that node. The exact
change of measure is determined by a heuristic which turns out to be equivalent to
Algorithm 3.2 of Chang et al. [1994] if we neglect all the tra#c not directed towards
the target node and the corresponding sources, and if we assume that the average
arrival rate after IS does not exceed the service rate at the nodes upstream of the
target node. The latter condition is always satisfied in our examples, but if it was
not, we could simply enforce it by choosing a change of measure that gives an arrival
rate equal to the service rate at those nodes where the condition is violated, exactly
as in Algorithm 3.2 of Chang et al. [1994]. Falkner et al. [1999] have proposed an
additional heuristic, based on the concept of decoupling bandwidth introduced in
de Veciana et al. [1994], to reduce the change of measure when the e#ect of the
additional tra#c not directed towards the target node may be non-negligible.
For both small and large networks, we obtain large e#ciency improvements (em-
pirically), especially when the CLR is small. When the size of the bu#er at the
target node increases, IS seems to improve the e#ciency by a factor that increases
exponentially with the bu#er size. This is consistent with the conjecture that the
change of measure is asymptotically optimal.
It must be emphasized that this IS methodology is heuristic because (i) CLR is
not the same as p; (ii) it is still unproven if the IS estimator proposed in Chang et
al. [1994] is asymptotically optimal for p, and proving it for CLR in the presence
of finite bu#ers at all nodes appears quite messy; (iii) even if the IS estimator was
asymptotically optimal for the CLR, this would only be an asymptotic result. It
would not necessarily imply that the variance is reduced su#ciently and that the
estimator is practically acceptable for the concrete models that one is interested in.
P. L'Ecuyer and Y. Champoux
Actual experimentation with this IS methodology, with examples similar to real-life
problems of interest, is therefore needed to see how practical the method is, and
whether or not additional heuristics can improve on it. The aim of the present
paper is to report on such experimentation.
Our examples also illustrate the fact, typical of rare-event simulation contexts,
that the variance estimators tend to be more noisy (in terms of relative error) than
the mean estimators, even after IS is applied. One must therefore be careful when
comparing empirical variances, or with the interpretation of confidence intervals.
(This important fact is typically left unmentioned in papers that report empirical
results on rare-event simulation.) The noise in the variance estimator is not caused
by the IS methodology, but is due to the fact that the cell losses are rare events. In
fact, in our examples, IS improves the variance estimator by a much larger factor
than the mean estimator. This is along the lines of the results of Sadowsky [1993],
who has shown, in the context of estimating large deviations probabilities for a
sum of independent and identically distributed random variables (a much simpler
situation than our model), that an optimal exponentially twisted change of measure
stabilizes all the moments of the estimator, in addition to being asymptotically
optimal. These results have been generalized to a more abstract setup in Sadowsky
[1996].
Beck et al. [1999] and Dabrowski et al. [1998] also study the application of IS to
a discrete-time queueing network model of an ATM switch. Their model is very
general. They obtain the asymptotics of the tail of the queue size distribution in
steady-state, and they use that to propose a change of measure for estimating the
CLR at a given node. Their IS methodology is related to (but di#erent from) that
proposed by Chang et al. [1994] (see also Section 6.3).
We have limited the generality of the model in order to avoid excessive notation
and unnecessary complexity. The model is flexible enough to illustrate how the IS
strategy behaves in several important situations. In addition to our basic model,
we also made several experiments with IS for generalizations and other variants
of this model. The proposed IS scheme works fine for some variants and fails for
others. This is summarized in Section 5.5. In particular, the fact that each source
is assigned a fixed destination (in contrast to having random routing for each cell
or group of cells) is an important factor for the success of the method.
The geometric distribution of the ON and OFF periods for the sources (implied
by the Markov-modulated model) is certainly not always realistic. Statistical analyses
of broadband tra#c traces (for Ethernet, video, etc.) indicate that the tra#c
is often long-range dependent, and that a more representative model in this case
is obtained when the sojourn-time distributions in the ON/OFF states are heavy-tailed
(e.g., Pareto-type, with infinite variance) [Beran et al. 1995; Leland et al.
1994; Willinger et al. 1995; Neame et al. 1999]. Unfortunately, the IS methodology
used in this paper relies heavily on the exponential tails of the ON/OFF
distributions, and it is not obvious how (and if) it can be adapted to heavy-tailed
distributions. On the other hand, despite the empirical evidence of the limited
representativity of the Markov-modulated model in certain contexts, there remain
supporters of this model who argue, based on experiments, that in networks with
finite bu#ers the correlations of very long range are usually not important to model
anyway [Heyman and Lakshman 1996; Krunz and Makowski 1998]. In any case,
Small Cell Loss Ratios in ATM via IS - 5
we think that the Markov-modulated model is still an important model, and that a
better understanding of the e#ciency improvement methods proposed for it remains
worthwhile.
The model is defined is Section 2. Section 3 recalls the A-cycles method and the
batch-means method, which are used jointly to compute confidence intervals. In
Section 4 we explain how IS is applied to estimate the CLR at a given target node.
The change of measure is the same as in Chang et al. [1994] for intree networks.
However, our explanation of the heuristic di#ers from that of these authors. It is
given directly in terms of the CLR and is along the lines of the arguments given
in Heidelberger [1995], pages 62-63. Numerical results are reported in Section 5.
We also summarize our numerical experiments with several other variants of the
model. In Section 6, we consider various refinements of the basic IS scheme, and
test them empirically to see how much additional variance reduction they can bring.
We report both positive and negative results. Section 7 explains how the CLR can
be estimated in functional form, as a function of certain parameters of the model.
Section 8 gives a conclusion. A preliminary report of this work was presented in
L'Ecuyer and Champoux [1996] and additional numerical results and details can
be found in Champoux [1998].
2. THE MODEL
We consider an acyclic queueing network with 4 layers of nodes , as illustrated in

Figure

1. Each node is a single-server FIFO queue with finite bu#er size. The #-th
layer is called level # and the nodes at level 4 transmit cells to destinations . Levels
2 and 3 have m 2 nodes each, while levels 1 and 4 have nodes each. Each
level-2 node is fed by m 1 level-1 nodes, while each level-3 node feeds m 1 nodes
at level 4. Cells (i.e., packets of information) arrive at level 1, visit one node of
each level, in succession, then leave the network. Each node at level 1 is fed by
sources . These are assigned to specific destinations;
i.e., all the cells produced by a given source follow exactly the same path. The
trajectory of a cell at the switching stage between levels 2 and 3 is thus determined
by its fixed destination. The arrival sources are time-synchronized, but otherwise
independent, stochastically identical, discrete-time ON/OFF Markov modulated
processes. A source is OFF for a while, then ON for a while, then OFF again, and
so on. The source produces one cell per unit of time during an ON period, and none
during an OFF period. The durations of OFF and ON periods are independent
geometric random variables with means # 0 and # 1 , respectively, so the arrival rate
called the average burst size. If we denote
OFF and ON by 0 and 1, respectively, these assumptions imply that the state of
a source evolves as a discrete-time Markov chain with two states, 0 and 1, with
transition probability matrix
. (1)
These Markov chains comprise all the stochasticity of the model; everything else is
deterministic. The arrival sources are numbered from 1 to and the nodes
are numbered from 1 to 2m 2 (1+m 1 ), level by level. When two or more cells reach a
given node simultaneously, they are placed in the queue (the bu#er) by order of the
P. L'Ecuyer and Y. Champoux
number of the node or source where they come from. This deterministic ordering
rule is for simplification and tends to favor the cells coming from certain sources
and nodes. One could order the cells randomly instead, but that would have no
major qualitative impact on our results.
Fig. 1. An ATM Switch Modeled as 4 Layers of Queues with Finite Bu#er Sizes
All the nodes at level # have the same bu#er size B # and the same constant service
time 1/c # (so c # is the service rate). Whenever a cell arrives at a node where the
bu#er is full, it is lost and disappears from the network. Our aim is to estimate
the CLR at a given node of the network, say node q # at level # , where the CLR is
defined as
where NT (#) is the total number of cells reaching node q # during the time interval
(#) is the number of those cells that are lost due to bu#er overflow at
node q # . We assume that the total arrival rate is less than the service rate at each
node. That is, if the cells from m sources pass through a given node at level #, then
m# < c # , and this holds for all nodes.
To simplify the discussion, we assume that each c # is an integer. Since the bu#ers
are finite, the entire model is then a discrete-time Markov chain with finite state
space. It is also aperiodic, and the zero state (the state where all sources are OFF
and all the nodes are empty) is positive recurent and is accessible from every other
state. As a consequence, there exists a limiting distribution # over the state space
of that chain, defined as
Small Cell Loss Ratios in ATM via IS - 7
where # is an integer and S(#) is the state of the chain at time # .
This model could of course be generalized in several directions and our approach
would be easy to adapt for certain types of generalizations (see also Section 5.5).
For example, the bu#er sizes and constant service times can di#er between nodes
at a given level, di#erent sources can have di#erent transition probability matrices
R, a source could produce a cell only with some probability when it is ON, and
additional sources may feed the nodes at levels higher than 1. IS still works in
these situations. We keep our simpler model to avoid burying the key ideas under
a complicated notation. On the other hand, if the destinations were determined
randomly and independently for each cell, or for each ON period at each source,
finding an e#cient way of applying IS would be more di#cult. The fixed source-destination
assignment model is reasonable if we assume that a typical connection
between a source and a destination lasts for several orders of magnitude longer than
the service times 1/c # , and much longer than the time required for bu#er overflow
in the most likely path to overflow. This assumption is commonly made in the
literature [Dabrowski et al. 1999; Falkner et al. 1999; Giroux and Ganti 1999] and
by the switch manufacturers.
In the case where # = 3 or 4 and all the nodes at level # are statistically
identical, an alternative estimator would take the global (empirical) average CLR
for all nodes at level # . With a straightforward simulation approach (no IS), this
would yield a better estimator than concentrating on a single node q # . But with
IS, according to our empirical investigations, it is much better to concentrate on a
single node and increase only the tra#c to that node. With the latter, the likelihood
ratio is less noisy and there is less (undesired) overflow at the nodes of level # .
3. A REGENERATION APPROACH FOR CONFIDENCE INTERVALS
IS is generally easier to apply to a model defined over a short time horizon or when
the model's evolution can be decomposed into short regenerative cycles. Here, the
model is over an infinite horizon, and to decompose its trajectory into cycles, we
apply a generalization of the classical regenerative method introduced in Nicola
et al. [1993] and Chang et al. [1994], and called the A-cycle method. Let A be a
subset of the state space of the system. Here we take A as the set of states for
which the bu#er at q # is empty. Let #
the bu#er at q # is empty at time # but not at time # - 1}. These # i are the
successive hitting times of the set A by the Markov chain {S(# IN}. The pair
is also a positive recurrent aperiodic Markov chain over
a finite state space, so it has a pointwise limit distribution - #, say, as #. Then
the state of the chain at the hitting times # i has the pointwise limit distribution #,
over A, defined by:
where A is the complement of the set A.
The process over the time interval (# i-1 , # i ] is called the ith A-cycle. Let X i be
the number of cells reaching node q # during the ith A-cycle, and Y i be the number
of those X i cells that are lost due to bu#er overflow at q # . Let E # denote the
mathematical expectation over an A-cycle when the initial state (at the beginning
P. L'Ecuyer and Y. Champoux
of the A-cycle) has distribution #. One has:
In the limit, as the number of A-cycles increases, the average distribution of the
system states at the times # i approaches #. By taking the average of the Y i and X i
over the first n A-cycles, one obtains the consistent estimator of -:
. (4)
This estimator is biased because the initial state at time 0 is normally not generated
from # (this would be too di#cult) and because it is a ratio estimator. However,
the bias can be reduced by warming up the system, e.g., by running n 0 +n A-cycles
and discarding the first n 0 from the statistics.
The A-cycles are asymptotically identically distributed (with probability law #
for their initial state) but they are dependent . To reduce the dependence, and
also improve the normality, one can batch the cycles, as in the usual batch means
method. One then applies the standard methodology for computing a confidence
interval for a ratio of expectations, using the batch means as observations, and
obtain a confidence interval for - [Law and Kelton 2000]. In our experiments, we
fix the number of cycles per batch (and not the simulation time per batch), as
explained in Section 5.1.
4. APPLYING IMPORTANCE SAMPLING
When - is very small, the vast majority of the Y i 's in (4) are 0 and the relative error
of -
- (i.e., its standard deviation divided by -) blows up. In (3), the denominator
is easy to estimate, but the numerator is hard to estimate because it depends
on rare events. In fact, denoting - observing that Y 1 is a non-negative
integer, one has Var # [Y 1
Y
the squared relative error satisfies
Y
as Following Chang et al. [1994], we will use IS for the numerator of (3)
and not for the denominator.
denote the set of sources feeding q # . The IS strategy for increasing the
tra#c towards q # is to increase r 01 and r 11 in the matrix R, for all the sources that
belong to S # and only those, so that the total long run arrival rate at q # becomes
larger than the service rate. The system starts with an empty bu#er at q # (a
state in A) and the change remains in e#ect until the bu#er at q # empties again or
overflows. When the bu#er overflows, R is set back to its original for all the sources
until the bu#er at q # empties again, which marks the end of the A-cycle. We call
this an A-cycle with IS . Under this strategy, if the tra#c to q # can be increased
su#ciently, cell losses are no longer rare events. This can certainly be achieved if
is the cardinality of S # and c is the service rate at the
target node.
Small Cell Loss Ratios in ATM via IS - 9
It remains to decide how to change R. For this, we proceed in a standard way
[Heidelberger 1995]. For a real-valued parameter #, define
let #) be the spectral radius (largest eigenvalue) of #), and let (f 0 (#), f 1 (#
be the corresponding right eigenvector (the prime means transposed), so that
The eigenvalue #) can be written explicitly as
For IS, we will change R to the stochastic matrix
#)
This formulation is quite flexible, because the mean arrival rate from a source can
be set to an arbitrary value between 0 and 1 by choosing an appropriate #, and it
leads to important simplifications in the likelihood ratio over an A-cycle, as we will
see.
During a given A-cycle, let N ij be the number of times a source in S # goes from
state i to state j while using the probabilities 1. The
total number of transitions generated from -
R is then N 00 +N
where t is the number of time steps where IS is on. The state transitions of the
sources are assumed to occur right before the (discrete) times of cell production.
The number of cells generated for q # during the time interval (0, t] is thus N 01 +N 11 .
By a simple counting argument, assuming that the bu#er overflows at time t, one
can decompose
is the bu#er size at q # , Q t is the number of cells already generated and
on their way to node q # at time t, Q 0 is the corresponding number at time 0, F t
is the di#erence between the total capacity of service c # t of the server at q # during
(0, t] and the actual number of cells served at q # during that interval of time, and
H t is the number of cells headed to q # but lost due to bu#er overflow either at q #
or upstream during (0, t]. In (6), B # is the number of cells filling up the bu#er
and c # t - F t is the total number of cells served at q # during (0, t]. We denote
The transition probability of a source from state i to state j is r ij without IS,
so each such transition under IS contributes a factor r ij /-r ij to
the likelihood ratio associated with this change of probabilities. Since there are N ij
transitions from i to j, the likelihood ratio becomes
N00
r 01
r 01
r 11
N11
P. L'Ecuyer and Y. Champoux
#)
N11
where
It is well known that if V is a random variable defined over an A-cycle with initial
state that has distribution #, then E # [V
denotes the expectation
under the probabilities -
R, over an A-cycle with IS, with initial state drawn
from #. Thus, computing LY 1 over the A-cycle with IS yields an unbiased estimator
of One always has -
1. IS is e#ective (roughly) if L tends to have a
small variance conditional on Y 1 > 0. In other words, we would like to match the
small values of L with the positive values of Y 1 (i.e., the bu#er overflows) and the
large values of L with the no-overflow situations.
We still have the choice of # at our disposal. Since t is an unbounded random
variable likely to have significant variance, a standard strategy in this situation
(e.g., [Heidelberger 1995]) is to choose # so that t disappears from L, i.e, take
This # gives the same IS strategy as that provided in Chang et al. [1994] if we
consider only the sources that feed q # and if we neglect the possibility that the
arrival rate exceeds the service rate at upstream nodes. The latter condition
was satisfied in all the examples that we have considered, but when it is not, one
can reduce # to meet the condition exactly as in Chang et al. [1994]. Note that
#) is a strictly increasing and di#erentiable function of #
(see, e.g., Chang et al. [1994], Example 2.6), and that # -1
Therefore, this # exists if and only if m # > c # , which we assume (otherwise, one
cannot overload the node q # and overflows should be negligible anyway unless B #
is near 0). With # , the likelihood ratio becomes
The variance of the estimator of - Y is -
Y where
The squared relative error of the IS estimator satisfies
Y
Y
. (10)
The ratio of expectations in (10) is # 1, by the Cauchy-Schwarz inequality. Bounding
this ratio by a constant independent of B # would prove that the relative error
under IS is bounded, but we do not have a proof. Sadowsky [1991] was able to prove
that IS based on an optimal exponential twisting found by solving an equation similar
to (8) is asymptotically optimal for estimating p in a single GI/GI/m queue
Small Cell Loss Ratios in ATM via IS - 11
with finite bu#er. The results of Chang et al. [1994] generalize this to dependent
interarrival times, still for a single queue.
We continue the discussion with intuitive heuristic arguments suggesting that L
is likely to be small when Y 1 > 0, and that the variance of LY 1 is also likely to be
small. Firstly, since q # is stable without IS and since IS is stopped as soon as the
bu#er overflows, Y 1 is unlikely to take very large values (it should be possible to
prove that its distribution has an exponentially bounded tail under IS). Concerning
for the matrices R and the values of m # /c # that correspond to typical
examples (see Section 5), one usually has f 0
it is unlikely that less sources are ON under IS than without IS, because
the ON state has a higher steady state probability under IS, so N 01
extremely unlikely. Therefore, it seems reasonable to expect that W (# ) is less
than 1 with large probability and has a small variance, bounded independently of
. (In fact, W (# ) is bounded independently of B # because |N
by m # , but this bound can be very large unless m # is small.) Now it su#ces to
hope that the large negative values of # are rare enough to be negligible when
This hope appears reasonable for the following reasons. Since the sources
are overproducing cells headed towards q # when IS is ON, and that the initial state
at time 0 follows approximately the limiting distribution # without IS, one should
normally expect Q t # Q 0 , and it should be (intuitively) extremely unlikely that
. (Remember that Q 0 and Q t are the total number of cells in the upstream
nodes, in steady-state (roughly), without IS and under IS, respectively.) Moreover,
since the production rate of the sources exceeds the capacity at q # , the server at
q # is not expected to remain idle very long, so a large positive value of F t is very
unlikely. Since exp(-# B # ) is a constant, the above arguments suggest that L
should tend to be very small when Y 1 > 0. From another viewpoint, when Y
(no overflow), the value of # is likely to be slightly smaller than -B # , i.e., the
total number of cells produced during (0, t] should be slightly less than the total
service capacity c # t, because the bu#er at q # is empty when IS is turned OFF. This
suggests that L should tend to be larger than 1 when Y which is good news
in view of the fact that -
As a very crude heuristic, suppose we assume that e -# Y 1 W (# ) in (9) is
bounded by a constant K 1 independent of B # . This would give
in which case IS would provide the approximate variance reduction factor
Y
Y
Again, these arguments are only hand-waving heuristics and the real test is to see
how the variance is actually reduced on concrete examples. One may be tempted
to modify the IS scheme adaptively (e.g., by stopping IS earlier or later) in order
to reduce the variability of the quantity e -# Y 1 W (# ). We will return to this in
Section 6.
This discussion assumed implicitly that m # is not too large (e.g., for bounding
may expect the method to eventually break down as m # for
fixed B # , in view of the fact that the asymptotic overflow frequency for fixed B #
P. L'Ecuyer and Y. Champoux
generally di#ers from that when B # (e.g., [Shwartz and Weiss
1995]). Courcoubetis and Weber [1996] compare these asymptotics and give numerical
examples, for fluid source models, where the asymptotic based on large m #
gives a better approximation than the one based on large B # to the steady-state
probability that the queue exceeds B # . (The asymptotic based on large B # overestimated
the probability.) This suggests that the change of measure that we use might
perform poorly when m # is large and B # is small. On the other hand, the large
converges to the large B # asymptotic as B # for
their model. Empirically, we tried examples with m # up to 50 and the method still
worked fine. We also observed that taking # slightly less than # in our experiments
tended to improve the e#ciency when m # is large (see Section 6.1).
What about the variance of the variance estimator, with and without IS? They
can be compared by comparing -
Using the same crude
argument as in (11) above, one can conjecture the approximation
for some constant K 2 independent of B # . We expect not only the estimator itself
to be less noisy with IS than without, but also its sample variance to be less noisy,
and by a larger factor (at least if B # is large).
We now explain how the A-cycles are simulated to estimate both the numerator
and the denominator in (3), in the IS case. One simulates two versions of each A-
cycle, one with IS and the other without, both starting from the same initial state.
Thus, the A-cycles come in pairs. For the ith A-cycle pair, one first simulates an
A-cycle with IS, which provides an estimation L i Y i of the numerator, where L i and
Y i are the value of the L and the number of cell losses for this cycle. Then, the
state of the system is reset to what it was at the beginning of this A-cycle with IS,
and a second A-cycle is simulated to obtain an estimator X i of the denominator.
The final state of the no-IS A-cycle, which obeys approximately the distribution #,
is then saved and is taken as the initial state for the next pair of A-cycles. After a
warmup of n 0 cycles without IS, n pairs of A-cycles are thus simulated and the IS
estimator of - is
A confidence interval is computed using batch means as explained in Section 3.
Starting the two A-cycles of each pair from the same state means that we must
save or reset the entire state of the system after each cycle. This means copying
how many cells are at each node of the network, the destinations of these cells,
and the state (ON or OFF) of each source. We also memorize/reset the state of
each random number generator, so that the two A-cycles of a pair use common
random numbers. This tends to increase the correlation between L
to decrease the variance of - as a result.
Small Cell Loss Ratios in ATM via IS - 13
5. SIMULATION EXPERIMENTS
5.1 The Setup
For several examples and parameter sets, we ran the simulation first using the
standard approach without IS, for C A-cycles, and then with IS for C # pairs of
A-cycles. In each case, the values of C and C # were chosen so that the total CPU
time was about the same for both IS and no-IS, and these A-cycles were regrouped
batches. (For sensitivity analysis with respect to b, we tried di#erent
values of b ranging from 50 to 3200, for several examples, and found that the
variance estimates were practically independent of b, in that range, for the values
of C and C # that we use). For
Y j denote the samples
means of the X i and Y i (or of the X i and L i Y i , for IS), respectively, within batch j.
The tables that follow report the value of the CLR estimator - and of its variance
estimator
X, and -
Y
, and SXY are the sample means, sample
variances, and sample covariance of the -
respectively. The tables also
report the relative half-width -
2.57- of a 99% confidence interval on - (under
the normality assumption), the CPU time TCPU (in seconds) required to perform
the simulation, and the relative e#ciency (e#.), defined as -
These
values are all noisy estimates but give a good indication of what happens.
For the cases where no cell loss was observed in all the A-cycles simulated, we
put -
and the entries for the variance and e#ciency are left blank. The
simulation with IS takes more CPU time than no-IS for the same total number of
simulated cells, but the relative e#ciency takes both the variance reduction and
the overhead into account. Beware: E#ciencies and CPU times can be compared
within a given table, but not across the tables, because the models are di#erent
and the experiments were run on di#erent machines (SUN SparcStations 4, 5, and
20). Within each table, common random numbers were used for the corresponding
A-cycles across the di#erent lines of the table.
5.2 CLR Estimation at Level 2
Example 1. Let
vary the bu#er
. There are 50 sources feeding the target node q # (i.e.,
so the average arrival rate at q # is 50/101 # 0.495, while the service rate is 3.
With these numbers, we compute
increases the total arrival rate at q # from
0.495 to 14.48.
We took (note that the IS cycles
are much longer than the no-IS on the average, and their average length increases
most of them fill up the bu#er before emptying it again, whereas
for most of the no-IS cycles the bu#er empties after just a few cell arrivals). Table 1
gives the results. For not a single cell loss was observed, so
the estimates are useless. On the other hand, the relative error of the IS estimators
does not increase significantly as a function of B 2 , and these estimators can estimate
14 - P. L'Ecuyer and Y. Champoux
very small CLRs. The e#ciency decreases slowly with B 2 . (The outlier at
will be discussed later on.)
Example 2. Same as the preceding example, except that B 2 is now fixed at 512
and we vary the average burst size # 1 . For large # 1 , - is large and easy to estimate,
but not for small # 1 (the other parameters remaining the same). The results are in

Table

2. Without IS, cell losses were observed only for # 1 # 100, and even in that
case IS is more e#cient. The total arrival rate with IS decreases with
from 22.5 for 150. The squared relative error with IS
(not show in the table) is approximately constant as a function of # 1 .

Table

1. CLR estimation at level 2 for di#erent bu#er sizes
no-IS
128 2.8E-5 2.5E-11 45% 2828 0.0113
IS
128 3.0E-5 6.3E-13 7% 1675 0.838
512 2.5E-9 5.4E-20 24% 2593 0.043
768 3.7E-11 5.9E-22 170% 3108 0.001

Table

2. CLR estimation at level 2 for di#erent average burst sizes
no-IS
IS
50 2.5E-9 5.4E-20 24% 2593 0.043
100 7.2E-5 3.2E-12 6% 2445 0.659
An important question now arises: How noisy are the variance and e#ciency
estimates given in the tables? One way of estimating the distribution of the variance
and e#ciency estimators is to bootstrap from the b batch means, as follows. Put
the b pairs ( -
in a table. Draw b random pairs from that table,
with replacement, and compute the quantities - # 2 and e#. that correspond to this
sample of size b. Repeat this N times and compute the empirical distributions of
Small Cell Loss Ratios in ATM via IS - 15
the N values of - # 2 and of e#. thus obtained. These empirical distributions are
bootstrap estimators of the distributions of - # 2 and e#., and the interval between
the 2.5th and 97.5th percentiles of the empirical distribution is a 95% bootstrap
confidence interval for the variance of -
- or for the e#ciency. Table 3 gives the xth
percentiles Q x of the bootstrap distributions obtained from the results of Example 1,
2.5, 50, and 97.5, with

Table

3. Bootstrap quantile estimates for Example 1
128 4.3E-13 6.2E-13 8.8E-13 0.23 0.31 0.42
768 3.6E-25 5.9E-22 1.8E-21 2.4E-4 3.7E-4 3.7E-3
We already pointed out the very low empirical e#ciency of the IS estimator with

Table

1. A closer look at the 200 batch means -
that one
of the -
Y j in that case is all others are less than 10 -9 , except
one which is It seems that a rare event has happened within that
particular batch. We did not observe such outliers for the other values of B 2 , but
we found some in other examples, although rarely as excessive. The presence of
these outliers indicates that there remains an important tail in the distribution of
Y j with IS, despite the large reduction in the variance of -
This outlier has an
important e#ect not only on the variance and e#ciency estimators, but also on
the bootstrap distributions, as can be seen from Table 3 (compare the behavior
of the quantiles for those for the other values of B 2 ). To assess
the e#ect of the outlier, we repeated the bootstrap after removing it from the
sample (i.e., with the 199 remaining pairs), and obtained the following quantiles
4.5E-24). The e#ect is significant.
The numbers suggest that for the variance is highly overestimated, that
the e#ciency is underestimated, and that the bootstrap distribution is more widely
spread than the true distribution. To confirm these suspicions, we made 5 additional
replications of the entire experiment, independently, with IS. The
results, in Table 4, give an idea of the variability. Table 5 provides similar results
512. One can see that the e#ciency estimator is (unfortunately) noisy. On
the other hand, - is (fortunately) much less noisy, and this is reassuring. Another
reassuring empirical observation is that the distribution of the variance estimator is
skewed towards conservatism, in the sense that we observed large overestimation of
the variance (like here) from time to time, but did not observe large underestimation
of the variance. Moreover, the variance estimators tended to be more noisy when
large. Often, in practice, only crude estimates of - are sought (e.g., up
to a factor of 10 or so), so these variance estimators may su#ce. Otherwise, one
must increase the sample size. We warn the reader that in a situation where the
change of measure is not very e#ective, -
- might be noisy as well. This is in fact
what happens without IS. However, this never happened in our examples.
P. L'Ecuyer and Y. Champoux

Table

4. Five additional independent replications for
1.2E-11 1.1E-24 22% 3132 0.044
1.2E-11 2.1E-24 30% 3100 0.023
1.1E-11 8.0E-25 21% 3108 0.048
1.1E-11 1.8E-24 31% 3119 0.022
9.5E-12 2.9E-25 14% 3100 0.101

Table

5. Five additional independent replications for
3.1E-9 1.8E-19 35% 2592 0.020
2.5E-9 2.0E-20 15% 2587 0.115
2.4E-9 1.0E-20 11% 2588 0.223
5.3 CLR Estimation at Level 3
Example 3. Let
and we vary the bu#er size B
of the 60 sources to q # (so One node at level 2 is fed by 2 of these 6 hot
sources, while no other node at levels 1 and 2 is fed by more than 1 of them. Here,
and the total arrival rate at
q # is 6/21 without IS and 5.0 with IS. We take
The results appear in Table 6. IS works nicely while the no-IS observes no cell
loss except at the smallest bu#er size. With IS, the relative error and the relative
e#ciency are almost constant with respect to B # .
Example 4. Same as the preceding example, except that B 3 is fixed at 256
and we vary the average burst size # 1 . Table 7 gives the results. While no-IS has
di#culty to observe cell losses, IS gives reasonable estimations.

Table

6. CLR estimation at level 3 for di#erent bu#er sizes
no-IS
128 2.4E-5 1.1E-10 112% 7036 0.002
IS
128 4.1E-5 5.3E-12 14% 5779 0.056
768 2.5E-13 1.7E-28 13% 12930 0.029
Small Cell Loss Ratios in ATM via IS - 17

Table

7. CLR estimation at level 3 for di#erent average burst sizes
no-IS
100 2.1E-5 2.1E-10 178% 1134 0.007
IS
50 6.0E-7 7.2E-16 11% 1042 0.48
100 4.1E-5 6.7E-12 16% 881 0.28
5.4 CLR Estimation at Level 4
Example 5. Let
and we vary the bu#er size B
We assign 6 of the 300 sources to q # . They are distributed as in Example 3. Here,
and the total arrival rate at
q # is 6/41 without IS and 3.692 with IS. We take 000. The
results are in Table 8 and they resemble what was observed at level 3 (the e#ciencies
are smaller because we have a larger network and more sources to simulate). For this
example, we also varied # 1 with B 2 fixed at 512, and the results were qualitatively
similar to those of Table 7.

Table

8. CLR estimation at level 4 for di#erent bu#er sizes
no-IS
128 1.6E-3 7.5E-8 44% 3580 0.004
IS
128 1.1E-3 4.0E-9 15% 1881 0.15
5.5 Other Variants of the Model
We made several experiments with variants of the model to explore the e#ectiveness
of the proposed IS strategy in other (sometimes more realistic) situations.
The original model is called variant A. For variant B , the sources are no longer
a#ected to fixed destinations, but the destination of each cell is chosen randomly,
independently of other cells, uniformly over all destinations. Variant C is similar
P. L'Ecuyer and Y. Champoux
except that each burst (i.e., all the cells from a source during a given ON period) has
a random destination. The IS approach of Section 4 did very badly for variant B,
and gave improvement for variant C only when - was very small. An appropriate IS
strategy for these models should also change the probabilities over the destinations
to increase the tra#c towards q # . Variants B and C are not very realistic for ATM
switches, but the next 2 variants have higher practical interest.
In variant D, each node at level 3 has k bu#ers, the first one receiving the cells
originating from the sources 1 to the second one taking those from the
sources A server at level 3 takes cells
from those bu#ers according to either a round-robin or longest-queue-first policy.
In variant E, the sources produce two classes of cells: High priority constant bit
rate (CBR) cells and low priority variable bit rate (VBR) cells. The VBR sources
are Markov modulated as before, whereas the CBR sources have constant ON and
OFF periods (they are completely deterministic). Each node has two bu#ers, one
for the CBR cells and one for the VBR cells, and the CBR cells are always served
before the VBR ones.
The IS strategy of Section 4 works fine for the variants D and E: It provides
reasonable estimates for values of - that standard simulation cannot handle. We
also observed in our empirical results that the longest-queue-first policy gives a
CLR generally smaller than round robin.
6. REFINING THE IMPORTANCE SAMPLING SCHEME
6.1 Optimizing #
The IS approach of Section 4 provides a good change of measure, but based only
on a heuristic and asymptotic argument, not necessarily the optimal value of #
for a given bu#er size. Moreover, when choosing #, the approach does not take
into account the computational costs which may depend on #. To evaluate the
sensitivity with respect to #, we performed additional experiments where # was
varied around # , and the variance and e#ciency were estimated. As a general
rule, for the class of examples examined, we found that the optimal # was around
20% to 25% less than # , and increased the e#ciency by a factor between 2 to
compared with # , at level 2 or 3 where m # is typically large. At level 1 or 4, where
m # is usually small, the optimal # tends to be much closer to (and no significantly
better than) # . We emphasize that there is noise in these estimated factors, due to
the variance of the e#ciency estimators. However, the tendency persisted when we
replicated the experiments. Moreover, the # that gives the best e#ciency also tends
to reduce significantly the variance of the variance estimator. These improvements
are important, so it seems reasonable to use, e.g., of # at
levels 2 and 3, for networks similar to those that we have considered here, and
perhaps try to optimize # adaptively in a small neighborhood around that value,
during the simulation. Taking a smaller # is more conservative, in the sense that
it produces a smaller change of measure. Using # , which is more aggressive,
is very dangerous because the variance increases very fast with # in that area, and
may even become infinite for finite #. The next examples illustrate typical behavior
at levels 3 and 4.
Example 6. Let
Small Cell Loss Ratios in ATM via IS - 19
1/21. The node q # is fed by 6 sources, whose
tra#c passes through as in example 3. We take
and the results for di#erent values of # around # are in Table 9. Taking
improves the empirical e#ciency by a factor of approximately 25 compared with # .
By examining the data more closely, we found that the e#ciency improves because
the smaller # gives a smaller value of S 2
which is the dominant term in - # 2 .
Further replications showed similar results, with registering e#ciencies
15 to 60 times higher than # . The variance estimator was also much less noisy at
(where the value of - # was always between 3.6% and 4.2%) than at # .

Table

9. Comparing di#erent values of #, for
Example 7. Let
Only 2 sources feed the node q # .
Both sources feed the same node at level 3, but di#erent nodes at levels 1 and 2. We
take In this case, # = 0.0394, and we found empirically that taking
# brings no significant e#ciency improvement. We made similar experiments
with exactly the same data as in Example 5, with observed an
e#ciency improvement by a factor between 1.5 and 2.
6.2 Defining the A-Cycles Di#erently
Instead of starting the A-cycles when the bu#er at q # becomes empty, one can start
them when the number of cells in the bu#er crosses # upward, where # is a fixed
integer. There is essentially nothing to gain in that direction, however, because
when increasing # the no-IS A-cycles tend to become excessively long (typically,
the bu#er at q # remains nearly empty most of the time).
Another idea is to impose a lower bound, say t 0 , on the length of the A-cycles, to
get rid of the extremely short (and wasteful) A-cycles which tend to occur frequently
under both the IS and no-IS setup. The A-cycle ends at the maximum time between
t 0 and the first time when node q # becomes empty. How to choose t 0 ? We want to
choose it large enough to make sure that most A-cycles under IS see some overflow,
but not too large, so that the A-cycles end at the first return to the empty state
after overflow. According to our arguments in Section 4, if overflow occurs at
then the total production by the twisted sources up to time t 1 should be
approximately equal to the number of cells required to keep the server busy until
fill up the bu#er at node q # , that is, m # -
# is the
P. L'Ecuyer and Y. Champoux
average production rate of a twisted source. The additional time t 2 to empty the
bu#er (with IS turned o#) should satisfy . We want (roughly)
We suggest taking t 0 somewhere between 20% and 50% of the value of that upper
bound. In our experiments, this always gave e#ciency improvement. Since the
variance associated with the IS cycles is the dominant term in the variance of -
-,
a good strategy is to choose t 0 just large enough so that most of the IS cycles fill
up the bu#er. Taking t 0 too large (close to not a good idea because it
makes us spend too much time on the no-IS cycles without bringing much additional
variance reduction. Beyond a certain point, increasing t 0 eventually decreases the
e#ciency.
Example 8. We used the same data as in Example 5 (for
and with IS. For # , we have t 1 # 95 and t 2 # 300. For
we have a total arrival rate of 1.82 with IS, which give t 1 # 312 and

Table

give the results. With # , raising t 0 from 0 to 75 increases
the (empirical) e#ciency approximately by a factor of 4. With raising
t 0 from 0 to 150 improves the (empirical) e#ciency by a factor of more than 10.
This gain is related to the rapid increase of -
X, which decreases -
when t 0 is small. We made 2 additional replications of this experiment and the
results were similar, although the empirical e#ciency for # and t
and which suggests that the factor of e#ciency improvement from this
setup to more around 20 to 30 instead of 10. The variance
estimator is also (empirically) much less noisy with

Table

10. Imposing Lower Bounds on the A-cycle lengths, for
50 4.8E-5 8.8E-13 5.0 % 8.1 1.25 7089 0.37
100 4.8E-5 5.0E-13 3.8 % 15.8 2.43 10523 0.43
200 4.9E-5 7.6E-13 4.5 % 30.7 3.25 13890 0.23
50 4.8E-5 6.6E-13 4.4 % 8.1 1.24 4646 0.75
100 4.8E-5 3.0E-13 2.9 % 15.7 2.42 7886 0.98
200 4.8E-5 2.5E-13 2.6 % 30.7 3.24 12645 0.75
Example 9. Let
sources feed the target node
, as in Example 3, which gives an average arrival rate of 6/21 # 0.286 to that
Small Cell Loss Ratios in ATM via IS - 21
node. We run simulations for di#erent values of t 0 both with the - r ij associated to
arrival rate of 5.00, t 1 # 85, and t 2 # 150)
and
# 0.54, a total arrival rate of 3.26, t 1 # 200,
and t 2 # 150). Using t gave the best empirical
e#ciency in this case, about 20 times the empirical e#ciency observed with t
and # .
6.3 Stopping IS Earlier
Suppose that and that we use IS. When the target bu#er at q # overflows and
IS is turned o#, there may be several cells already in the network at previous levels,
and this may produce more overflow than necessary. Because of that, it could make
sense to turn o# IS earlier, e.g., when the total number of cells in bu#er q # or at
previous nodes but on their way to q # , reaches some threshold N 0 . Beck et al. [1999]
and Dabrowski et al. [1998] use this criterion for turning o# IS, with N Our
experiments with this idea showed no significant improvement compared with the
method which turns o# IS when q # overflows. With N seems to
reduce the e#ciency instead. Here is a typical illustration.
Example 10. Let
sources feed the node q # ,
which gives an arrival rate at q # of When IS is applied the arrival rate
increases to 1.5887. These 2 hot sources feed di#erent nodes at level 2. In Table 11,
Y is the average number of cell losses per cycle with IS and N
to turning o# IS when q # overflows. Taking N 0 between 520 and 600 appears to be
about as good as our usual method, but N 0 # 510 is definitely worse.

Table

11. Di#erent stopping criteria for IS
500 4.72E-9 6.7E-18 19.8% 1.51 0.188 3308 1.1
6.4 Retroactive Manipulations to Control the Overflow
The criterion for turning o# IS earlier, considered in the previous subsection, is
rather blind. Remember that all the randomness in our model is in the state
transitions of the sources. It is therefore possible, in principle, to compute at any
given point # in time whether or not there will be overflow at q # caused only by the
cells generated so far, and turn o# IS as soon as this happens (this is a stopping
time with respect to the filtration generated by the trajectory of S(-)). In this way,
IS is turned o# before the target bu#er fills up, but only when overflow is guaranteed
to occur. In practice, this can be implemented by actually running the simulation
until there is overflow, and then turning o# IS retroactively right after the time #
22 - P. L'Ecuyer and Y. Champoux
when all the cells having reached q # when the first cell overflows (at time #,
say) were already produced by a source. (This # remains a stopping time in our
model; retroaction would be used only to avoid too much computation at each time
step.) This is complicated to implement and implies significant overhead. Despite
spending a lot of time on experimenting with this idea, we were unsuccessful in
improving the e#ciency with it.
6.5 Combining IS with Indirect Estimation
Srikant and Whitt [1999] proposed the following indirect estimator of the CLR.
(This approach was presented by Ward Whitt during the keynote address of the
1997 Winter Simulation Conference.) The CLR at node q # satisfies
is the total (average) production rate of the m 0 sources
feeding node q # is the (average) output rate from node q # , 1/c # is the service
time at node q # , and # is the steady-state fraction of time where the server is busy
at node q # . The second equality follows from the conservation equation #/c # .
Using (13), - can be estimated indirectly by estimating #. Srikant and Whitt [1999]
showed that the indirect estimator brings substantial variance reduction in heavy
tra#c situations, especially for queues with several servers and random service
times, but not in light tra#c. In our context, the tra#c at q # is light, but becomes
heavy when IS is applied, so it was not clear to us a priori if the indirect estimator
combined with IS could help.
The results of our extensive numerical experiments can be summarized as follows.
For a single queue with several servers, without IS, the indirect estimator reduces
the variance by large factors when the total arrival rate exceeds the service capacity,
and increases the variance by large factors when the total arrival rate is much less
than the service capacity. This is true even for constant service times and single-server
queues, but less servers or less variability in the service times favors the
direct estimator. A larger bu#er at q # tends to accentuate the factor of variance
reduction or variance increase. When the indirect estimator was combined with IS,
we observed a variance increase instead of a variance reduction, even if the total
arrival rate after IS was larger than the service rate. An intuitive explanation seems
to be that because IS is turned o# as soon as the bu#er overflows, the condition
favoring the indirect estimator (sustained overloading at q # ) does not hold for a
large enough fraction of the time.
7. FUNCTIONAL ESTIMATION
So far we have considered the problem of estimating the CLR for fixed values of the
model parameters. But in real life one is often interested in a wide range of values
of the r ij 's and of the bu#er sizes. We now examine how the CLR can be estimated
in functional form, as a function of the matrix R, from a single simulation, and also
as a function of B # by re-using certain portions of the simulation.
Let R and -
R be as before, where -
R is the twisted version of R determined as in
Section 4, but suppose that we now want to estimate the CLR - for R replaced
by -
R, for several -
R in some neighborhood of R, by simulating pairs of A-cycles
with -
R and R only. This can be achieved as follows. One simulates pairs of A-
Small Cell Loss Ratios in ATM via IS - 23
cycles and computes X i , Y i , and the likelihood ratio L i for each pair just as before.
Afterwards, the estimators L i Y i and X i of the numerator and the denominator are
multiplied by the likelihood ratios
r 00
r 01
r 11
r 11
r 00
r 01
r 11
N #,
respectively, where N #
kl
and N #
kl
are the total number of transitions of the sources
from state k to state l during the A-cycle with IS and without IS, respectively. The
functional estimator of - is then
and it can be evaluated a posteriori for as many di#erent matrices -
R as desired, as
long as -
R is not too far away from R. The additional overhead during the simulation
amounts only to storing the values of N #
kl
and N #
kl
, together with those of X i and
all the pairs of A-cycles. This type of functional estimator based on a
likelihood ratio is discussed in a more general context by Rubinstein and Shapiro
[1993] and L'Ecuyer [1993], for example.
Example 11. We give an example of functional estimation at level 4. Let
sources to the node q # and take
R and run the simulation as usual, and then compute two functional estimators.
For the first one, # 1 is fixed and - is estimated as a function of # 0 , whereas for
the second one, # 1 fixed and - is estimated as a function of # 1 .

Tables

12 gives a partial view of the results for the first estimator.
The relative half-widths of pointwise 99% confidence interval, - #, remain reason-able
for a good range of values of # 0 and # 1 . If one is interested in a wider region,
that region can be partitioned into a few subintervals and a di#erent -
R can be used
for each subinterval.
For the estimation of - as a function of B # , one cannot use the likelihood ratio
approach, because B # is not a parameter of a probability distribution in the
model. However, observe that when an A-cycle is simulated, the sample path of
the system is independent of B # as long as there is no overflow at q # . Therefore,
when estimating - for several large values of B # , the initial part of the simulation
(until overflow occurs) does not have to be repeated for each value. One can start
with a single simulation (or sample path) and create a new subpath (or branch)
each time the number of cells at q # exceeds one of the bu#er sizes of interest. If
one is interested in N distinct values of B # , one eventually ends up with N parallel
simulations, but a lot of work is saved by starting these parallel simulations only
when needed. This type of approach is studied in more generality in L'Ecuyer and
V-azquez-Abad [1997]. In our experiments with this method, the savings in CPU
time were typically around 50%.
P. L'Ecuyer and Y. Champoux

Table

12. Functional estimation at level 4, for fixed # 1
500 2.4E-6 6.5E-14 27%
588 9.8E-7 4.0E-15 17%
714 3.5E-7 1.9E-16 10%
800 2.0E-7 4.3E-17 7%
909 1.1E-7 1.0E-17 7%
The development of Section 4 suggests an approximately linear relationship between
, at least asymptotically. Our empirical experiments confirm that
the linear model fits very well for large enough B # . We can therefore
recommend, for estimating - as a function of B # , to perform simulations at 4
or 5 large values of B # only, and fit a linear model to the observations (B # ,
-) by
least squares regression.
8. CONCLUSION
We have discussed how to implement IS for estimating the CLR in an ATM switch
model, by a direct adaptation of the approach proposed in Chang et al. [1994]. Our
extensive empirical experiments with large networks indicate that the method is
viable, at least for the type of model and examples that we have considered. IS
improves the CLR estimator and also improves, by a larger factor, the corresponding
variance estimator (although this variance estimator often remains noisy even with
IS, in which case the natural way of improving the reliability of the confidence
intervals would be to increase the sample size). We have also explored heuristic
refinements of the IS scheme. Some gave improvements with respect to the basic
scheme, others did not. Those that improved the e#ciency of the CLR estimator,
namely taking # slightly less than # and imposing a lower bound on the A-cycle
length, also (empirically) reduced the noise of the variance estimator.
This IS methodology can be adapted to other variants of the model, in addition
to those considered in Section 5.5, e.g., by replacing the geometric sojourn time
distributions in the ON and OFF states by other distributions with an exponential
tail. Statistical studies of tra#c traces suggest that heavy-tailed distributions,
with infinite variance, are more appropriate, although some people still argue that
the exponential-tailed distributions do the job well-enough in many cases. The IS
scheme considered in this paper does not appear easily adaptable to heavy-tailed
distributions. Finding an e#ective way of applying IS for such distributions remains
a challenging problem.



--R

A unified approach to fast teller queues and ATM.

Accelerated simulation of a leaky bucket controller.
A Guide to Simulation (Second
Estimation du taux de perte de r-eseaux ATM via la simulation et le changement de mesure
Sample path large deviations and intree networks.

Fast simulation of packet loss rates in a shared bu

Accelerated simulation of ATM switching fabrics.
Decoupling bandwidths for networks: A decomposition approach to resource management.
Fast simulation of networks of queues with e
Monte Carlo: Concepts
Quality Service in ATM Networks.
Multilevel splitting for estimating rare event probabilities.

Importance sampling for stochastic simulations.
Management Science

Fast simulation of rare events in queueing and reliability models.
ACM Transactions on Modeling and Computer Simulation
What are the implications of long-range dependence for VBR video tra#c engineering
Modeling video tra
Simulation Modeling and Analysis (Third
Two approaches for estimating the gradient in a functional form.

Importance sampling for large ATM-type queueing networks


On the self-similar nature of ethernet tra#c (extended version)
Application of the M/Pareto process to modeling broadband tra
Fast simulation of steady-state availability in non-markovian highly dependable systems
A quick simulation method for excessive backlogs in networks of queues.
Discrete Event Systems: Sensitivity Analysis and Stochastic Optimization by the Score Function Method.
Large deviations and e
On the optimality and stability of exponential twisting in Monte Carlo estimation.
On Monte Carlo estimation of large deviation probabilities.
Importance sampling for the simulation of highly reliable markovian systems.
Large Deviations for Performance Analysis.
Variance reduction in simulations of loss models.

--TR
A guide to simulation (2nd ed.)
Importance sampling for stochastic simulations
On the self-similar nature of Ethernet traffic (extended version)
Importance sampling for the simulation of highly reliable Markovian systems
Efficiency improvement and variance reduction
Effective bandwidth and fast simulation of ATM intree networks
Fast simulation of rare events in queueing and reliability models
through high-variability
Fast simulation of packet loss rates in a shared buffer communications switch
What are the implications of long-range dependence for VBR-video traffic engineering?
Importance sampling for large ATM-type queueing networks
Two approaches for estimating the gradient in functional form
Fast simulation of networks of queues with effective and decoupling bandwidths
Simulation Modeling and Analysis
Quality of Service in ATM Networks
Functional Estimation with Respect to a Threshold Parameter via Dynamic Split-and-Merge
Multilevel Splitting for Estimating Rare Event Probabilities
Variance Reduction in Simulations of Loss Models
Application of the M/Pareto Process to Modeling Broadband Traffic Streams

--CTR
P. T. de Boer , D. P. Kroese , R. Y. Rubinstein, Rare event simulation and combinatorial optimization using cross entropy: estimating buffer overflows in three stages using cross-entropy, Proceedings of the 34th conference on Winter simulation: exploring new frontiers, December 08-11, 2002, San Diego, California
Rick Siow Mong Goh , Ian Li-Jin Thng, Twol-amalgamated priority queues, Journal of Experimental Algorithmics (JEA), v.9 n.es, 2004
Victor F. Nicola , Tatiana S. Zaburnenko, Efficient heuristics for the simulation of population overflow in series and parallel queues, Proceedings of the 1st international conference on Performance evaluation methodolgies and tools, October 11-13, 2006, Pisa, Italy
Sandeep Juneja , Victor Nicola, Efficient simulation of buffer overflow probabilities in jackson networks with feedback, ACM Transactions on Modeling and Computer Simulation (TOMACS), v.15 n.4, p.281-315, October 2005
Victor F. Nicola , Tatiana S. Zaburnenko, Efficient importance sampling heuristics for the simulation of population overflow in Jackson networks, ACM Transactions on Modeling and Computer Simulation (TOMACS), v.17 n.2, p.10-es, April 2007
variance reduction;importance sampling;ATM;rare events