Likelihood Ratio Derivative Estimation for Finite-Time Performance Measures in Generalized Semi-Markov Processes.
This paper investigates the likelihood ratio method for estimating derivatives of finite-time performance measures in generalized semi-Markov processes (GSMPs). We develop readily verifiable conditions for the applicability of this method. Our conditions mainly place restrictions on the basic building blocks (i.e., the transition probabilities, the distribution and density functions of the event lifetimes, and the initial distribution) of the GSMP, which is in contrast to the structural conditions needed for infinitesimal perturbation analysis. We explicitly show that our conditions hold in many practical settings, and in particular, for large classes of queueing and reliability models. One intermediate result we obtain in this study, which is of independent value, is to formally show that the random variable representing the number of occurring events in a GSMP in a finite time horizon, has finite exponential moments in a neighborhood of zero.
Introduction
When running a simulation, one is often interested in estimating derivatives of a performance
measure with respect to parameters of the input distributions (e.g., routing probabilities,
time distributions and service time distributions in queueing networks; lifetime
distributions and repair time distributions of components in reliability models). One reason
for computing derivative estimates is that they are useful for sensitivity analysis. For example,
in the reliability context, a designer of a system may want to determine where to focus design
efforts to improve overall system performance. In addition, efficient derivative estimation
plays an important role in simulation-based optimization algorithms (e.g., see Glynn 1986,
1989a).
One approach for estimating derivatives via simulation is the likelihood ratio method (also
called the score function method); e.g., see Aleksandrov, Sysoyev, and Shemeneva (1968),
Arsham et al. (1989), Glynn (1986,1989a,1990), Reiman and Weiss (1986,1989), Rubinstein
(1986,1989), L'Ecuyer (1990,1995), Nakayama, Goyal, and Glynn (1994), and Nakayama
(1995). Overviews of this technique may be found in Bratley, Fox and Schrage (1987),
Glynn (1987), L'Ecuyer (1991) and Rubinstein and Shapiro (1993). The validity of this technique
has been established explicitly for certain performance measures of classes of discrete-
space, discrete-time Markov chains (Glynn 1986), discrete-space, continuous-time Markov
chains (Reiman and Weiss 1989, and Nakayama, Goyal, and Glynn 1994), and general-state-
space Markov chains that have a regenerative structure, also called Harris chains (Glasserman
and Glynn 1992). The performance measures that have been studied have mainly been transient
ones: performance measures that are based on the behavior of the stochastic process
in some time interval from the beginning, the length of which is either fixed or random (a
stopping-time). Steady-state performance measures have been studied insofar as they can
be expressed as functions of stopping-time based transient performance measures using some
regenerative structure. The main impediment in verifying the validity lies in justifying the interchange
of derivative and expectation that is required in this method. L'Ecuyer (1990,1995),
Reiman and Weiss (1989) and Rubinstein and Shapiro (1993) have general conditions that
ensure when this interchange is justified. However, these conditions are difficult to verify in
many specific contexts. In this paper, we establish conditions that are readily verifiable in
many contexts for the validity of the likelihood ratio method in the setting of generalized
semi-Markov processes (GSMPs).
A GSMP is a general mathematical framework for modeling many discrete-event systems,
and there has been a considerable amount of work studying different aspects of this class
of stochastic processes; e.g., see Shassberger (1976), Whitt (1980), Glynn (1989b), Haas and
Shedler (1987), and Glasserman and Yao (1992a, 1992b). Also, because their dynamics closely
follow those of event-driven simulations, GSMPs have been useful in modeling systems analyzed
through simulation and for studying various simulation methodologies; e.g., see Glynn
and Iglehart (1988), Glasserman (1991a), and Glasserman and Yao (1992c). In this paper we
use the GSMP framework to study the likelihood ratio derivative estimation method applied to
finite-time performance measures, i.e., performance measures that are expectations of random
variables that can be determined by a fixed time t. The generality of the GSMP framework
enables us to verify the applicability of the likelihood ratio derivative estimation method for
large classes of reliability and queueing models. In fact, this whole project was initially motivated
by the need to do simulation-based derivative estimation in a large class of reliability
models of the type in Heildelberger, Shahabuddin, Nicola (1994), Nicola, Nakayama, Heildel-
beger and Goyal (1993) and Nicola, Shahabuddin, Heidelberger, Glynn (1993). These are
non-Markovian versions of models constructed and analyzed by the SAVE software package
(Blum et al. 1993).
An alternative method for estimating derivatives using simulation is infinitesimal perturbation
analysis (IPA); e.g., see Ho and Cao (1983), Heidelberger et al. (1988), and Glasserman
(1991a). Glasserman (1991a,1991b,1991c) established conditions under which this method
will give rise to unbiased estimates of derivatives of finite-time performance measures of
GSMPs. One main difference between the conditions for IPA and our conditions for the
likelihood ratio method is that ours impose restrictions on the basic building blocks (i.e., the
transition probabilities, the distribution and density functions of the event lifetimes, and the
initial distribution) of the GSMP, whereas those for IPA relate to the underlying structure of
the GSMP. We give examples of large classes of stochastic models for which one can justify
the validity of the likelihood ratio method but not IPA, due to the limitations posed by the
structural conditions. On the other hand, the IPA method admits a larger class of event-time
distributions than the likelihood ratio method. Also, the variance of likelihood ratio method
derivative estimators grows quite fast as the length of the observation increases; e.g., see
Glynn (1987) for results in the setting of discrete-time Markov chains. Thus, for large time
horizons t, the likelihood ratio method is statistically inefficient for estimating derivatives.
IPA derivative estimators do not suffer from this drawback.
The rest of the paper is organized as follows. Section 2 contains a description of the GSMP
model with which we will work. Section 3 reviews the basic idea of the likelihood ratio method
for estimating derivatives. The section also provides a proof of the validity of the method in
the GSMP setting under a set of technical conditions. One intermediate result that is derived
here, shows that the random variable representing the number of occurring events in the
GSMP in [0; t], has finite exponential moments around the neighborhood of zero. This is
done by bounding the GSMP by an age-dependent branching process, and using results on
finite exponential moments of such branching processes. Section 4 gives some simple sufficient
conditions for our conditions in Section 3 to hold. In particular, we show that the conditions
hold in many settings arising in practice. Finally, in Section 5 we explicitly demonstrate that
the conditions in Section 3 hold for certain classes of reliability models and queueing network
models.
Mathematical Model
We model the evolution of our stochastic system as a generalized semi-Markov process (GSMP),
which we now describe. Our development of the model closely follows Glasserman's (1991a)
approach. We denote the (physical) state space of the process by S, where we assume that
S is either finite or countably infinite. For example, in the context of multi-class queueing
networks, a state in S contains information about the number of customers of each type at
each station, along with any description of the queueing at the various stations. We let A
denote the set of events. In the example of multi-class queueing networks, the events are the
completion of service and external customer arrival at each station. We make the following
assumption:
A1 The set A is finite.
A case where this assumption is not satisfied is the infinite-server queue. For a state s 2 S,
we define E(s) as the set of active events in state s; i.e., these are the events that may occur
when the system is in state s. In our queueing-system example, a service-completion event at
a particular queue is only possible in states in which the server at that queue is busy.
We introduce a scalar parameter ' with respect to which derivatives are computed. For
example, ' might represent the mean arrival time of customers or a routing probability in a
queueing network. After computing the derivative, we then want to evaluate the derivative at
the parameter value and \Theta is some open set. We can easily extend our
setting of calculating a derivative with respect to a scalar-valued parameter ' to work with
computing the gradient with respect to a vector-valued parameter
considering each ' (i) separately.
For each event e 2 A, we let F ( denote the distribution used to generate event
times for e under parameter value '. We assume that each distribution function F (
has a density f( (with respect to Lebesgue measure) and that F (0;
initial state is chosen using the probability mass function -( \Delta ; '), and given that a state s is
selected as the initial state, we generate event times for each event e 2 E(s) from F (
For each ' 2 \Theta, s 2 S, and e 2 A, let p( \Delta ; s; e; ') be a probability mass function on S. Under
parameter value ', if the system is currently in a state s 2 S and some event e 2 E(s) occurs,
then the system immediately moves to a new state s 0 2 S with probability p(s
the are called the transition probabilities. (The general
theory of GSMPs allows for more generality; e.g., see Glynn 1989b.)
For any state s 2 S, we define the set of possible clock readings as
We identify with an jAj-dimensional vector that we write as c. When the system
is in state s, the clock c(e) for any active event e 2 E(s) is run down at rate 0 -(s; e) ! 1.
We assume that -(s; e) ? 0 for some e 2 E(s). In most applications -(s;
letting -(s; e) assume values other than 1 allows greater modeling capabilities. For example,
we can model pre-emptive service by setting -(s; is a service-completion event that
is pre-empted in state s. We assume that the clock rates are bounded; i.e.,
There exists some constant -   ! 1 such that the clock rates -(s; e) -   for all s 2 S
and all e 2 A.
In most applications, A2 holds, but one example where it does not is an infinite-capacity open
queueing system in which the rate at which the server processes customers is proportional to
the queue length.
We assume that in each state s 2 S, there is at least one active event e 2 E(s). Given that
the system is currently in state s with clock vector c, we define
as the holding time in the current state. Also, we define
which is the event triggering the transition out of the current state. (Since all of the clock
times have continuous distributions, two clocks are never equal with probability one for all
parameter values '.)
For states s; s 0 2 S and an event e 2 A, let C(s denote the set of events that are active
in state s but are cancelled upon the transition to state s 0 triggered by the event e. Then
we define the set of new events N (s and the set of old events O(s
respectively. After a
transition from state s to state s 0 triggered by event e, for each (new) event e 0
a new clock value c(e 0 ) is generated from the distribution F ( Also, for each (old)
event we update its new clock value to be c(e immediately
after entering state s 0 . Each (cancelled) event e 0 2 C(s originally active in state s
but may no longer be active in state s 0 . Note that we allow for the possibility that a cancelled
event e 0 is immediately reactivated (i.e., e e)), in which case a
new clock value is immediately generated for the event e 0 .
We define the discrete-time Markov process f(Y is the n-th state
visited and c n is the vector of clock readings just after the n-th transition. The state space of
the process is
s2S (fsg \Theta C(s)); and its transition probability is as follows: for A ae \Sigma of
the form fixed (possibly infinite) real
numbers the transition probability from (s; c) to A under parameter value ' is
-((s; c); A;
Y
Y
and 1fBg is the indicator function of the probability
event fBg.
We define our continuous-time GSMP as a piecewise-constant process that stays in state Y n
for an interval of length t   (Y precisely, let -
which is the epoch of the n-th transition. For notational simplicity, define e
which is the number of transitions that
occur by time t - 0. Then, we define our GSMP
Let P ' denote the probability measure governing our GSMP under parameter value ', and
be the expectation operator induced by P ' .
space of the GSMP, and let F t denote the filtration of the GSMP up to a (deterministic)
loosely speaking, the history of the process up to time t); e.g., see Section 36
of Billingsley (1986).
Let ff(t; our performance measure, where V (t; ') is some F t -measurable
random variable. We will be interested in estimating the derivative d
d' ff(t; ') at the parameter
be the probability measure of the GSMP when we fix
is
the expectation operator induced by P ' 0
We now impose some restrictions on the random variable V (t; '):
A3 There exists some constant h ? 0 such that V 0 (t; '; !) j d
d'
-almost
2\Omega and all ' 2 (' h), and there exist random variables W and W 0 such that
-probability 1
for all ' 2 ('
We show in Section 4 that many performance measures of interest satisfy the above condition.
3 Likelihood Ratio Derivative Method
In this section we review the likelihood ratio method for estimating d
d'
d'
which is to be evaluated at the parameter value To compute d
like to bring the derivative operator inside the expectation. However, the expectation operator
depends on ', and so some care must be taken. Assuming that P ' is absolutely continuous
with respect to P ' 0
restricted to the filtration F t , we have that
Z
Z
where L(t;
is the Radon-Nykodym derivative of P ' with respect to P '0 restricted
to the filtration F t , or simply the likelihood ratio. absolutely
continuous with respect to another probability measure Q 2 if the sets of Q 2 -measure 0 are
also of Q 1 -measure 0; e.g., see Section 32 of Billingsley 1986.) This is known as a "change
of measure," and Glynn (1989) and Damerdji (1994) established its validity for GSMPs in
certain contexts. Thus we have expressed the performance measure ff(t; ') as an expectation
of some random quantity taken with respect to a probability measure that is independent of
the parameter '. The key question that remains now is whether we can bring the derivative
operator inside the expectation; i.e., whether
d
d'
d
d'
Assumption A1 of L'Ecuyer (1990) and Assumption A1 0 (k) of L'Ecuyer (1995) (see also
Rubinstein and Shapiro 1993) provide conditions when (2) holds on a general probability
space. However, in the specific context of the GSMP, these are difficult to verify directly.
Before answering the interchange question in the GSMP setting, let us first determine the
exact form of the likelihood ratio L(t; '). Heuristically, the likelihood ratio L(t; ') is given by
Y
Y
see Glynn and Iglehart (1989) for further details. However, if the above version of the likelihood
ratio is used, we would need very restrictive conditions on the distribution functions F
to justify the interchange of derivative and expectation. (For example, we would have to
change A5 below to alternatively require f( uniformly on [0; 1) as
which is not even satisfied when F is the exponential distribution with mean
1='.) To get around this, we instead work with a slightly modified GSMP, which leads to
a different likelihood ratio but does not alter the distribution of V ('; t). We do this by
changing the distributions from which clock times are generated: for all e 2 A, define the new
distribution function
where we recall that t is the time horizon. Note that F   (x; e; ') is a mixed (continuous and
discrete) distribution that is the same as F (x; e; ') for x -   t, and has a point mass of
probability -
F (-   t; e; ') at the point -
before, now more than one event can occur at the same time. We can use any rule we want to
deal with this. For example, we can order the set of events in some manner, and if more than
one event occurs at the same time, we take the "lowest" of these events as the "occurring"
event (for determining the next physical state visited, etc. In any case, both the original
and the modified GSMP (i.e., the process Z defined in the previous section) have the same
distribution of sample paths in [0; t]. Since V ('; t) is only dependent on the sample path of the
GSMP from time 0 to t (i.e., it is F t -measurable), using the new distributions for generating
clock times does not change the distribution of V ('; t). By applying this approach, we will
only require mild conditions on the distribution functions F to prove the interchange.
The new likelihood ratio is given by
Y
Y
Y
Y
Now assuming that (2) holds, we see that applying the product rule of differentiation yields
d
d'
where
with
d'
d'
d'
d'
gives the derivative of a likelihood ratio expression
similar to (3), but for the case of a fixed number of transitions (in contrast to a fixed time
horizon). Also, Glynn did not establish the validity of the interchange of the derivative and
expectation in the likelihood ratio derivative method applied to GSMPs.
Now we develop conditions on the basic building blocks of the GSMP under which (2) is
valid. (We later give in Section 4 some easily verifiable sufficient conditions for the following
conditions to hold.) First, we impose some regularity conditions on the densities generating
the event lifetimes.
A4 For each event e 2 A, there exists h ? 0 such that the set   f (e)
independent of ' for ' 2 ('
A5 For each event e 2 A such that   f (e) 6= ;, for each ffl ? 0 there exists h ? 0 such that
uniformly on   f (e) as
A6 For each event e 2 A such that   f (e) 6= ;, there exists h ? 0 such that f( \Delta ; e; ') is
differentiable in ' on   f (e) for all ' 2 (' \Delta) is uniformly
bounded on   f (e) \Theta ('
Note that Condition A4 does not permit the support of the distribution of the lifetime
of an event e to depend on '. Thus, for example, we disallow an event e to have a uniform
distribution with support (0; ').
We will now present some properties of -
are implied by the previous conditions
on f( The proof of this lemma may be found in the Appendix. Let   -
Assume that Condition A4-A6 hold. Then
(i) there exists h ? 0, such that   -
F (e) is independent of ' for ' 2 ('
(ii) for each event e 2 A,
uniformly on [0; -   t] as
(iii) for each event e 2 A such that -
0 such that
(iv) for each event e 2 A such that -
is differentiable in ' for all ' 2 ('
F (-   t; e; \Delta; ) is bounded on
We now impose some regularity conditions on the initial distribution.
A7 There exists h ? 0 such that the set   - independent of ' for
A8 For each ffl ? 0 there exists h ? 0 such that
uniformly on   - as
There exists h ? 0 such that -( \Delta ; ') is differentiable in ' on   - for all ' 2 ('
uniformly bounded on   - \Theta ('
Finally, we impose some regularity conditions on the transition probabilities.
There exists h ? 0 such that the set
0g is independent of ' for ' 2 ('
A11 For each ffl ? 0 there exists h ? 0 such that
all uniformly on   p
as
There exists h ? 0 such that p( \Delta ; \Delta; \Delta; ') is differentiable in ' on   p for all ' 2 ('
\Delta) is uniformly bounded on   p \Theta ('
In some sense, Conditions A4-A12 ensure that the basic building blocks of our GSMP are
well-behaved functions of ' in a neighborhood of ' 0 . Now we establish the validity of (2).
Theorem 1 If Conditions A1-A12 hold, then
[V (t; ')L(t; ')] is differentiable and (2) is
valid.
To show that Theorem 1 holds, we will first need to prove the following result which is of
independent value.
Conditions A1, A2, A4 and A5 hold. Then there exists constants z 0 ? 0
Proof. Observe that N(t) has the same distribution whether we use F   (
to generate lifetimes for the event e 2 A. Hence, we work with F (
For the case of systems with no event cancellation, we can show the result using the
methods employed in the proofs of Corollary 3.3 of Glasserman (1991a) and Theorem 2.1 on
page 155 of Prabhu (1965), that use results from renewal theory. For systems with event
cancellation, we need a new approach, which is given below and uses results from the theory
of branching processes.
We will establish the result by first constructing a new process that bounds N(t). Then
we will show that this new process has the same distribution as a multi-type, age-dependent
branching process (see, e.g., Athreya and Ney (1972) for a definition). We then bound this
multi-type age-dependent branching process by a single-type, age-dependent branching pro-
cess, and finally use a result of Nakayama, Shahabuddin, and Sigman (1996). For simplicity,
we will assume that all of the clock rates -(s; 1. The proof can be easily generalized
when the clock rates are unequal and bounded, as in Condition A2. For each e 2 A, recall
that lifetimes of the event e follow a distribution F ( '). Some of the lifetimes may not
expire in the GSMP because of event cancellation.
In what follows, we frequently will be using the following multi-type, age-dependent
branching process. To simplify notation, let 1. The
branching process has m different types of objects, indexed by e. Objects of type e have
lifetime distributions F ( '). The branching process begins at time 0, at which time
one object of each type is born, i.e., a total of m objects are born. Each object lives for a
random time (specified by its corresponding lifetime distribution) and then dies. When any
object dies, it instantaneously gives birth to exactly one object of each type. The process thus
continues.
We now construct a sample path of our first bounding process as follows. We begin by
generating a sample path of the original GSMP in [0; t]. Then we do the following:
ffl In the sample-path generation of the original GSMP, at each transition point (equivalent
to an event expiration point of a non-cancelled event, but also includes
of events (which may also be a null-set) in A are scheduled. In the new process, at each
transition point of the GSMP we schedule each of the other events in A that have not
already been scheduled at that transition point. These new events are never cancelled
and whenever any of them expire, we start a multi-type branching process of the sort
mentioned before.
ffl In the original GSMP, events are cancelled. In the new process, we let each of the
cancelled events run to expiration. And whenever any of them expires, we again start a
multi-type branching process of the kind mentioned before.
Let B(t) be the total number of events scheduled in [0; t] in this constructed process.
We now claim that this constructed
process is a sample path of a multi-type, age-dependent branching process.
On first glance this claim may seem incorrect. This is because the event expiration times
of cancellable events (if not cancelled) depend on the event expiration times of events that
cause the cancellation (by their expiration). Then the independent lifetimes assumption of
a branching process seems to be violated. However, a careful argument reveals otherwise.
Note that the only difference between the constructed process and an ordinarily generated,
multi-type, age-dependent branching process is the order of generation of the lifetimes. In the
ordinary case, the lifetimes may have been generated sequentially, i.e., we first generate the
lifetimes of each of the m objects in the first generation, then of each of the m 2 objects in
the second generation, and so on. In the constructed process, the order in which we generate
lifetimes is random. Not only is it random, but the next lifetime one chooses to generate
depends on the position and length of the lifetimes already generated. However, the lifetimes
themselves are independent of anything in the system and hence the order of generation will
not change the distribution of the branching process.
Now consider a new branching process defined as follows. The branching process begins at
time 0, at which time m objects are born. When any object dies, it gives birth to exactly m
objects. The lifetime of each object follows a distribution H, which depends only on ' 0 . We
define it as follows. First, define a distribution function G with G(x;
for x - 0, and 0 otherwise. Note that if X e has distribution F (
using Condition A4, A5, and part(ii) of
Lemma 1 we can easily show that G(x; uniformly on the set [0; t], as
Hence, we can choose a sufficiently small h ? 0 such that jG(x;
note that ffl - 1=(2m). Now define the new distribution function H as follows: H(x; ' 0
otherwise. Observe that the distribution H( stochastically smaller (see Ross 1983)
than G( which implies that H( stochastically smaller
than
for all ' 2 (' because we assumed that F (0;
Note that the above defines a single-type, age-dependent branching process. Let B 0 (t) be
the number of objects born up to and including time t, Also, since H( stochastically
smaller than G( we have that B 0 (t) is stochastically larger
than B(t) for all t - 0 and ' 2 ('
Finally, since H(0; ' 0 )m ! 1, we can use Theorem 3.1 of
Nakayama, Shahabuddin, and Sigman (1996) to show that there exists some z 0 ? 0 such that
so the proof is complete.
Proof of Theorem 1.
We need to show that for all sufficiently small ffi, the absolute values of the difference
quotients are dominated uniformly by a P '0 -
integrable random variable. By the mean value theorem, for each sufficiently small ffi, there
exists
which is finite since jAj ! 1. Choose ffl ? 0 such that
Clearly, by Lemma 2, this is possible. Now choose a sufficiently small h ? 0 so that Conditions
A3, A4, A6, A7, A9, A10, A12 and part (iv) of Lemma 1 hold; Conditions A5, A8, A11
and part (iii) of Lemma 1 hold for the above ffl; and (7) holds for all jffij - h. Hence in (7), the
j's corresponding to these ffi's always lie in ('
First consider the term V (t; j)L 0 (t; j). From (5),
Y
Y
F (-   t; e;
Y
Y
F (-   t; e;
\Delta6 6 4
F (-   t; e;
All of the density functions f( in the above expression are evaluated at points no
greater than -   t. Now we define
sup
sup
sup
By Conditions A1, A4-A12, and part(iii) and (iv) of Lemma 1, M 2 (h), M 4 (h), M 6 (h), and
8 (h) are finite, and M 1 (h), M 3 (h), M 5 (h), and M 7 (h) lie in the interval
Then, letting OE 1
with
-probability 1 by Condition A3. Repeated applications of the Schwarz inequality
yields
In addition, using the fact that OE 1 we have that
implies that all moments of N(t) exist under P ' 0
and hence
Hence, OE(h) is
-integrable. Similarly, we can show that jV 0 (t; j)L(t; j)j is uniformly dominated
-integrable random variable. Finally, applying the dominated
convergence theorem completes the proof.
It will be worthwhile comparing our conditions with similar ones in the literature. First,
consider the "amiability" condition used in Theorem 3 of Reiman and Weiss (1989). One may
use Theorem 3 of Reiman and Weiss to validate the exchange of derivative and expectation
only in special cases of our setup, i.e., where ' is present only in one distribution and it
never appears as a variable in the performance measure. However, even if we impose these
restrictions in our setting, the amiability condition is difficult to verify directly. One can
show, only with the help of Lemma 2 and a mean value theorem, that our assumptions (with
Condition A3 replaced by the weaker and more tractable assumption used in Proposition 1
below) imply the amiability condition.
Secondly, consider assumption (1995). Though the setting
used there is different from ours, we can still make some comparisons. For example, our
Condition A3 is similar to his conditions on the performance measure. Also, in some sense,
our type of distributional conditions imply the distributional conditions in A3(k) of L'Ecuyer;
i.e, if our type of distributional conditions hold, then one can use a constant close to 1 as his
bounding random variable \Gamma 1i (i) and some constant as his bounding random variable \Gamma 2i (i).
If we use N(t) as an analogue to his discrete stopping time - , then an analogous result to
Lemma 2 given below will prove that the conditions specified in A3(k) on \Gamma 1i (i) and \Gamma 2i (i)
hold.
We now compare our conditions to those employed by Glasserman (1991a,1991b) to show
the unbiasedness of IPA derivative estimators in the GSMP context. First, note that our
conditions mainly impose restrictions on the basic building blocks of the GSMP. On the other
hand, the conditions for IPA are primarily on the structure of the GSMP. In particular,
Glasserman requires a "commuting condition" which places limitations on the relationship
between active events in certain states and the transition probabilities from those states and
from possible immediate successor states. The condition essentially ensures that if slightly
altering the value of the parameter ' results in a change in the order of two events on a sample
path, then the sample path can (in some sense) correct itself on the next event. Moreover,
Glasserman does not permit the initial distribution nor the transition probabilities to depend
on the parameter '.
Another difference is the type of performance measure ff(t; In
our results, we restrict V (t; ') to be any F t -measurable random variable that satisfies Condition
A3, where t is deterministic. Glasserman allows V (t; ') to be an additive functional up
to a time T , but T is allowed to be either deterministic, the time of the n-th transition, or
the time of the n-th occurrence of a particular event.
Sufficient Conditions for A1-A12 to Hold
We now show that our assumptions hold in many settings arising in practice. (In Section 5
we explicitly demonstrate that the following sufficient conditions are satisfied for large classes
of reliability and queueing models.) First consider Condition A3.
Proposition 1 Suppose Conditions A1, A2, A4 and A5 hold. Also, suppose ff(t;
is an F t -measurable random variable such that V 0 (t; '; !) exists for P '0 -almost
2\Omega and all ' 2 (' that satisfies the following:
(i) there exist finite constants
-probability 1 for all ' 2 ('
(ii) there exist finite constants K 3
with
-probability 1 for all ' 2 ('
Then Condition A3 holds.
Proof. We need to show that W j K 1
second moments when the parameter value . But this follows from Lemma 2.
We now show directly that many performance measures arising in practice satisfy Condition
A3.
Proposition 2 Suppose the performance measure ff(t; one of the
real-valued function on S with g(s; \Delta )
continuously differentiable for each s 2 S, and either (a) jSj ! 1, or (b) g(
and are uniformly bounded on S \Theta ('
d'
Then Condition A3 holds.
Note that V 0 (t;
has form (ii).
Proof. First note that V (t; ') in (i) and (ii) are F t -measurable. When condition (i) holds,
trivially satisfies Condition A3. Now consider form (ii). Assuming
that either condition (ii)(a) or (ii)(b) holds and since g(s; \Delta ) is continuously differentiable,
we have that for sufficiently small h ? 0, -
t, which satisfies Condition A3.
Now we show that two rich classes of distributions, the Gamma distribution and the
Weibull distribution, with the ' in each case being the scale parameter, satisfy the distributional
conditions given in A4-A6. The Gamma distribution is widely used in all types of
stochastic modeling; the Weibull distribution is one of the most frequently employed distributions
in reliability modeling.
Proposition 3 Suppose the lifetime of an event e has a gamma density with shape parameter
a ? 0 and scale parameter ' ? 0; i.e., f(x;
for the given e,
(i) the set fx independent of ' for ' ? 0;
uniformly on [0; t] as
uniformly bounded on [0; t] \Theta (' 0 =2; 3' 0 =2) .
Proof. Part (i) is trivial. For part (ii), note that f(x; e; ')=f(x;
which converges to 1 uniformly on [0; t] as
For part (iii), note that f 0 (x;
a
2a
for all x 2 [0; t] and all ' 2 (' 0 =2; 3' 0 =2).
Proposition 4 Suppose the lifetime of an event e has a Weibull density with shape parameter
a ? 0 and scale parameter ' ? 0; i.e., f(x;
Then for the given e,
(i) the set fx independent of ' for ' ? 0;
uniformly on [0; t] as
uniformly bounded on [0; t] \Theta (' 0 =2; 3' 0 =2) .
We omit the proof since it is similar to that of Proposition 3.
Note that we do not need to consider any distributions that do not depend on the parameter
' when checking if Conditions A4-A6 are satisfied. The same is true for the initial
distribution -( Conditions A7-A9, and the transition probabilities p(
Conditions A10-A12.
The following result shows that if the state space is finite and the convergence and bounds
in Conditions A8, A9, A11, and A12 hold, then they are automatically uniform. The proof is
trivial and is therefore omitted.
Proposition 5 Suppose that jSj ! 1. Then
(i) Condition A8 holds if -(s; ')=-(s; '
(ii) Condition A9 holds if for each s 2 S, there exists some h ? 0 (which may depend on s)
such that - 0 (s; \Delta)=-(s; \Delta) is bounded on ('
(iii) Condition A11 holds if for each (s
(iv) Condition A12 holds if for each (s 0 ; s; e) 2   p , there exists some h ? 0 (which may
depend on s) such that p 0 is bounded on ('
5 Examples
We now give some examples of systems for which Conditions A1-A12 hold. The first example
deals with a large class of reliability models, and the second example considers a large class
of multiclass queueing networks.
5.1 Reliability Models
Consider the class of reliability models studied in Heidelberger, Shahabuddin, and Nicola (1994),
Nicola, Shahabuddin, Heidelberger, and Glynn (1993), and Nicola, Nakayama, Heidelberger,
and Goyal (1993). These are systems composed of K components, where each component can
fail and get repaired. Component lifetimes and repair times are generally distributed. There
are R classes of repairmen in the system, each consisiting of a finite number of repairmen. The
repair of each component is assigned to one repairman class. The repair strategy used by a
repairman class may be pre-emptive or non-pre-emptive priority with FCFS or random order
service used between members of the same priority class. We allow for the following interpen-
dencies among the components. First, the operation of an operational component may depend
on certain other components being operational. In addition, the repair of a failed component
may depend on certain other components being operational. Also, the failure of a component
may cause other components to fail instantaneously with certain probabilities. This is
called failure propagation, which is an example of how cancelled events arise. The system is
considered operational if certain sets of components are operational. The physical state of
the system after the nth transition is given by Y
n is 1 or 0 depending on whether the ith component is operational or failed after the
n-th transition, and Q (j)
n is the repair queue at the j-th repairman class. Note that jSj ! 1.
It is easy to see that the process fZ As we
previously mentioned, the system is considered operational if certain sets of components are
operational. More formally, U is the set of operational
states and F is the set of failed states. The system is considered operational at time t if
otherwise, it is considered down. We assume that at time
in the system are operational and new (i.e., all components have age = 0). Thus, the initial
distribution is independent of the parameter ', and so Conditions A7-A9 are automatically
The events in the system are either component failures or component repairs. Clearly,
1 as there are a maximum of 2K events in the entire system since each of the K
components can have a failure and repair event associated with it. Hence, Condition A2 is
satisfied since jSj ! 1. If there is no failure propagation in the system then the outcome of
any event e, given that the system is currently in some state s, is fixed; i.e., p(s
for exactly one state s 0 and 0 otherwise. However, if we have failure propagation in the system
then there may be some randomness in the final state reached from any given state after an
event, and Proposition 5 provides sufficient conditions for Conditions A11 and A12 to hold.
The transient performance measures in which one may be interested are the unreliability
and the expected interval unavailability. The unreliability is the probability of a system failure
before some fixed time t. Hence, in this case ff(t; is of form (i)
in Proposition 2 and A is the event fZ s 2 F for some s - tg. Thus, Condition A3 holds.
The interval unavailability is the fraction of time in [0; t] that the system is down. In this
case, ff(t; is of form (ii) in Proposition 2 with g(Z s
Fg=t. Note that g(Z s ; ') is independent of ' and therefore V 0 (t; ';
-probability 1, and so Condition A3 is satisfied.
In the reliability setting one is usually interested in the derivatives of the performance
measures with respect to (i) parameters of the failure time distribution of individual compo-
nents; (ii) parameters of the repair time distribution of individual components; (iii) the failure
propagation probabilities. For the first two cases, the likelihood ratio derivative estimation
method is applicable if the distributions that depend on the parameter with respect to which
the derivative is taken satisfy Conditions A4-A6. Some special cases where these conditions
are satisfied are when the distribution and the parameter are as in Propositions 3 and 4.
In the third case, the density functions f( are independent of the parameter ', and
so Conditions A4-A6 are automatically satisfied. Also, p(s typically will be of the form
are constants. Clearly, Condition A10 is satisfied for ' 2 (0; 1).
Since jSj and jAj are finite, we only have to check the conditions given in Proposition 5.
For the particular form of p(s these conditions are satisfied for all
For this class of reliability systems, it turns out that Glasserman's (1991a,1991b) conditions
for the unbiasedness of IPA derivatives typically are not satisfied. In particular, Glasserman
assumed that events are never cancelled, but this could occur in this setting since we allowed for
failure propagation. Also, his "commuting condition" will not hold for many repair strategies.
For example, if there is only one repairman who fixes all failed components in a first-come
first-served fashion, then this condition will be violated. Moreover, Glasserman's work does
not cover performance measures such as the unreliability. However, IPA can admit many more
event-time distributions. For example, one can use a repair time that is uniformly distributed
over (0; ') (the ' also being the parameter) in IPA, but it is not allowed by the results in this
paper.
5.2 Queueing Networks
Consider a multiclass queueing network where service times for each customer class at each
station are generally distributed random variables. We assume that there is a finite number
of servers at each station, and that each customer class uses Markovian routing; i.e., when
a customer of type i completes service at station j, he immediately goes to station k with
probability R i (j; k). We consider both open and closed networks. For the case of open
networks, we assume that the arrival process of a class at each station is a renewal process
with generally distributed inter-renewal times. The events in the system are customer arrivals
(one for each customer class at each station) and customer service completions (one for each
customer class for each server at each station). Hence, again the event set is finite.
In most applications Condition A2 is satisfied. As mentioned earlier, one situation in which
it is not is in an open network with the service rate of customers at a station is proportional
to the number of customers in the system.
If we assume that the network is closed and always starts in the same customer configura-
tion, then the initial distribution is independent of the parameter ', and so Conditions A7-A9
are automatically satisfied. Similarly, if we assume that the network is open and always starts
with no customers present, then the initial distribution is independent of the parameter ',
and Conditions A7-A9 automatically hold.
Some of the performance measures of interest in such systems are (i) the expected time-average
queue length over [0; t] at a given station; (ii) for the case of open networks, the
expected time-average number of customers over [0; t] in the entire system; (iii) for the case
of open networks, the probability that total number of customers in the system exceeds some
threshold in time interval [0; t] (if the network has a finite combined buffer, then this is the
probability that there is at least one customer loss due to buffer overflow in [0; t]).
We now show that Condition A3 holds for each of the above performance measures. Since
in all cases V (t; ') is independent of ' and V 0 (t; '; only have to show
that E '0 [(V (t; 1. For case (iii), this is obvious. Case (i) for closed networks is also
obvious. For case (i) for open networks and case (ii) and no customers are initially in the
system, this follows from Lemma 2 when Conditions A2, A4 and A5 hold. This is because,
say for case (ii), if we let Q(s) be the total number of customers in the system at time s, then
Q(s)ds -t
N(s)ds -t
and so Proposition 1 implies that Condition A3 holds under Conditions A2, A4 and A5.
One may be interested in the derivatives of the performance measures with respect to (i)
a parameter of the service time distribution of a class at a particular station; (ii) in the case
of open networks, a parameter of the interarrival-time distribution of a certain customer class
at a particular station; (iii) the routing probability of a certain customer class at a particular
station. As in the case with the reliability models, for the first two cases, the likelihood ratio
derivative estimation method is applicable if the distributions that depend on the parameter
with respect to which the derivative is taken satisfy Conditions A4-A6. Some widely used
cases when these conditions are satisfied are when the distribution and the parameter are of
the type given in Proposition 3 and Proposition 4.
In the third case, for the case of closed networks, due to the finiteness of the state space,
we can again apply Proposition 5. For the case of open networks, Proposition 5 is no longer
useful as now the state space is no longer finite. However, since the number of stations is finite
and since there is Markovian routing, there are only finitely many different values that the
assume. Hence, to show that Conditions A11 and A12 hold,
we only have to make sure that each of the transition probabilities is continuous at ' 0 , and
that each of the in some neighborhood of ' 0 . These
conditions hold if p(s is of the form c 1 are constants.
For the multi-class queueing systems considered here, Glasserman's (1991a,1991b) conditions
for unbiasedness of IPA derivatives typically are not satisfied. Specifically, his "com-
muting condition" will not hold for multi-class networks of our generality. For example, if
a particular station that is visited by more than one customer class is fed by more than
one source, then this condition is violated. Also, Glasserman's work does not cover performance
measures such as the probability of a buffer overflow occurring before time t. Finally,
Glasserman disallows the routing probabilities (more generally, the state-dependent transition
probabilities) to depend on the parameter '. However, as we mentioned before, IPA admits
many more event-time distributions.
We should mention that Conditions A1-A12 will hold for more general queueing networks.
For example, it can be shown that they typically are valid when there is state-dependent
routing. 1
A

Appendix


Proof of Lemma 1:
1 The work of the first author was supported in part by NSF CAREER Award DMI-96-24469 and by NJIT
SBR Grant No. 421180. The work of the second author was supported in part by NSF Grant DMS-95-08709
and by NSF CAREER Award DMI-96-25291. The authors would like to thank the area editor, associate
editor, and two referees for their helpful and detailed comments.
(i) Obvious.
(ii) We wish to show that for any ~ ffl ? 0 there exists ffi ? 0, such that for all ' satisfying
From Conditions A4 and A5, for any ~ ffl ? 0, there exists than the h is
Condition A4), such that for all ' satisfying
or
for all x 2   f (e). From Condition A4, for x 2 [0; -
all ' satisfying Therefore (10) still holds. Choose . Hence from (10), for
any ~ ffl ? 0 there exists ffi ? 0, such that for all ' satisfying we have that
or -
for all y 2 [0; -   t]. Subtracting -
from each term and using the fact F (y;
we get (9).
(iii) This follows immediately from part (ii) of the Lemma.
(iv) By Conditions A4 and A5, for any ~ffl ? 0, there exists than the h in
Condition A4) such that f(x; e; ') - f(x;
and x 2   f (e), define D(ffi; Conditions A4 and
A6, there exists a h 2 smaller than h 1 , such that for all ' 2 ('
exists and jf 0 (x;
ffl)Kf
some constant K ? 0. From part(iii) of Lemma 1, there exists h 3 such that for all ' 2
Choose the h in part (iv) of Lemma 1 to be min(h )=2. By the mean value theorem,
for all ' 2 (' as shown in the previous paragraph, f 0 (x; e; ') exists), for all
sufficiently small ffi there exists
all we have that
Also
Z
Then by the dominated convergence theorem, for all ' 2 (' 0 \Gammah; ' 0 +h),
R
R
R
R



--R


"Stochastic Optimization,"
"Sensitivity Analysis and the "
Probability and Measure
System Availability Estimator (SAVE) Language Reference and User's Manual
A Guide to Simulation (2nd
"Maximum Likelihood Ratio Estimation for Generalized Semi-Markov Pro- cesses,"
Gradient Estimation Via Perturbation Analysis
"Structural Conditions for Perturbation Analysis Derivative Estimates: Finite Time Performance Indices,"
"Structural Conditions for Perturbation Analysis of Queuing Systems,"
"Gradient Estimation for Regenerative Processes,"
and Wilson
"Monotonicity in Generalized Semi-Markov Processes,"
"Generalized Semi-Markov Processes: Antimatroid Structure and Second-Order Properties,"
"Some Guidelines and Guarantees for Common Random Numbers,"
"Stochastic Approximation for Monte Carlo Optimization,"
Likelihood ratio gradient estimation: an overview.
"Optimization of stochastic systems via simulation,"
IEEE Press
"A GSMP Formalism for Discrete-Event Systems,"
"Likelihood Ratio Derivative Estimators for Stochastic Systems,"
"Simulation Methods for Queues: An Overview,"
"Importance Sampling for Stochastic Simulations,"
"Regenerative Generalized Semi-Markov Processes,"
"Convergence Properties of Infinitesimal Perturbation Analysis Estimates,"
"Bounded Relative Error in Estimating Transient Measures of Highly Dependable Non-Markovian Systems,"
"Optimization and Perturbation Analysis of Queueing Networks,"


"Note: On the Interchange of Derivative and Expectation for Likelihood Ratio Derivative Estimators,"
"Asymptotics of Likelihood Ratio Derivative Estimators in Simulations of Highly Reliable Markovian Systems,"
"Likelihood Ratio Sensitivity Analysis for Markovian Models of Highly Dependable Systems,"
"A Note on Exponential Moments for
"Fast Simulation of Highly Dependable Systems with General Failure and Repair Processes.,"
"Fast Simulation of Steady State Availability in Non-Markovian Highly Dependable Systems,"

"Sensitivity Analysis via Likelihood Ratios,"
IEEE Press
"Sensitivity Analysis for Simulations via Likelihood Ratios,"
New York
"The Score Function Approach for Sensitivity Analysis of Computer Simulation Models,"
Sensitivity analysis and performance extrapolation for computer simulation models.
Discrete Event Systems: Sensitivity Analysis and Stochastic Optimization by the Score Function Method
"On the Equilibrium Distribution of a Class of Finite-State Generalized Semi-Markov Processes,"
"Continuity of Generalized Semi-Markov Processes,"
--TR

--CTR
Sigrn Andradttir , Paul Glasserman , Peter W. Glynn , Philip Heidelberger , Sandeep Juneja, Perwez Shahabuddin, 1962--2005: A professional appreciation, ACM Transactions on Modeling and Computer Simulation (TOMACS), v.17 n.2, p.6-es, April 2007
Pierre L'Ecuyer , Valrie Demers , Bruno Tuffin, Rare events, splitting, and quasi-Monte Carlo, ACM Transactions on Modeling and Computer Simulation (TOMACS), v.17 n.2, p.9-es, April 2007
gradient estimation;Generalized Semi-Markov Processes;simulation