Fast algorithms for generating discrete random variates with changing distributions.
One of the most fundamental operations when simulating a stochastic discrete-event dynamic system is the generation of a nonuniform discrete random variate. The simplest form of this operation can be stated as follows: Generate a random variable X that is distributed over the integers 1,2,,n such that are fixed nonnegative numbers. The well-known alias algorithm is available to accomplish this task in O(1) time. A more difficult problem is to generate variates for X when the ai's are changing with time. We  present three rejection-based algorithms for this task, and for each algorithm we characterize the performance in terms of acceptance probability and the expected effort to generate a variate. We show that, under fairly unrestrictive conditions, the long-run expected effort is O(1). Applications to Markovian queuing networks are discussed. We also compare the three algorithms with competing schemes appearing in the literature.
Introduction
The problem of generating a nonuniform discrete random variate is fundamental to the
simulation of any discrete event stochastic system. The simplest version of this problem
is to generate a variate for a random variable X such that P
given that the p i 's do not change with time. The well-known `alias algorithm' (e.g., see [2]
pp. 158-160) takes O(n) preprocessing time, after which it can generate a variate in O(1)
time. In this paper we are interested in developing efficient algorithms for the case of the
changing with time. Such a procedure is desirable in the simulation of multi-dimensional
Markov processes, such as those used to model queueing and telephone networks.
A definition of the problem is as follows. Suppose that for each are
given n non-negative numbers a 1 be a sequence of
independent random variables such that
a
a
The problem is to generate a variate for X(s) for each realization
of fX(s)g). We refer to the integers in ng as the 'outcomes' of X. Call a i (s) the rate
for outcome i at time s. In the static version of the problem, the rates do not change with
time. In the dynamic version, zero or more of the a i (s)'s change each time s is incremented.
In discrete event simulation the number of a i (s)'s that change each time s is incremented
is typically small. The algorithms developed in this paper attempt to exploit this special
property.
In Section 2 we first show how the rejection algorithm can be employed to generate
variates with changing distributions. The performance of the algorithm can be characterized
by its acceptance probability. We then develop two additional algorithms which, at the
expense of additional memory, can improve significantly the acceptance probability. In
Section 3 we suppose that the rates are random and are determined by a (discrete state
space) Markov process, as is typically the case in applications. We explicitly characterize
the long-run acceptance probability and the long-run expected effort of the three algorithms
in terms of the steady state rates of the underlying Markov process. In Section 4, several
classes of queueing networks are considered, and for each class we characterize the long-run
performance of the algorithms in terms of the defining parameters of the network. In Section
5 we compare the algorithms of this paper with the event-list approach, assuming that the
event list is organized as a heap, and with the algorithm TRANSIT proposed by Fox [8]. In
Section 6 we discuss several methods for improving the performance of the three algorithms
discussed in this paper.
Three Algorithms
The following assumptions will be in force throughout the paper.
(A1) For each exists a finite -a i such that a i
.
There exists a finite b ? 0 such that
a
We refer to (-a 1 as the majorizing n-vector. All the three algorithms presented
in this paper are rejection based. Any rejection algorithm that operates with respect to a
majorizing n-vector runs in O(1) average time, provided the acceptance probabilities are
bounded away from zero as n increases. Traditionally and correctly, however, a major theme
in the literature on rejection algorithms is making the implicit constant (depending on the
reciprocal of the acceptance probability) small. We quantify this constant in terms of the
a i 's. In particular we quantify this constant for product form instances of the queueing-
network models in terms of the traffic conditions. Clearly, though, the heavier the traffic at
all nodes, the smaller will be the constant for essentially any network.
We call our three algorithms Generate1, Generate2, and Generate3. Generate1 is an
implementation of a discrete version of standard rejection using the dynamic rate vector
with respect to the uniform majorizing n-vector (-a;
erate3 is an implementation of a discrete version of standard rejection using the dynamic
rate vector with respect to the non-uniform majorizing n-vector (-a
erate1 and Generate3 are based on assigning a single 'bucket' for each event. Generate2
combines the principle in Generate1 with assignment of (possibly) more than one bucket to
any outcome. Bucket schemes have been considered before in the literature, but apparently
not linked directly to rejection algorithms.
2.1 Algorithm Generate1
We now present our first algorithm for generating variates with a changing distribution. Our
algorithm makes use of n buckets, i.e., an array B of size n. Each bucket has a 'capacity' of 1.
For each bucket i we initially set B[i] = a i (0)=-a. This completes preprocessing. Clearly this
processing can be completed in O(n) time. Notice that the initial contents of each bucket
does not exceed unity. Here is an algorithm to generate a variate for X(0): (Hereafter let U
stand for the uniform random variable in the interval (0; 1).)
Algorithm Generate1
Step 1. Generate U , and compute I = dnUe.
Step 2. Let
Step 3. If R - B[I], accept and output I and quit;
Else go to Step 1.
Observe that Generate1 picks a bucket at random (uniformly). If I is the bucket chosen,
then I is output with probability B[I]. (Notice that R is a uniform variable in the interval
(0; 1) and is independent of I.)
The algorithm Generate1 is the rejection algorithm (e.g., see [2] pp. 151-152) for generating
discrete random variates with fixed distributions. (The rejection algorithm is typically
discussed in the context of continuous random variables.) It follows that the output of the
algorithm Generate1 has the same distribution as that of X(0).
The key observation we make here is that Generate1 can also be used to generate variates
with changing distributions. Indeed, suppose Generate1 is used to generate a variate for
X(s). (Thus (s)=-a.) If at time s
we simply reset 1)=-a. The procedure is valid since, by Assumption (A1), the
contents of each bucket never exceeds unity.
In order to characterize the performance of Generate1, let "acceptance at time s" be the
event that R - B[I] the first time Step 3 is invoked when generating a variate for X(s).
We refer to P ("acceptance at time s") as the acceptance probability at time s. Note that the
expected number of uniform random variates U needed to generate a variate for X(s) by the
algorithm Generate1 is equal to the reciprocal of the acceptance probability.
Theorem 2.1 The acceptance probability at time s for Generate1 is given by
P ("acceptance at time
n-a
Proof: We have
P ("acceptance at time (R - B[I])
(R - B[i])
a
-a
establishing the result. 2
From Theorem 2.1 we know that the performance of Generate1 will be poor at time s when
the average of the a i (s)'s is significantly less than the uniform upper bound -
a. In particular,
Generate1 does not exploit the fact that we may have - a i !! -a for many outcomes i, as is
often the case in applications. The algorithms discussed below, Generate2 and Generate3,
take advantage of this special structure.
2.2 Algorithm Generate2
In Generate1 one bucket is assigned to each outcome i, n. The idea behind
Generate2 is to assign one or more buckets to any outcome i. The number of buckets
assigned to i is proportional to the upper bound of the ith rate. Specifically, let d be a fixed
positive number and let the number of buckets assigned to outcome i, denoted by l i , be given
by
- a i
d
Also let be the total number of buckets. Note that we have the following
bounds on the total number of buckets:
d
d
Thus if
In Generate2, the number of buckets assigned to any i is held fixed even when the a i (s)'s
change. Once having chosen at random one of the l buckets, if that bucket is assigned to i,
then we will output i with a certain probability. This probability, however, will change as
the a i (s)'s change.
To be more precise, we make use of two arrays C and B of size l and n, respectively. For
B we set
We set C[j] equal to the outcome i assigned to bucket j (for 1
Preprocessing
Step 1. For each outcome i compute l i .
Step 2. Compute the prefix sums of (l the sums be
Step 3. Initialize C as follows: Fill cells 1 through m 1 with 1, cells
Step 4. Set
Next we present the procedure for generating a variate for X(0).
Algorithm Generate2
Step 1. Generate U and compute
Step 2. Let
Step 3. Let I = C[J ]. If R - B[I], output I and quit;
Else go to Step 1.
As with Generate1, this algorithm can also be used to generate a variate for X(s). After
having generated a variate for X(s), if at time s
we simply reset . Note that, since a i never exceeds unity.
Also note that we can always set - a
a, in which case Generate2
reduces to Generate1. In a manner entirely analogous to that for Generate1, we can define
the event "acceptance at time s" for the algorithm Generate2.
Theorem 2.2 At time s, the output of Generate2 has the same distribution as that of X(s).
Moreover, for Generate2 we have
P ("acceptance at time
a
d
Proof: We have
P ("acceptance at time
l i
l
from which (2) follows. Also note that
"acceptance at time
l
B[i]:
Combining the above two expressions gives P ij"acceptance at time
which completes the proof. 2
Theorem 2.2 directly gives the following result.
Corollary 2.1 We have the following bounds on the acceptance probability for Generate2:
P ("acceptance at time s") -
Furthermore, if d is a divisor of - a i , ("acceptance at time
a
Since
can be significantly less than n-a, the algorithm Generate2 can give a much
higher acceptance probability as compared with Generate1.
2.3 Algorithm Generate3
Our last algorithm is similar in spirit to Generate2. It has the advantage of always achieving
the upper bound in Corollary 2.1 for the acceptance probability. It has the disadvantage
that the effort required to pass through one iteration of the algorithm is somewhat greater
than that required by Generate2. Generate3 makes use of the alias algorithm and its data
structures. It also makes use of an array B of size n. Initially we set . The
preprocessing, including that required by the alias algorithm, is O(n). An equivalent version
of this algorithm has been independently discovered by Fishman [5].
Algorithm Generate3
Step 1. Use the alias algorithm to generate a variate I with distribution
a
Step 2. Generate a uniform variate U .
Step 3. If U - B[I], accept and output I and quit;
Else go to Step 1.
We leave it to the reader to establish the following analog to Theorem 2.2.
Theorem 2.3 At time s, the output of Generate3 has the same distribution as that of X(s).
Moreover, for Generate3 we have
P ("acceptance at time
We now give a brief comparison of the above three algorithms. Throughout this comparison
we suppose that
Generate2; hence, from Corollary 2.1, we know
that the acceptance probability for Generate2 is greater than 1=2 of the acceptance probability
for Generate3. The chart below gives a comparison of the memory requirements and
the operation counts for the three algorithms. We have assumed the worst case memory
count for Generate2. For the operation counts we have assumed that each of the following
operations count as 1: multiplication, subtraction, upper floor operation, table look up,
comparison, uniform variate generation. We define the operation count as the maximum
number of operations performed in one pass through the steps of the algorithm.
Memory Operation
Count
2n integers
Generate3
integers
In particular, if the acceptance probability of Generate2 is close to that of Generate3, we
may conclude from the above chart that Generate2 is preferable to Generate3.
We also mention that if d is a divisor of each rate and if the rates do not change with
time, Generate2 can be streamlined because streamlined
version of Generate2 will run faster than the alias algorithm. The streamlined algorithm is
similar, but not identical, to the algorithm developed by Fox (see [2], pp. 157-158). (Fox's
algorithm has more flexibility due to its choice of the parameter m.)
2.4 A Performance Measure for Randomized Algorithms
?From Theorem 2.3 we know that the expected number of U 's needed to generate a variate
for X(s) is not greater than
=b. However, Theorem 2.3 does not exclude the possibility
that the number of U 's needed to generate X(s) exceeds
probability. Analogous statements can be made for Theorem 2.1 and Theorem 2.2.
In order to study this question further let
are independent and
identically distributed geometric random variables with parameter p. (V can be thought of
as the number of times a coin has to be flipped before a head appears for the mth time, p
being the probability that a head appears in a single flip). Chernoff bounds (see e.g., [11],
pp. 388-393) can be used to obtain tight bounds on probabilities in the tail ends of V . In
particular we can show the following:
for any fixed 0 ! ffl !! 1 (see [15]). We will make use of these bounds in our analysis.
Just like the O() function is used to specify the asymptotic resource bounds of deterministic
algorithms, e
O() is used to specify resource (like time, space etc.) bounds of randomized
algorithms. e
O() is defined as follows.
Definition. [15] We say a function f(:) is e
O(g(:)) if there exist constants c and n 0 such that
any input of size n - n 0 , for any ff ? 0.
We say a function f(:) is !(g(:)) if lim n!1
Theorem 2.4 Let h be a lower bound on the acceptance probability of any of the above three
algorithms. If this algorithm is called times, then the total number of U's generated by
this algorithm is no more than (1
with probability fixed
Proof. Follows from the Chernoff bound and the observation that the number of U 's generated
by any of the three algorithms to obtain a variate for X is upper bounded by a geometric
random variable with parameter h. 2
Corollary 2.2 If n), the number of U's generated is no more than (1
with
probability -
for any constant 0 ! ffl !! 1 and any fixed ff - 1. Thus the number
of U's generated is e
O( m
3 Random Rates
In discrete event simulation the rates typically evolve in a random manner. In this context
we shall write A i (s) for the ith rate at time s in order to emphasize the fact that the
rates are random. Also write
for each s. In the context of random rates, we no longer require that the variables X(s),
independent. We do require, however, that be independent of
given A(s). It is easily seen that Generate1, Generate2, and Generate3
satisfy this requirement.
Assumptions (A1) and (A2) translate directly to the current context with a i replaced
by A i (s). The acceptance probabilities in Theorems 2.1, 2.2, and 2.3 have their analogs in
the context of random rates: replace a i (s) by E[A i (s)].
In most applications (see Section 4), there is an underlying (discrete state space) Markov
process which defines the rates A 1 (s). In this case Y (s) is the total jump rate
of the Markov process just before the sth jump. It is convenient to make the following
assumption.
(A3) The underlying Markov process is irreducible, aperiodic, and positive recurrent

With Assumption (A3) in force, Y (s) will converge in distribution to a random variable Y
that is the total rate of the process just before a jump in steady state. Let Z be the total
rate of the process at an arbitrary time instant in steady state (which is also well defined due
to Assumption A3). The random variables Y and Z take on the same discrete values but
typically have different distributions. In fact, it follows from Theorem 2.10.6 (a) of Walrand
[16] that
y
The following theorems enable us to characterize the long-run performance of the algorithms
in terms of the steady-state behavior of the underlying Markov process. (Of course,
not all simulations are steady state.) For each of the three algorithms, let V (s) be the number
of random variables U that are called in order to generate a variate for X(s). We shall
refer to V (s) as the effort required to generate a variate for X(s).
Theorem 3.1 Under assumptions (A1)-(A3), Generate1 has the following long-run acceptance
probability and long-run expected effort:
lim
P ("acceptance at time
lim
Proof: From Theorem 2.1 we have
P ("acceptance at time
n-a
and
Taking expectations and applying Assumption (A3) to both of these equalities gives
lim
P ("acceptance at time
n-a
and
lim
n-a
Y
But from (4) we have E[Y combined
with (5) and (6), give the desired result. 2
The proofs of the following two theorems are entirely analogous to the proof of Theorem
3.1.
Theorem 3.2 Under assumptions (A1)-(A3), Generate2 has the following long-run acceptance
probability and long-run expected effort:
lim
P ("acceptance at time
E[Z]d
lim
d
Theorem 3.3 Under assumptions (A1)-(A3), Generate3 has the following long-run acceptance
probability and long-run expected effort:
lim
P ("acceptance at time
lim
4 Application to Queueing Networks
In this section we give examples for which it is possible to explicitly calculate the moments
appearing in the theorems of section 3. Although some of these examples are better studied
by analytical methods than by simulation, the results give insight into the run times of the
generate algorithms.
Queueing networks are used extensively in the modeling and analysis of computer systems
and networks (e.g., see [12]). To illustrate our ideas for variate generation, consider the open
Jackson network (e.g., see [16]), one of the most fundamental queueing networks. There are n
single-server queues. Customers arrive from the outside according to a Poisson process with
rate - and are routed to queue i with probability r 0i , Service times at queue i
are exponentially distributed with parameter - i , When a customer completes
service at queue i, it is routed to queue j with probability r ij ; it leaves the network with
probability r i0
Service times are assumed to be independent of each other
and of the arrival process.
We will assume that there exists a nonnegative solution (- to the 'traffic equa-
such that ae i := - i =- This assumption assures that the underlying
Markov process is irreducible, aperiodic, and positive recurrent.
Now consider simulating the Markov process associated with this queueing network. Note
that the rate associated with queue i, A i (s), is either - i or 0 depending on whether queue i
is occupied or empty. Note that the rate associated with external arrivals, A 0 (s), is -. Given
that the current rates are A 0 (s), the elapsed time (the inter-event time) until the
next arrival or service-completion event is exponentially distributed with parameter Y
This event is a service completion at queue i with probability A i (s)=Y (s),
an arrival with probability A 0 (s)=Y (s). After having determined the
elapsed time and event type, the alias algorithm can be used to determine, in O(1) time, the
queue to which the customer is to be routed (including the possibility of being routed to the
outside). A new iteration then begins (i.e., s is incremented), during which we determine the
next inter-event time, the next event-type, and the the next queue to which the customer is
routed. Note that at most two of the rates change from iteration to iteration.
Consider the central component of the above simulation procedure: determining the event
type according to the probabilities A i (s)=Y (s), any one of the three
algorithms can be applied. We will now characterize the performance of Generate3. (It will
then be a simple exercise for the reader to do the same for Generate1 and Generate2.) Note
that in this application we have - a
Theorem 4.1 The long-run acceptance probability for Generate3 for the Jackson network
is given by
lim
P ("acceptance at time
The long-run expected effort for Generate3 for the Jackson network is given by
lim
Proof: Let L i denote the random variable for the number of customers present at queue i
in steady state. It is well known that
Y
Since
it follows from (7) that
and
Combining these last two expressions for E[Z] and var[Z] with Theorem 3.3 completes the
proof. 2
We see from Theorem 4.1 that if the time average utilization, ae i , at each queue i is at least
1=2 (as is the case in most applications of practical interest), then the long-run acceptance
probability is greater than 1=2. Theorem 4.1 also tells us that it is particularly desirable to
have high values of ae i for those queues i with high service rates.
Let us now attempt to characterize the performance for some other queueing networks.
First consider the Jackson network described above, but with m i servers at queue i. Now
suppose that ae these conditions, the associated Markov process is
again irreducible, positive recurrent, and aperiodic. In this case A 0
on values 0; - Following the proof of Theorem 4.1 it can be
shown that the long-run expected effort of Generate3 is given by
lim
We leave it to the reader to derive the corresponding expression for the long-run acceptance
probability. It is easily shown from (8) that if
then
lim
for Generate3. Now consider a sequence of queueing networks, where the nth network in the
sequence has n stations. If each network in the sequence satisfies (9) for some ffl ? 0, then
the long-run average expected effort is O(1).
Second consider the case of multiple classes, where each class c, exogenous
arrival rate - c and has routing probabilities r c
Suppose that there
is one server at each queue, and the service rate is given by - i which does not depend on
the class. Also suppose that the service discipline at each queue is First-Come-First-Serve
(FCFS). There are now C sets of traffic equations; suppose each equation has a nonnegative
solution (- c
suppose that
modified definitions and assumptions, it is
well known that (7) continues to hold true for this multiclass network. Hence, Theorem 4.1
holds unchanged.
Third consider a multiclass network with a single server at each node, as discussed above,
but now suppose that the network is closed, i.e., customers neither enter nor leave the network
so that the population size is fixed for each class. The finite state Markov process associated
with this network is irreducible and aperiodic. We now have
lim
where L i is again the number of customers present at node i in steady state. The quantity
can be expressed in terms on the defining parameters through 'normalization
constants'.
The networks discussed above are 'product-form' queueing networks. Let us now consider
a non-product-form network. In particular, consider the closed, multiclass, single-server
network discussed in the paragraph above, but now suppose that the service rate for class
c customers at queue i is - c
i (i.e., the service rates now depend on the class as well as the
queue). Thus A i (s) takes on values 0; - 1
n. The maximum service rate
at queue i is now given by -
Cg. Let
class-c customer is being served at queue i in steady state"):
We now have
lim
Unfortunately, the current tools of queueing theory do not offer a means for expressing fl c
in
terms of the defining parameters for this non-product-form network. But there are numerous
techniques available to approximate these quantities (e.g., see [12]).
Competing Methods
We now compare our algorithm with several competing schemes. We carry out the comparison
in two contexts: multiserver queueing networks and systems with similarity.
5.1 Multiserver Queueing Networks
Consider a network of n queues in series, where each queue is equipped with m servers,
each server operates at rate -, arrivals from the outside arrive at rate -, and each customer
traverses the n queues in sequence. We assume that m is large (say greater than 100). We
shall assume that the traffic is moderate to heavy (i.e., -=m- is between .5 and .9). We must
stress that the complexity bounds given in this Section are derived for this specific example;
they do not necessarily extrapolate to other examples.
The first competing method that we discuss makes use of a list of future events. (We
assume that the reader is familiar with this standard approach of simulating a discrete
event system.) For this method there is a data structure which contains the time of the
next external arrival and the times of the service completions of all the customers currently
in the network. Thus this data structure contains at least 1 and no more than nm
events. Assuming that the data structure is organized as a heap, this method requires
O(nm) memory and O(log n log m) worst-case effort to generate an event. Under the
assumed traffic conditions, the effort will be close to the worst case.
The second competing method also employs a list of future events. But instead of being
implemented across the servers, the future event list is implemented across queues. In this
case the data structure contains the time of the next arrival and, for each nonempty queue,
the time of the next service completion. Again assuming that a heap is employed, this
method requires O(n) memory and O(log n) worst-case effort. Note that the complexity
bounds for heaps hold for all traffic conditions, not just for the case of moderate to heavy
traffic.
The third competing scheme is the algorithm TRANSIT proposed by Fox [8]. Employing
Fox's notation, the fi s vector takes the form fi
is the number of busy servers at queue i at time s. The memory
required by this approach is O(n). It can be shown that the long-run expected effort of
TRANSIT is O(n). (Indeed, employing the notation of Fox, for any choice of fi 0 , jS i j will
typically be close to n and w i will typically be close to zero; it follows then that the average
complexity in this typical situations is O(n).) However, TRANSIT can be modified to have
a run time of O(log n) as is shown in [9] - though with a high implicit constant.
The papers [8] [9] discuss TRANSIT in the context of being implemented across the
queues, as discussed in the paragraph above. However, analogous to the first competing
method, TRANSIT can also be implemented across servers [10]. The fourth competing
scheme is therefore TRANSIT implemented over servers, which we will refer to as TRAN-
SIT(tailored). In this case the fi s vector takes the form fi
depending on whether the jth server at time s is busy
or not. (We assume that the nm servers are numbered.) We also suppose that fi 0 is set to
the maximum rates, i.e., fi 0
TRANSIT simplifies:
Algorithm TRANSIT(tailored)
Step 1. Use the alias algorithm to generate a variate J with distribution
Step 2. If fi s (J) ? 0, accept and output J and quit;
Else go to Step 1.
The above algorithm TRANSIT(tailored) was independently discovered by Fox [10] [7]. Note
that TRANSIT(tailored) is just TRANSIT with respect to a majorizing vector. Clearly,
TRANSIT(tailored) is an implementation of standard rejection with respect to a non-uniform
n-vector. In remarks T10 and T11 of [9], this case of TRANSIT was mentioned; however its
implementation across servers (except for networks of single-server queues) was not noted
there but, independently of this paper, was noted in [7]. Note that in the case
TRANSIT(tailored) and Generate3 are identical. In an analogous fashion, we can define
"tailored" algorithms for Generate2 and Generate3. It is easily seen that Generate3(tailored)
is equivalent to TRANSIT(tailored), whereas Generate2(tailored) is different.
Mimicking the analysis in Section 4, we see that TRANSIT(tailored) has a long-run
expected effort given by
lim
-n
Here V (s) should be interpreted as the number of times Step 1 is invoked in TRAN-
SIT(tailored) when generating an event at time s. It follows from (10) that the long-run
expected effort of TRANSIT(tailored) is O(1) for this sequence of models. However, the
memory requirements are O(mn).
We therefore arrive at the following table:
Memory Long-Run
Expected Effort
Future O(n) O(log n)
TRANSIT O(n) O(log n)
TRANSIT O(mn) O(1)
Generate2, O(n) O(1)
Generate3
It is therefore clear that, for these traffic conditions and for a large n, the future event list
and TRANSIT are not competitive.
We must stress that in the complexities of the above chart there is an implicit constant
factor depending on the reciprocal of the acceptance probability. The constants for the algorithms
in this chart are comparable except for the effort of TRANSIT, whose constant is
substantially larger. Let us compare these constants for TRANSIT(tailored) and Generate3.
We first observe that the acceptance probabilities for the two algorithms are identical. How-
ever, TRANSIT(tailored) requires one less operation than Generate3 for each pass through
the steps of the algorithm. Therefore the run times of the two algorithms are going to be
very close, with TRANSIT(tailored) being slightly faster. The acceptance probability of
Generate3 and of TRANSIT(tailored) is greater than that of Generate2, but never more
than twice as much if
Based on these observations and the discussion at the end of Section 2, we can make the
following conclusions. First suppose that memory is not an issue. Then TRANSIT(tailored)
gives slightly better performance than Generate3. As compared with Generate2, TRAN-
SIT(tailored) requires a few more operations per pass through the steps but has a higher
acceptance probability. However, the acceptance probability of Generate2 can be made arbitrarily
close to the acceptance probability of TRANSIT(tailored) by decreasing d. Hence,
for sufficiently small d, Generate2 should be slightly faster than TRANSIT(tailored).
Now let us suppose that memory is an issue. For Generate2 we suppose that
so that the memory requirements of Generate2 and Generate3 are roughly
the same. Because TRANSIT(tailored) requires approximately m times the memory required
by Generate2 and Generate3, and m - 100, we eliminate it from the competition. As
mentioned at the end of Section 2, there is a trade-off between operations per pass and
acceptance probability when comparing Generate2 and Generate3 (although the acceptance
probability of Generate2 is greater than one half of that of Generate3). We recommend that
the user perform pilot runs when choosing between the two algorithms.
In this example we have assumed that the servers and customers are homogeneous. If
this is not the case, e.g., if the various servers were in various phases of a PH-distribution,
keeping track of the individual servers is needed. For the inhomogeneous case, either a heap
implemented across servers (i.e., the first competing scheme) or one of the tailored algorithms
would have to be implemented.
5.2 Systems with Similarity
For the definition of similarity please refer to [8, 9]. Under similarity, the general version of
TRANSIT has O(1) complexity - even when implemented across queues in queueing network
settings and/or when similarity is with respect to an n-vector that does not majorize a
dynamic rate vector. However, under similarity TRANSIT(tailored), Generate1, Generate2,
and Generate3 perform badly for queueing networks where some nodes have light traffic and
for the (common) reliability model mentioned in [8]; for the later, arguably, the acceptance
probabilities for Generate1, Generate2, and Generate3 are not bounded away from zero as n
increases unless the growth in repair rates and the decrease in failure rates are unnaturally
restricted. The point of the general version of TRANSIT is that similarity does not imply
that the reference vector globally majorizes.
6 Improving Performance
Although the algorithms Generate2 and Generate3 should be useful for many applications,
there are situations for which their performance will be less than satisfactory. In particular,
unsatisfactory performance will occur when the expected rates are typically far from their
respective upper bounds, i.e., when
. In this section we discuss
three modifications of Generate2 which are designed to alleviate or even overcome this prob-
lem. In the paper ([14]) we discuss how some of these modifications can be employed in
the efficient simulation of large-scale telephone networks. Throughout this Section we shall
simplify the notation by suppressing all references to the time variable s.
Before discussing these modifications, it will be beneficial to introduce yet another algorithm
for generating variates with changing distributions. The same algorithm has also
been independently discovered by Fishman and Yarberry [6] and an equivalent version can
be found in Devroye [4]. A similar algorithm has also been given in [13] (Ex. 4.27 and page
422) in a different context. An informal description of the algorithm is as follows. (For the
sake of convenience, suppose that log 2 n is an integer in this discussion.) The algorithm is
based on a binary tree with 1 levels and n leaf nodes at the bottom level. The value
associated with the ith leaf node is A i . The value associated with any other node is the sum
of the values associated with the two sons of that node. Thus, the value of the node at the
root of the binary tree will be Y := A 1 . Note that O(n) preprocessing time is
needed to set up the binary tree. Now to generate a variate for X, we first generate a uniform
number, call it Z, over (0; Y ). We then compare this number with the value associated with
the left child of the root. If it is less, we know that the variate is in f1; so we
proceed with the algorithm on the left side of the binary tree. If it is more, we proceed on
the right side of the tree with a value Z minus the value stored in the left child. Note that
the number of comparisons needed to generate a variate for X is log 2 n. If an A i changes to
after generating a variate for X, we can update the binary tree with O(log n) operations
as follows. We first reset the value associated with leaf node i to A 0
i . We then move up the
path between this node and the root and reset the sums accordingly.
In summary, the above 'binary-tree algorithm' requires O(n) preprocessing time and
O(log n) time per variate to generate variates for X with changing distributions. The advantage
of this algorithm, as compared to Generate2, is that its performance does not depend
on how close the A i 's are to their upper bounds. Its disadvantage is that it can take significantly
more time to generate a variate for X when n is large and the acceptance probability
for Generate2 is not small.
6.1 Method I: Partitioning
For any subset S of
d
and
Y (S) :=
Now suppose there exists a partition (S 1 ng with the following properties: (i)
is in the 'vicinity of' or larger than
Under the above conditions, the following scheme makes sense. First draw a uniform
random variate U . If U then we declare that the variate for X
belongs to S 1 and we use Generate2 (across S 1 ) to determine it. If U - Y (S 1 )=[Y (S 1
then we declare that the variate for X belongs to S 2 and we use the binary-tree
method (across S 2 ) to find it. Of course, in order to utilize this method, we need to get a
handle on E[A i n. This can perhaps be done with analytical analysis or with
pilot runs.
In the context of queueing networks, the above method may be suitable when a fraction
of the queues are in heavy traffic and the remaining queues are in light traffic. We mention
here that another partitioning scheme was noted in [7].
6.2 Method II: Pseudo Upper Bounds
Here, for those indices i such that E[A i ] !! - a i , we replace - a i with a smaller value, perhaps
with E[A i is the long-run standard deviation of A i (s). (This, of course,
assumes that one can get a handle on well as on E[A i ]). Now we may have A i ? - a i
for some indices i. When this occurs, the algorithm Generate2 is no longer correct since we
may have B[i] ? 1 for some i. In order to rectify the algorithm, we keep track of the set
and of the counter ffi := maxfA i \Phig. (If \Phi is empty, set
Then in Step 3 of Generate2 we replace the test 'R - B[I]' with the test `R - B[I]=ffi'. We
leave it to the reader to verify the correctness of the procedure.
But when is it necessary to update \Phi and ffi, and how much effort is required for each
update? Suppose that an A i changes to an A 0
i . It is only necessary to perform an update
in the following circumstances: (i) A i =-a
a
In order to minimize the effort
to update \Phi and ffi, \Phi can be implemented as a priority queue. With such a data structure,
for any one of the four events occurs, O(log j\Phij) operations are sufficient for the update.
The pseudo upper bounds should be chosen large enough so that j\Phij and ffi are typically
small. However, they should not be chosen so small that the acceptance probability becomes
undesirably low.
6.3 Method III: Global Updates
For Method II it may not be possible to determine, a priori, good pseudo upper bounds - a i ,
may be the case that the appropriate choice of pseudo upper bounds changes
with the evolution of the underlying process that we are simulating. In these cases we may
want to consider occasional global updates of our pseudo upper bounds and reinitializations
of the arrays C and B. The reinitialization will make sure that all the entries in B are close
to 1. Of course, the process of global update is very costly (requiring O(n) time). But if this
has to be done very rarely, it may be advantageous.
For when compared to -a
a n . It may
be worthwhile to perform a global update. But how do we detect at any given time if Y is
low or not? When Y is small, many entries in the array B will be much less than 1. We
can make use of the following sampling process: We pick say 10 log n random outcomes. If a
major fraction of these outcomes have an entry much less than (say) 1=2 in the B array, we
perform a global update. Using Chernoff bounds we can show that if a major fraction in the
sample has a low B[] value, then a major fraction of all the B[] values will be low with high
probability. Also, we perform this checking only every 20 log n samples (thus making sure
that the total additional cost per sample due to this checking is no more than a fraction).
A similar reinitialization was mentioned in [8].
Conclusions
In this paper we have presented three simple algorithms for generating a nonuniform discrete
random variate with changing distributions. Under fairly unrestrictive assumptions, these
algorithms have an expected O(1) run time. Simulation results [14] confirm the competitiveness
of our algorithms in relation to future event schedule and other algorithms now widely
in use.
Though our algorithms are useful for many applications, an important open problem is
as follows. Does there exist a constant time algorithm for discrete random variate generation
which does not make any assumptions on the way the rates of X(s) change? Our algorithms
are a positive step in this direction.

Acknowledgements

We are grateful to Professor George Fishman, Professor Bennet Fox, Professor Kurt
Mehlhorn and Dr. Marty Reimann for their valuable comments.



--R

'Fast Probabilistic Algorithms for Hamiltonian Circuits and Matchings,' Journal of Computer Systems and Science 18(2)
A Guide to Simulation
'A Measure of Asymptotic Efficiency for Tests of a Hypothesis Based on the Sum of Observations,' Annals of Mathematical Statistics 23
'Non-Uniform Random Variate Generation,' Springer
'Exploiting Special Structure to Improve Future Event Set Management in Simulation,' Technical Report No.

'Shortening Future Event Lists,' to appear in ORSA J.
'Generating Markov-Chain Transitions Quickly: I,' Operations Research Society of America Journal on Computing
'Generating Markov-Chain Transitions Quickly: II,' Operations Research Society of America Journal on Computing

Queueing Systems
'Computer Performance Evaluation Methodol- ogy,' IEEE Transactions on Computers 33
Introduction to Algorithms: A Creative Approach
'Efficient Simulation of Large-Scale Loss Networks,' Technical Report MS-CIS-91-62
'Derivation of Randomized Sorting and Selection Algo- rithms,' Technical Report
An Introduction to Queueing Networks
--TR
A guide to simulation (2nd ed.)
Introduction to Algorithms

--CTR
Paul B. Callahan, Output-sensitive generation of random events, Proceedings of the ninth annual ACM-SIAM symposium on Discrete algorithms, p.374-383, January 25-27, 1998, San Francisco, California, United States
Keith W. Ross , Danny H. K. Tsang , Jie Wang, Monte Carlo summation and integration applied to multiclass queuing networks, Journal of the ACM (JACM), v.41 n.6, p.1110-1135, Nov. 1994
randomized algorithms;queuing networks