Computing Densities for Markov Chains via Simulation.
We introduce a new class of density estimators, termed look-ahead density estimators, for performance measures associated with a Markov chain. Look-ahead density estimators are given for both transient and steady-state quantities. Look-ahead density estimators converge faster (especially in multidimensional problems) and empirically give visually superior results relative to more standard estimators, such as kernel density estimators. Several numerical examples that demonstrate the potential applicability of look-ahead density estimation are given.
Introduction
Visualization is becoming increasingly popular as a means of enhancing one's understanding of a stochastic
system. In particular, rather than just reporting the mean of a distribution, one often finds that more
useful conclusions may be drawn by seeing the density of the underlying random variable.
We will consider the problem of computing the densities of performance measures associated with
a Markov chain. For chains on a finite state space, this typically amounts to computing or estimating
a finite number of probabilities, and standard methods may be applied easily in this case (see below).
When the chain evolves on a general state space, however, the problem is not so straight-forward.
General state-space Markov chains arise naturally in the simulation of discrete-event systems (Hen-
derson and Glynn 1998). As a simple example, consider the customer waiting time in the single-server
queue with traffic intensity ae ! 1 (see Section 6). The sequence of customer waiting times forms a
Markov chain that evolves on the state space [0; 1). More generally, many discrete-event systems may
The research of the second author was supported by the U.S. Army Research Office under Contract No. DAAG55-97-
1-0377 and by the National Science Foundation under Grant No. DMS-9704732.
be described by a generalized semi-Markov process, and such processes can be viewed as Markov chains
on a general state space (see, e.g., Henderson and Glynn 1998). General state-space Markov chains are
also prevalent in the theory of control systems; see Chapter 2 of Meyn and Tweedie (1993).
This paper is an outgrowth of, and considerably extends, Glynn and Henderson (1998), in which
we introduced a new methodology for stationary density estimation. For a general overview of density
estimation from i.i.d. observations, see Prakasa Rao (1983), Devroye (1985) or Devroye (1987). Yakowitz
(1985), (1989) has considered the stationary density estimation problem for Markov chains on state space
ae IR d where the stationary distribution has a density with respect to Lebesgue measure. He showed
that under certain conditions, the kernel density estimator at any point x is asymptotically normally
distributed with error proportional to (nh d
is the so-called "bandwidth", and n is the
simulation runlength. One of the conditions needed to establish this result is that hn ! 0 as n ! 1.
Hence, the rate of convergence for kernel density estimators is typically strictly slower than n \Gamma1=2 , and
depends on the dimension d (see Remarks 5 and 7). In contrast, the estimator we propose converges at
rate n \Gamma1=2 independent of the dimension d.
In fact, the estimator that we propose has several appealing features.
1. It is relatively easy to compute (compared, say, to nearest-neighbour or kernel density estimators).
2. No tuning parameters need to be selected (unlike the "bandwidth" for kernel density estimators,
for example).
3. Well-established steady-state simulation output analysis techniques may be applied to analyze the
estimator.
4. The error in the estimator converges to 0 at rate n \Gamma1=2 independent of the dimension of the state
space, where n is the simulation runlength.
5. Under relatively mild assumptions, look-ahead density estimators consistently estimate not only
the density itself, but also the derivatives of the density; see Theorem 9.
6. The estimator can be used to obtain a new quantile estimator. The variance estimator for the
corresponding quantile estimator has a rigorous convergence theory, and converges at rate n \Gamma1=2
(Section 5).
7. Empirically, the estimator yields superior representations of stationary densities compared with
other methods (Example 1 of Section 6).
We first introduce the central ideas behind look-ahead density estimation in a familiar context. Although
this problem is subsumed by the treatment of Section 3, a separate development should prove
helpful in understanding the look-ahead approach. Let be an irreducible positive
recurrent Markov chain on finite state space S, and -(y) be the stationary probability of a point y 2 S.
Our goal is to estimate the stationary "density" -(\Delta); in the finite state space context, the stationary
"density" coincides with the stationary probabilities -(y), for y 2 S. To estimate -(y), the standard
estimator is
~
where I(\Delta) is the indicator function that is 1 if its argument is true, and 0 otherwise. The estimator ~ -n (y)
is simply the proportion of time the Markov chain X spends in the state y.
Notice however, that one could also estimate -(y) by
where P (\Delta; \Delta) is the transition matrix of X . The estimator -n (y) is a (strongly) consistent estimator of
-(y) as can be seen by noting that
as by the strong law for positive recurrent Markov chains on discrete state space. Notice that
so that the quantity P (X(i); y) is, in effect, "looking ahead" to
see whether the next iterate of the Markov chain will equal y. This is the motivation for the name
"look-ahead" density estimator.
In the remainder of this paper we assume a general state-space (not necessarily discrete) unless
otherwise specified. We refer to the density we are trying to estimate as the target density, and the
associated distribution as the target distribution.
In Section 2, look-ahead density estimators are developed for several performance measures associated
with transient simulations, and their pointwise asymptotic behaviour is derived. Steady-state performance
measures are similarly considered in Section 3. In Section 4, we turn to the global convergence behaviour
of look-ahead density estimators. In particular, we give conditions under which the look-ahead density
estimator converges to the target density in an L q sense (Theorem 5), is uniformly convergent (Theorem
7), and is differentiable (Theorem 9).
In Section 5 we consider the computation of several features of the target distribution, including the
mode of the target density, and quantiles of the target distribution. Finally, in Section 6 we give three
examples of look-ahead density estimation.
Computing Densities for Transient Performance Measures
be a Markov chain taking values in a state space S. Since our focus in this
section is on transient performance measures, we will permit our chain to possess transition probabilities
that are non-stationary.
Recall that is a transition kernel if Q(x; \Delta) is a probability measure on S
for each x 2 S, and if Q(\Delta; dy) is suitably measurable. (If S is a discrete state space, Q corresponds to
a transition matrix.) By permitting X to have non-stationary transition probabilities, we are asserting
that there exists a sequence of transition kernels such that
a.s.
for Our basic assumption is:
A1. There exists a (oe-finite) measure fl on S and a function
for
Remark 1: Assumption A1 is automatically satisfied when S is finite or countably infinite.
Remark 2: Given that this paper is concerned with density estimation, the case where fl is Lebesgue
measure and S is a subset of IR d is of the most interest to us. However, it is important to note that A1
does not restrict us to this context. In fact, Example 1 in Section 6 shows that this apparent subtlety
can in fact be very useful.
Remark 3: If X has stationary transition probabilities, P
In our discussion of steady-state density estimation (see Section 3), we will clearly wish to restrict
ourselves to such chains.
We will now describe several different computational settings to which the ideas of this paper apply.
In what follows, we will adopt the generic notation pZ (\Delta) to denote the fl-density of the r.v. Z. In other
words, pZ (\Delta) is a function with the property that
for all y in the range of Z. Also, for a given initial distribution - on S, let P - (\Delta) be the probability
distribution on the path-space of X under which X has initial distribution -.
Problem 1: Compute the density of X(r).
For r - 1, let p X(r) (\Delta) be the fl-density of X(r). Note that
Z
Z
so that
Z
is the expectation operator corresponding to P - . To compute the density p X(r) (y), simulate n
i.i.d. replicates of X under P - . Then, A1 and the strong law of large numbers together
guarantee that
a.s.
as n !1, so that p X(r) (y) can indeed be computed by our look-ahead estimator p 1n (y).
Remark 4: Suppose that A1 is weakened to
for so that now we are assuming the existence of a density only for the m-step
transition probability distribution. Provided r - m, we can write
so that p X(r) (y) can again be easily computed via independent replication of X .
For a given subset A ' S, let Ag be the first entrance time to A.
Problem 2: Compute the density of X(T ).
Suppose that P - starts in A c under initial distribution -. Then, for
Z
Z
so that for y 2 A,
Z
Again, A1 and the strong law of large numbers ensure that
as n !1, where the X i 's are independent replicates of X under P - , and T Ag.
An important class of transient performance measures is concerned with cumulative costs. Specifically,
be a sequence of real-valued r.v.'s in which \Gamma(n) may be interpreted as the ``cost''
associated with running X over
\Gamma(i)
is the cumulative cost corresponding to the time interval [0; n). We assume that
Y
so that, conditional on X , the \Gamma(i)'s are independent r.v.'s and the conditional distribution of \Gamma(i) depends
on X only through X(i \Gamma 1) and X(i). An important special case arises when
IR. In this case, (1) is automatically satisfied and f(x)
may be viewed as the cost associated with spending a unit amount of time in x 2 S. (We permit the
additional generality of (1) because such cost structures are a standard ingredient in the general theory
of "additive functionals" for Markov chains, and create no difficulties for our theory.)
Before proceeding to a discussion of the cumulative cost C(n), we note that Problems 1 and 2 have
natural analogues here. However, we will need to replace A1
A2. There exists a (oe-finite) measure fl on S and a function ~
for
Problem 3: Compute the density of \Gamma(r).
For y 2 IR, the density p \Gamma(r) (y) can be consistently estimated by
=n
~
are independent replicates of X .
Problem 4: Compute the density of \Gamma(T ).
Here, the density p \Gamma(T ) (y) can be consistently estimated via
=n
~
As usual, are independent replicates of X under P - , and T i is the first entrance time of
X i to the set A.
In addition to consistency of p 3n (y) and p 4n (y), A2 permits us to solve a couple of additional computational
problems that relate the to the density of the cumulative cost r.v. introduced earlier.
Problem 5: Compute the density of C(r).
We assume here that fl is Lebesgue measure. Then, if r - 1, we may use A2 to write
Z
Z
Z y\Gammaz
so that
Z
Z
Evidently, A2 and the strong law of large numbers together guarantee that
~
a.s.
as n !1, so that p 5n (y) is a consistent estimator of p C(r) (y), the (Lebesgue) density of C(r).
Problem Compute the density of C(T ).
As we did earlier, we assume that arguments to those
used above establish the identity
~
Thus A2 and the strong law prove that
is a consistent estimator for
To this point, we have constructed unbiased density estimators for each of the six density computation
problems described above. We now turn to the development of asymptotically valid confidence regions
for these densities. The key is to recognize that each of the six estimators may be represented as
p in
1. (Here, the index set   is either S or IR, depending on which of the
estimators is under consideration.) For y 2  , let p(y; i) 4
(y). For d points y
in (~y) (~p(~y; i)) to be a d-dimensional vector with jth component p in (y j
A straightforward application of the multivariate central limit theorem (CLT) yields the following result.
points in  , selected so that E- ij (y k
as is a d-dimensional multivariate normal random vector with mean vector
zero and covariance matrix \Sigma i (~y) having (j; k)th element given by cov(- i1 (y j ); - i1 (y k )).
Proposition 1 suggests the approximation
~p in (~y) D
for n large, where D
- denotes the (non-rigorous) relation "has approximately the same distribution as",
and \Sigma 1=2
i (~y) is a square root (Cholesky factor) of the non-negative definite symmetric matrix \Sigma i (~y). Since
easily estimated consistently from X by the sample covariance matrix, it follows that
(2) may be used to construct asymptotically valid confidence regions for ~ p(~y; i).
Remark 5: Equation (2) implies that the error in the look-ahead density estimator decreases at rate
n \Gamma1=2 . This dimension-independent rate stands in sharp contrast to the heavily dimension-dependent
rate exhibited by other density estimators, including kernel density estimators; see Prakasa Rao (1983).
The convergence rate for such estimators is typically (nh d
is the bandwidth parameter
and d is the dimension of  . To minimize mean squared error, the bandwidth hn is typically chosen to
be of the order n \Gamma1=(d+4) , and then the error in the kernel density estimators decreases at rate n \Gamma2=(d+4) .
Even in one dimension, this asymptotic rate is slower than that exhibited by the look-ahead density
estimator, and in higher dimensions, the difference is even more apparent; see Example 2 of Section 6.
Computing Densities for Steady-State Performance Measures
We now extend our look-ahead estimation methodology to the steady-state context. In order for the
concept of steady-state to be well-defined, we assume that X has stationary transition probabilities, so
that the transition kernels introduced in Section 2 are independent of n. In other words,
we assume that there exists a transition kernel P such that P
for defined as in A1.
Proposition 2 Under A1, any stationary distribution of X possesses a density - with respect to fl.
Proof: Let - be a stationary distribution of X . (Note that we are using - to represent both the
stationary distribution and its density with respect to fl. The appropriate interpretation should be clear
from the context.) Then,
for all (suitably) measurable
Z
Z
Z
Z
Z
It follows from (3) and (4) that the stationary distribution - has a fl-density -(\Delta) having value
Z
at y 2 S. 2
According to Proposition 2, the density -(y) may be expressed as an expectation, namely
see (5). Relation (6) suggests using the estimator
-n (y) =n
to compute -(y); -n (y) requires simulating X up to time
To establish laws of large numbers and CLT's for -n (y), we require that X be positive recurrent in a
suitable sense.
A3. Assume that there exists a subset B ' S, positive scalars -; a and b, an integer m - 1, a probability
distribution '(\Delta) on S, and a (deterministic) function
1.
2. E[V
where I(x 2 B) is 1 or 0 depending on whether or not x 2 B.
In the language of general state-space Markov chain theory, A3 ensures that X is a geometrically
ergodic Harris recurrent Markov chain; see Meyn and Tweedie (1993) for details. Condition 1 of A3 is
typically satisfied for reasonably behaved Markov chains by choosing B to be a compact set; -; ', and m
are then determined so that 1 is satisfied. Condition 2 is known as a Lyapunov function condition. For
many chains, a potential choice for V is something of the form V
great deal of ingenuity may be necessary in order to construct such a V . See Example 1 of Section 6 for
an illustration of the verification of A3. In any case, A3 ensures that X possesses a unique stationary
distribution.
Remark Assumption A3 is a stronger condition than is necessary to obtain the laws of large numbers
and CLT's below. However, in most applications, A3 is a particularly straightforward sufficient condition
to verify, and we offer it in that spirit.
consist of d points selected from S, and let ~- n (~y) (~-(~y)) be a d-dimensional vector
in which the jth component is -n (y j
Theorem 3 Assume A1, A3, and suppose that for 1
as n !1. Also, there exists a non-negative definite symmetric matrix
as is a d-dimensional standard Brownian motion and ) denotes weak
convergence in D[0; 1).
Proof: The proof follows directly from results from Meyn and Tweedie (1993). The strong law is a
consequence of Theorem 17.0.1. Lemma 17.5.1 and Lemma 17.5.2 together imply the existence of a
square integrable (with respect to -) solution to Poisson's equation. This then enables an application of
Theorem 17.4.4 to yield the result. 2
Remark 7: Equation (2) implies that the error in the look-ahead density estimator for estimating
stationary densities decreases at rate n \Gamma1=2 . This is the same rate we observed in Remark 5 for the
case of independent observations. Furthermore, exactly as in the independent setting, other existing
estimators, including kernel density estimators, converge at a slower rate; see Yakowitz (1989). The
convergence rate for such estimators is typically (nh d
is the bandwidth parameter, and d
is the dimension of the (Euclidian) state space. Since hn ! 0 as convergence rate is slower
than n \Gamma1=2 .
Yakowitz (1989) does not give the optimal (in terms of minimizing mean squared error) choice of
bandwidth hn . However, an i.i.d. sequence is a special case of a Markov chain, and, as noted in Remark
5, the fastest possible root mean square error convergence rate in that setting is of the order n \Gamma2=d+4 .
This rate is heavily dimension dependent, so that in large-dimensional problems, one might expect very
slow convergence of kernel density estimators.
To obtain confidence regions for the density vector ~-(~y), several different approaches are possible. If
\Sigma(~y) is positive definite and there exists a consistent estimator \Sigma n (~y) for \Sigma(~y), then Theorem 3 asserts
that for D ' IR d ,
for large n, where \Sigma n (~y) 1=2 D is defined to be the set fx
w for some ~
confidence regions for ~-(~y) can then be easily obtained from (7). If X enjoys regenerative structure,
the regenerative method for steady-state simulation output analysis provides one means of constructing
such consistent estimators for \Sigma(~y); see, for example, Bratley, Fox, and Schrage (1987).
An alternative approach exploits the functional CLT provided by Theorem 3 to ensure the asymptotic
validity of the method of multivariate batch means; see Mu~noz and Glynn (1998) for details.
Remark 8: The discussion of this section generalizes to the computation of the density \Gamma(1) under the
stationary distribution -. In particular, suppose that X satisfies A2 and A3. Then for y 2 IR,
where ~
defined as in A2; the methodology
of this section then generalizes suitably.
Global Behaviour of the Look-ahead Density Estimator
In the previous section we focused on the pointwise convergence properties of the look-ahead density
estimator. Specifically, we showed that for any finite collection y of points in either S or IR
(depending on the estimator), the look-ahead density estimator converges a.s. and the rate of convergence
is described by a CLT in which the rate is dimension independent. In this section, we turn to the
estimator's global convergence properties. We assume throughout the remainder of this paper that S is
a complete separable metric space. (In particular, this includes any state space that is a "reasonable"
subset of IR k .)
Let fl be as in Sections 2 and 3, and let   be either S or IR (depending on the estimator considered).
Then, for any function f :   ! IR, we may define, for q - 1, the L q -norm
kfk q
'Z

For any two functions f 1 and f is a measure of the "distance" from f 1 to f 2 . We first analyze
the look-ahead density estimators introduced in Section 2.
Theorem 4 Suppose that E -
R

kp in (\Delta) \Gamma p(\Delta; i)k q
as n !1.
Proof: Evidently,
for y. Note that
-n
due to the convexity of jxj q . For each y satisfying (8), the right hand side of (9) converges a.s. and in
expectation to E - j- i1 (y) \Gamma p(y; i)j q . Consequently, the right-hand side of (9) is uniformly integrable.
Also, for each such y, the left-hand side of (9) converges to zero a.s. Since the left-hand side is dominated
by a uniformly integrable sequence, it follows that
Ejp in
as n !1 for fl a.e. y. Also, taking the expectation of both sides of (9) yields the inequality
Ejp in
the right-hand side is integrable in y, by hypothesis. The Dominated Convergence Theorem, applied to
(10), then gives Z

Ejp in
and hence
Ekp in (\Delta) \Gamma p(\Delta; i)k q
Consequently,
kp in (\Delta) \Gamma p(\Delta; i)k q
as n !1, from which the theorem follows. 2
We turn next to obtaining the analogous result for the steady-state density estimator -n (\Delta) of Section
3.
Theorem 5 Suppose Z
for x 2 S, with q - 1. If A3 holds and the initial distribution - has a density with respect to fl, then
as n !1.
Proof: Condition (11) guarantees that
Z
see Theorem 14.3.7 of Meyn and Tweedie (1993). So, y, and the proof
follows the same pattern as that for Theorem 2. That argument yields the conclusion that for ffl ? 0,
as almost every x. Therefore, the
Dominated Convergence Theorem allows us to conclude that
Z
as n !1, where u(\Delta) is the -density of -. This is the desired conclusion. 2
Remark 9: Note that the hypotheses of both Theorems 2 and 3 are automatically satisfied when
Convergence of the estimated density in L q ensures that for a given runlength n, errors of a given size
can only occur in a small (with respect to fl) set.
We now turn to the question of when the look-ahead density estimator converges to its limit uniformly.
Uniform convergence is especially important in a visualization context. If one can guarantee that the
error in the estimator is uniformly small, then graphs of the estimated density will be "close" to the
graph of the limit.
We will focus our attention here on the steady-state density estimator -n ; similar results can be
derived for our other density estimators through analogous arguments.
Theorem 6 Suppose that A1 is in force, and p : S \Theta S ! [0; 1) is continuous and bounded. If A3
holds, then for each compact set K,
sup
as n !1.
Proof: Fix ffl ? 0. Since - is tight (see Billingsley 1968), there exists a compact set K(ffl) for which -
assigns at most ffl mass to its complement. Write
1. The second term on the right-hand side of (12) may be bounded by
2 K(ffl)), which has an a.s. limit supremum of at most -ffl. As for the first term, note
that if K is compact, then K(ffl) \Theta K is compact and p is therefore uniformly continuous there. Because
of uniform continuity, there exists ffi(ffl) such that whenever \Theta K is within distance ffi(ffl)
of ffl. Since K is compact, we can find a finite collection
of points in K such that the open balls of radius ffi(ffl) centered at y
each y 2 K, there exists y i in our collection such that jp(X(j);
So, for y 2 K,
Letting 1)ffl. Hence, for y 2 K, we obtain the uniform
bound
By letting n ! 1, applying the strong law for Harris chains to n \Gamma1
sending ffl ! 0, we obtain the desired conclusion. 2
Remark 10: If S is compact, Theorem 6 yields uniform convergence of -n to - over S, under a
continuity hypothesis on p. (The boundedness is automatic in this setting.)
Our next result establishes uniform convergence of -n to - over all of S, provided that we assume
that p(x; \Delta) "vanishes at infinity".
Theorem 7 Suppose that A1 holds and uniformly continuous and bounded.
Assume that for each x 2 S and ffl ? 0, there exists a compact set K(x; ffl) such that whenever
ffl. If A3 holds, then
sup
as n !1.
Proof: Fix ffl ? 0, and choose ffi(ffl) so that whenever lies within distance ffi(ffl) of
choose K(ffl) as in the proof of Theorem 6 and let x K(ffl) be a finite
collection of points such that the open balls of radius ffi(ffl) centred at x K(ffl). For each x i ,
there exists K i
and note that K is compact. Theorem 6 establishes that
sup
as n !1. To deal with
construct the sequence (X
2 K(ffl) and X 0 (n) is the closest point within the collection fx
Then, for
Sending n !1 allows us to conclude that -(y) -
K. The inequality (14) then yields
lim sup
sup
together imply the theorem. 2
The following consequence of Theorem 7 improves Theorem 5 from convergence in probability to a.s.
convergence when basically Scheff'e's Theorem (see, for example, p. 17, Serfling 1980).
Corollary 8 Under the conditions of Theorem 7,
Z
as n !1.
Proof: The result is immediate if fl is a finite measure (since j- n (\Delta) \Gamma -(\Delta)j is uniformly bounded and
converges to zero a.s. by Theorem 7. If fl is an infinite measure (like Lebesgue measure), Theorem 7
asserts that -n (\Delta) ! -(\Delta) a.s. so that path-by-path, we may argue that
Z
Z
(since the integrand is dominated by -(\Delta), which integrates to one, thereby permitting the application of
the Dominated Convergence Theorem path-by-path). 2
A very important characteristic of the look-ahead density estimator is that it "smoothly approximates"
the density to be computed. To be specific, suppose that either or that we are considering the
density of one of the real-valued r.v.'s associated with the estimators p 3n (\Delta); p 4n (\Delta), p 5n (\Delta) or p 6n (\Delta). Since
we are then working in a subset of Euclidian space, it is reasonable to measure smoothness in terms of
the derivatives of the density.
any real loss of generality, assume so that y 2 IR. The look-ahead density estimators
we have developed take the form
for some sequence of random functions (G 1). (Both the estimators of Section 2 and Section
3 admit this representation.) To estimate the kth derivative of the target density to be computed, the
natural estimator is therefore
dy k pn
dy
Under quite weak conditions on the problem, it can be shown that the above estimator computes the
kth derivative of the target density consistently; see below for a discussion. Such a result proves that not
only does look-ahead density estimation compute the density, but it also approximates the derivatives of
the density in a consistent fashion. In other words, it "smoothly approximates" the target density.
As an illustration of the types of conditions needed in order to ensure that the look-ahead density
estimator smoothly approximates the target density, we consider the steady-state density estimator of
Section 3. Let p 0 (x;
dy p(x; y).
Theorem 9 Suppose A1 holds, continuously differentiable with bounded
derivative. If A3 holds, then - has a differentiable density -(\Delta), and
sup
as n !1 for each compact K ' S. Furthermore, if jp 0 (x; y)j - V 1=2 (x) for x 2 S, then there exists d(y)
such that
as n !1.
Proof: Note that
Z
lies between y and y Because the derivative is
assumed to be bounded, the Bounded Convergence Theorem then ensures that the limit in (18) exists
and bounded and continuous, exactly the same argument as that used
in proving Theorem 6 can be used here to obtain (16).
The CLT (17) is an immediate consequence of Theorems 17.0.1 and 17.5.4 of Meyn and Tweedie
(1993). 2
An important implication of Theorem 9 is that the look-ahead density estimator computes the derivative
accurately. In fact, the density estimator converges at rate n \Gamma1=2 , independent of the dimension of
the state space, and furthermore, independent of the order of the derivative being estimated.
Remark 11: It is also known that kernel density estimators smoothly approximate the target density;
see Prakasa Rao 1983 p. 237, and Scott 1992 p. 131. The choice of bandwidth that minimizes mean
squared error of the kernel density derivative estimator is larger than in the case of estimating the target
density itself. The resulting rates of convergence of kernel density derivative estimators are adversely
affected by both the order of the derivatives, and the dimension of the state space. For example, in
one dimension, kernel density derivative estimators of an rth order derivative converge at best at rate
This rate is fastest when estimating the first derivative, and even then,
is slower than the rate of convergence of the look-ahead density derivative estimator discussed above.
Computing Special Features of the Target Distribution Using
Look-ahead Density Estimators
As discussed earlier, computation of the density is a useful element in developing visualization tools
for computer simulation. In this section, we focus on the computation of certain features of the target
distribution, to which our look-ahead density estimator can be applied to advantage.
5.1 Computing the Relative Likelihood of Two Points
In Sections 2 and 3, we introduced a number of different look-ahead density estimators, each of which we
can write generically as pn (\Delta). The look-ahead density estimator pn (\Delta) is an estimator for a target density
p(\Delta) say. For each pair of points (y represents the likelihood of the point y 1
relative to that of y 2 .
The joint CLT's developed in Proposition 1 and Theorem 3 can be used to obtain a CLT (suit-
able for construction of large-sample confidence intervals for the relative likelihood) for the estimator
as
as n !1. If the covariance matrix of (N 1 can be consistently estimated (as with, for example, the
regenerative method), then confidence intervals for the relative likelihood (based on (19)) can easily be
obtained. Otherwise, one can turn to the batch means method to produce such confidence intervals; see
Mu~noz and Glynn (1997).
5.2 Computing the Mode of the Density
The mode of the density provides information as to the region within which the random variable of
interest attains its highest likelihood. Given that the target distribution here has density p(\Delta), our goal
is to compute the modal location y   and the modal value p(y   ). As discussed earlier in this section, we
our look-ahead density estimator generically as pn (\Delta). The obvious estimator of y   is, of course, any
y
which maximizes pn (\Delta), and the natural estimator for p(y   ) is then pn (y
can and will assume
that the maximizer y
n has been selected to be measurable.) We denote the domain of p(\Delta) by  . Because
our analysis involves using a Taylor expansion, we require that   2 IR d .
Theorem
1. p(\Delta) has a unique mode at location y   ;
2. sup
a.s. as n !1;
3. there exists an ffl-neighbourhood of y   with ffl ? 0 such that p(\Delta) and pn (\Delta) are twice continuously
differentiable there a.s.;
4. sup
ky\Gammay   k!ffl
a.s. as n !1;
5. sup
ky\Gammay   k!ffl
a.s. as n !1, where Hn (y) and H(y) are the Hessians of pn (\Delta) and
p(\Delta) at y, respectively;
6. H(y   ) is negative definite;
7. n 1=2 (p n (y
Then, y
a.s. as n !1 and
as n !1.
Proof: The almost sure convergence of y
n to y   is an immediate consequence of relations 1 and 2.
For the weak convergence statement, observe that rpn (y
are local
maxima of pn (\Delta) and p(\Delta), respectively. (To be precise, this is valid only for n so large that y
n lies in the
ffl-neighbourhood specified by relations 3 - 5.) So,
rpn (y
But
rpn (y
a.s. as n !1. It follows that n 1=2 (y
\Gammay
lies on the line segment joining y
n and y   . Since rpn (- n a.s. and we have
established weak convergence of n 1=2 (y
evidently the first term in (21) converges to zero in
probability. Consequently, (20), (21), relation 7, and the "converging together principle" (see Billingsley
1968, for example) imply the desired joint convergence result. 2
Remark 12: The uniform convergence theory of Section 3 for pn and its derivatives can be easily
applied to verify relations 2, 4 and 5.
Remark 13: The CLT established in Theorem 10 shows that the look-ahead estimator of the mode
converges at the asymptotic rate n \Gamma1=2 , independent of the dimension d of the state space. This compares
very favourably with the rate of convergence of kernel estimators of the mode. A kernel estimator of the
mode converges at rate (nh d+2
when the bandwidth hn is chosen appropriately; see Theorem 4.5.6
of Prakasa Rao 1983 p. 284.
Remark 14: To construct confidence intervals based on Theorem 10, there are again a couple of
alternatives. Assume first that for each fixed y 2  , one can consistently estimate the covariance matrix
that arises in the joint CLT of relation 7. (For example, this can be done in the transient context or
the setting of regenerative processes in steady-state simulation). To estimate the covariance structure
at y   , one can compute the corresponding covariance estimate evaluated at the point
n . In the
transient problems considered in Section 2, it is typically easy to verify that that the covariance matrix
is continuous in y, so that using the estimator associated with y
n in place of the covariance at y   is
asymptotically valid. In the steady-state context, it is not as straightforward to theoretically establish
the continuity of the covariance (although one suspects it is valid in great generality); one potential avenue
is to adopt the methods of Glynn and L'Ecuyer (1995). If consistent estimates of the covariance matrix
at each fixed y 2   are not available (as might occur in non-regenerative steady-state simulations), then
one can potentially appeal to the method of batch means; see Mu~noz (1998).
5.3 Computing Quantiles of the Target Distribution
We focus here on the special case in which   ' IR, so that the target distribution is that of a real-valued
r.v. In this setting, suppose that dy, and let
p(y) dy
be the target distribution. An important special feature of this distribution is the pth quantile of F .
Specifically, for each p 2 (0; 1), we define the pth quantile of F as the quantity
There is a significant literature on the computation of such quantiles. Iglehart (1976) considered
quantile estimation in the context of regenerative simulation, and proved a central limit theorem for the
standard estimator. Seila (1982) introduced the batch quantile method, again for regenerative processes,
that avoids some of the difficulties associated with the estimation procedure proposed by Iglehart (1976).
The approach suggested by Heidelberger and Lewis (1984) is based on the so-called "maximum transfor-
mation" and mixing assumptions of the underlying process, and does not require regenerative structure.
Hesterberg and Nelson (1998) and the references therein discuss the use of control variates to obtain
variance reduction in quantile estimation. Avramidis and Wilson (1998) obtain variance reduction in
estimating quantiles through the use of antithetic variates and Latin hypercube sampling.
Kappenman (1987) integrated and inverted a kernel density estimator for p(\Delta) in the case when the
observations are i.i.d. Our approach is similar to Kappenman's in that we invert the integrated look-ahead
density estimator. Let pn (\Delta) be the look-ahead density estimator, and set
pn (y) dy:
The natural estimator for the quantile q is then
(p).
Theorem 11 Suppose that
1. p(q) ? 0;
2. p(\Delta) is continuous in an ffl neighbourhood of q;
3. sup
a.s. as n !1;
4. n 1=2
Proof: Recall that kpn (\Delta) \Gamma p(\Delta)k 1 ) 0 as n !1; see Remark 9. Hence,
sup
x
as n !1, from which it follows that Qn ) q as n !1. But
and
lies between Qn and q. Relations 2 and 3 imply that pn (- n ) ) p(q) as n !1. The result then
follows from (22), (23), and relation 4. 2
Remark 15: Similar issues to those discussed in Remark 14 arise in constructing confidence intervals
based on Theorem 11. Once again, it is possible to consistently estimate the variance parameter that
arises in the CLT in relation 4 in either the transient context, or the setting of regenerative steady-state
simulation. To see this, recall that the look-ahead density estimator pn (\Delta) may be expressed as
some sequence of random functions (G
pn (y) dy
=n
Evidently, (24) is a sample mean over a sequence of real-valued r.v.'s, and the sequence is either i.i.d.
or regenerative, depending on the context. Therefore, standard methods may be applied to estimate the
variance parameter in the CLT of relation 4.
The comments in Remark 14 related to continuity of the variance parameter apply directly here. In
particular, one must establish that the variance of the r.v. N in relation 4 is continuous as a function of
q, so that estimating the variance at q by an estimate of the variance at q n is asymptotically valid.
It is natural to ask how the performance of the look-ahead quantile estimator compares with that of
a more standard quantile estimator. For ease of exposition, in the remainder of this section we specialise
to the case where is a Markov chain taking values in IR. Suppose that A1 and A3 are
in force with and we are interested in computing is the distribution
function of the stationary distribution - of X .
A natural approach to estimation of q is to first estimate F by the empirical distribution function ~
Fn ,
where
~
and then choose the estimator ~
Qn of q as ~
(p).
Alternatively, using look-ahead methodology, one could estimate q by F \Gamma1
The proof of the following proposition rests primarily on the observation that
so that the estimators Fn and ~
Fn are related through the principle of extended conditional Monte Carlo
(Bratley et al. 1987 p. 71, Glasserman 1993).
Let p be a density of the stationary distribution - with respect to Lebesgue measure, and let var -
denote the variance operator associated with the path space of X , where X has initial distribution -.
Proposition 12 Suppose that A1 and A3 hold, and conditions 1 and 2 of Theorem 11 are satisfied.
Then,
as
~
Fn (q), and N 1 (0; 1) and N 2 (0; 1) are
standard normal r.v.'s. In addition, if X is stochastically monotone, then oe 2 - ~
oe 2 .
Proof: Observe that for all w 2 IR,
Z q
so that Theorem 17.5.3 of Meyn and Tweedie implies that n 1=2 Similarly,
Theorem 17.5.3 also gives n 1=2 ( ~
If X is stochastically monotone, then in view of (25), we can apply Theorem 12 of Glynn and
Iglehart (1988) to achieve the result. (The required uniform integrability follows from the fact that
~
Combining the results of Theorem 11 and Proposition 12, we see that under reasonable conditions,
as n !1. It can also be shown, again under reasonable conditions, that
~
oe
as Henderson and Glynn (1999). Proposition 12 asserts that oe 2 - ~
so that in the
context of steady-state quantile estimation for stochastically monotone Markov chains, the look-ahead
quantile estimator may typically be expected to achieve variance reduction over a more standard quantile
estimator.
Remark 16: It is well-known that the waiting time sequence in the single-server queue is a stochastically
monotone Markov chain, and thus the results of this section may be applied in that context.
6 Examples
We present three examples of the application of look-ahead density estimators. Our first example is an
example of steady-state density estimation and illustrates how to establish A3.
Example 1: It is well-known that the sequence of customer waiting times
(excluding service) in the FIFO single-server queue is a Markov chain on state space In
particular, W satisfies the Lindley recursion (p. 181, Asmussen 1987)
is an i.i.d. sequence with Y (n
V (n) is the service time of the nth customer, and U(n + 1) is the interarrival time between the nth and
1)st customer.
To verify A3, we proceed as follows. Define some yet to be determined constant
note that
Let us assume that the moment generating function OE(t) 4 =Ee tY (1) of Y (1) exists in a neighbourhood
of zero, so that OE(t) is finite for sufficiently small t. For stability, we must have EY (1) ! 0, which
implies that OE 0 (0) ! 0. Hence, there exists an ff ? 0 such that OE(ff) ! 1. Now choose K ? 0 so that
then for all x ? K, we see from (26) that E[V (W (1))jW
where From (26) we also see that for x - K, E[V (W (1))jW
1. Thus, we have verified condition 2 of A3 for the set
To verify condition 1, note that EY (1) ! 0 implies that there exists
fi. It follows that conditional on W
Taking ' to be a point mass at 0, and we see that condition 1 of A3 is verified. We have
therefore established the following result.
Proposition 13 If the moment generating function of Y (1) exists in a neighbourhood of 0, and EY (1) !
0, then A3 is satisfied.
We now specialise to the M/M/1 queue with arrival rate -, service rate -, and traffic intensity
ae 4
1. The transition kernel for W is then given by
where p(x;
is the probability measure that assigns unit mass
to the origin. Noting that p(\Delta; \Delta) is bounded by max(-; 1), it follows (after possibly scaling the function
that the conditions of Theorem 3 are satisfied, and the look-ahead density estimator therefore
converges at rate n \Gamma1=2 to the stationary density of W .
Defining a suitable kernel density estimator is slightly more problematical, due to the presence of the
point mass at 0 in the stationary distribution and the need to select a kernel and bandwidth. To estimate
the point mass at 0 we use
the mean number of visits to 0 in a run of length n. For y ? 0, we estimate -(y) using
where
is the density of a standard normal r.v., and This choice of hn (modulo a multiplicative
constant) yields the optimal rate of mean-square convergence in the case where the observations are i.i.d.
(Prakasa Rao 1983, p. 182), and so it seems a reasonable choice in this context.
For this example we chose so that the traffic intensity ae = 0:5. To remove the
effect of initialization bias (note that both estimators are affected by this), we simulated a stationary
version of W by sampling W 0 from the stationary distribution.
exact
look-ahead
kernel
100.050.150.250.350.45waiting time
density
Density Estimators for the M/M/1 Queue

Figure

1: Density estimates from a run of length 100.
The density estimates for x ? 0, together with the exact density, are plotted for simulation runlengths
of

Figure

We observe the following.
1. Visually, the look-ahead density estimate appears to be a far better representation of the true
density than the kernel density estimate.
2. The kernel density estimate has several local modes, and its performance near the origin is particularly
poor, even for the run of length 1000.
The previous example is a one-dimensional density estimation problem. Our results suggest that the
rate of convergence of the look-ahead density estimators is insensitive to the underlying dimension of the
problem. However, the rate of convergence of kernel density estimators is known to be adversely affected
by the dimension; see Remarks 5 and 7. To assess the difference in performance in a multi-dimensional
setting, we provide the following example.
exact look-ahead
kernel
100.10.30.5waiting time
density
Density Estimators for the M/M/1 Queue

Figure

2: Density estimates from a run of length 1000.
Example 2: Let be a sequence of d dimensional i.i.d. normal random vectors with
zero mean and covariance matrix the identity matrix I . Define the Markov chain
inductively by
1. The Markov chain X is a (very) special case of the linear state space model defined
on p. 9 of Meyn and Tweedie (1993). We chose such a model for this example so that the steady-state
density is easily computed. In particular, the stationary distribution of X is normal with mean zero and
covariance matrix thus X has stationary density
exp
We estimate this density at dimensions using both a kernel density estimator
and a look-ahead density estimator, with estimators are constructed from simulated
sample paths of length 10, 100 and 1000. We sample X(0) from the stationary distribution to remove
any initialization bias. To estimate the mean squared error (MSE) of the density estimators at
repeat the simulations 100 times.
The kernel density estimator we chose uses a multivariate standard normal distribution as the kernel,
and a bandwidth for the rationale behind this choice of bandwidth).

Table

1 reports the root MSE for the two estimators as a percentage of the true density value -(0).
Observe that as the dimension increases, the rate of convergence of the kernel density estimator deteriorates
rapidly. In contrast, the rate of convergence of the look-ahead density estimator remains constant
(for each increase in runlength by a factor of 10, relative error decreases by a factor of approximately 3),
independent of the dimension of the problem.
Remark 17: It is possible to construct look-ahead density estimators for far more complicated linear
state space models than the one considered here. The critical ingredient is A1 which is easily satisfied,
d -(0) Estimator Runlength

Table

1: Root MSE of estimators of -(0) as a percentage of -(0).
for example, if the innovation vectors W (k) have a known density with respect to Lebesgue measure.
Our final example is an application to stochastic activity networks (SANs). This example is not easily
captured within our Markov chain framework, and therefore gives some idea of the potential applicability
of look-ahead density estimation methodology.
Example 3: In this example, we estimate the density of the network completion time (the length of the
longest path from the source to the sink) in a simple stochastic activity network taken from Avramidis
and Wilson (1998). Consider the SAN in Figure 3 with independent task durations, source node 1, and
sink node 9. The labels on the arcs give the mean task durations. We assume that tasks (6,
have densities (with respect to Lebesgue measure), so that the network completion time L has a density
p(\Delta) (with respect to Lebesgue measure).

Figure

3: Stochastic activity network with mean task duration shown beside each task.
Suppose that we sample all task durations except task (6, and compute the lengths
and L(8) of the longest paths from the source node to nodes 6 and 8 respectively. Then
F
where, for a given task ab, F ab denotes the task duration distribution function, -
F ab ab (\Delta), and
f ab (\Delta) is the (Lebesgue) density. Then, A1 and the strong law of large numbers ensure that the look-ahead
density estimator
F
is a strongly consistent estimator of p(y).
For the purposes of our simulation experiment, we assumed that all task durations were exponentially
distributed with means as indicated on Figure 3. The resulting density estimate is depicted in Figure 4
for a run of length 1000.
1500.0050.0150.025Completion Time
density
Estimated Network Completion Time Density

Figure

4: Estimate of the Network Completion Time Density.
Remark 18: The approach taken in this example clearly generalizes to other SANs where all of the
arcs entering the sink node have densities (with respect to Lebesgue measure).
Remark 19: One need not base a look-ahead density estimator on the arcs that are incident on the
sink. For example, one might instead focus on arcs that leave the source. In the above example, these
arcs correspond to tasks (1, 2) and (1, 3), and one would condition on the longest paths from nodes 2
and 3 to the sink.

Acknowledgments



--R

Applied Probability and Queues.

Convergence of Probability Measures.
A Guide to Simulation
Nonparametric Density Estimation: The L 1 View.
A Course in Density Estimation.
Filtered Monte Carlo.
Estimation of Stationary Densities of Markov Chains.
Likelihood ratio gradient estimation for stochastic recursions.
Quantile estimation in dependent sequences.
Regenerative steady-state simulation of discrete-event systems
Asymptotic results for steady-state quantile estimation in Markov chains
Control variates for probability and quantile estimation.
Simulating stable stochastic systems
Improved distribution quantile estimation.
Markov Chains and Stochastic Stability.
A batch means methodology for the estimation of quantiles of the steady-state distribution
Multivariate standarized time series for output analysis in simulation experiments.
Batch means methodology for estimation of a nonlinear function of a steady-state mean
Nonparametric Functional Estimation.
Multivariate Density Estimation: Theory
A batching approach to quantile estimation in regenerative simulations.
Approximation Theorems of Mathematical Statistics.
Nonparametric density estimation
Nonparametric density and regression estimation for Markov sequences without mixing assumptions.
--TR

--CTR
Shane G. Henderson, Simulation mathematics and random number generation: mathematics for simulation, Proceedings of the 33nd conference on Winter simulation, December 09-12, 2001, Arlington, Virginia
density estimator;markov chain;simulation