Using permutations in regenerative simulations to reduce variance.
We propose a new estimator for a large class of performance measures obtained from a regenerative simulation of a system having two distinct sequences of regeneration times.  To construct our new estimator, we first generate a sample path of a fixed number of cycles based on one sequence of regeneration times, divide the path into segments based on the second sequence of regeneration times, permute the segments, and calculate the performance on the new path using the first sequence of regeneration times.  We average over all possible permutations to construct the new estimator.  This strictly reduces variance when the original estimator is not simply an additive functional of the sample path.  To use the new estimator in practice, the extra computational effort is not large since all permutations do not actually have to be computed as we derive explicit formulas for our new estimators.  We examine the small-sample behavior of our estimators.  In particular, we prove that for any fixed number of cycles from the first regenerative sequence, our new estimator has smaller mean squared error than the standard estimator.  We show explicitly that our method can be used to derive new estimators for the expected cumulative reward until a certain set of states is hit and the time-average variance parameter of a regenerative simulation.
INTRODUCTION
The regenerative method is a simulation-output-analysis technique for estimating
certain performance measures of regenerative stochastic systems; see [Crane and
Iglehart 1975]. The basis of the approach is to divide the sample path into i.i.d.
segments (cycles), where the endpoints of the segments are determined by a sequence
of stopping times. Many stochastic systems have been shown to be regenerative
[Shedler 1993], and the regenerative method results in asymptotically valid
confidence intervals.
In this paper we propose a new simulation estimator for a performance measure
of a regenerative process having two different sequences of regeneration times, and
study its small-sample behavior. The idea of our approach is as follows. First simulate
a fixed number of regenerative cycles from the first sequence of regeneration
times, and compute one estimate. We construct another estimator by dividing up
the original sample path into segments with endpoints given by the second sequence
of regeneration times, and creating a new sample path by permuting the segments
(except for the initial and final segments). We then compute a second estimate of ff
from the new permuted path. We show that this estimate has the same distribution
as the original one. Our new estimator is finally constructed as the average of the
estimates over all possible permutations. This strictly reduces variance when the
estimator is not a purely additive function of the sample path. We show that to
compute our new estimators, one does not have to actually calculate all permutations
and the average over all of them. Instead, we derive formulas for the new
estimators, where the expressions can be easily computed by accumulating some
extra quantities during the simulation. The storage requirements of our methods
are fixed and do not grow as the simulation run length increases. Hence, there is
little extra computational effort or storage needed to construct our new estimators.
For a run length of any fixed number of cycles from the first regenerative se-
quences, the new estimator has the same expected value as the standard estimator
and lower variance; thus, it has lower mean squared error. While it turns out that
our method has no effect on the standard regenerative ratio estimator for certain
steady-state performance measures, the basic technique can still be beneficially applied
to a rich class of other performance measures, and in this paper, we consider
three specific examples.
First, we derive a new estimator for the second moment of the cumulative reward
during a regenerative cycle. We show that the standard regenerative variance
estimator fits into this framework. Hence, our estimator will result in a variance
estimator having no more variability than the standard one. This is important
because one measure of the quality of an output-analysis methodology is the variability
of the half-width of the resulting confidence interval [Glynn and Iglehart
1987], which is largely influenced by the variance of the variance estimator.
We also construct a new estimator for the cumulative reward until some set of
states is first hit, which includes the mean time to failure as a special case. Here,
the performance measure can be expressed as a ratio of expectations, and we apply
Using Permutations in Regenerative Simulations to Reduce Variance \Delta 3
our technique to the numerator and denominator separately.
In some sense our method reuses the collected data to construct a new estimator,
and as such, it is related to other statistical techniques. For example, the bootstrap
[Efron 1979] takes a given sample and resamples the data with replacement. In
contrast, one can think of our approach as resampling the data without replacement
(i.e., permuting the data), and then averaging over all possible resamples. Other
related methods include U-statistics (Chapter 5 of [Serfling 1980]), V -statistics [Sen
1977], and permutation tests (e.g., [Conover 1980]).
The rest of the paper is organized as follows. In Section 2, we discuss our assumptions
and the standard estimator of a generic performance measure ff. We present
the basic idea of how to construct our new estimator using a simple example in
Section 3. Section 4 contains a more formal description of our method. Section 5
describes the new estimator for the second moment of the cumulative reward over a
regenerative cycle and shows how these results can be used to derive a new estimator
of the variance parameter arising in a regenerative simulation. We also discuss
here the special case of continuous-time Markov chains. In Section 6 we derive
new estimators for the expected cumulative reward until some set of states is hit.
We analyze the storage and computational costs of our new estimator in Section 7.
We present in Section 8 the results of some simulation experiments comparing our
new estimators with the standard ones. Section 9 discusses directions for future re-
search. Most of the proofs are collected in Appendix A. Also, we give pseudo-code
for one of our estimators in Appendix B. (Calvin and Nakayama [1997] present the
basic ideas of our approach, without proofs, in the setting of discrete-time Markov
chains.)
2. GENERAL FRAMEWORK
Let X be a continuous-time stochastic process having sample paths that are right
continuous with left limits on a state space S ae ! d . Note that we can handle
discrete-time :g in this framework by letting
X btc for all t - 0, where bac is the greatest integer less than or equal to a.
be an increasing sequence of nonnegative finite
stopping times. Consider the random pair (X; T ) and the shift
We define the pair (X; T ) to be a regenerative process (in the classic sense) if
(i). f' T (i) (X; are identically distributed;
(ii). for any i - 0, ' T (i) (X; T ) does not depend on the "prehistory"
See p. 19 of [Kalashnikov 1994] for more details. This definition allows for so-called
delayed regenerative processes (e.g., Section 2.6 of [Kingman 1972]).
be two distinct increasing sequences of nonnegative finite stopping times such that
are both regenerative processes. For example, if X is an ir-
reducible, positive-recurrent, discrete-time or continuous-time Markov chain on a
countable state space S, then we can define T 1 and T 2 to be the sequences of hitting
M. Calvin and M. K. Nakayama
times to the states v 2 S and w 2 S, respectively, where we assume that
and w 6= v.
Our goal is to estimate some performance measure ff, which we will do by generating
a sample path segment ~
of a fixed number m 1
of regenerative 1-cycles of our regenerative process. Here, we use the terminology
"1-cycles" to denote cycles determined by the sequence T 1 ; i.e., the ith 1-cycle is
the path segment (i)g. We similarly define "2-cycles"
relative to the sequence T 2 . Now we define the standard estimator of ff based on
the sample path ~
1-cycles to be
where h j hm1 is some function. This general framework includes many performance
measures of interest.
Example 1. Suppose
"/ Z
for some function g : S ! !, where p - 1. Then we can define h( ~
X) by
h( ~
Y (g;
where
Z
for k - 1. Note that here b
X) is an unbiased estimator of ff. We will examine
this example with in Section 5.
Example 2. Suppose that
where for k - 1,
now is the variance parameter arising from a regenerative simulation. (More details
are given in Section 5.1.) Then we can define h( ~
X) by
h( ~
where
Using Permutations in Regenerative Simulations to Reduce Variance \Delta 5
Note that b
X) is the standard regenerative estimator of oe 2 . We will return to this
example in Section 5.1.
Example 3. Suppose we are interested in computing
set of states F ae S. Thus, j is the
expected cumulative reward until hitting F conditional on T 1 and the mean
time to failure is a special case. It can be shown that
where
and
with a - see [Goyal et al. 1992]. To estimate j, we generate
sample paths ~
each consisting of m 1 1-cycles, and we use ~
X 1 to estimate
- and ~
X 2 to estimate fl. We can either let ~
independent of ~
We examine the estimation of the numerator and denominator in (6) separately.
First, if we want to estimate ff = -, then we define the function h by
h( ~
-( ~
where
Z
with
Fg. On the other hand, if we want to
estimate then we define the function h by
h( ~
where
and 1f \Delta g is the indicator function of the event f \Delta g. Thus, the standard estimator
of j is
-( ~
We will return to this example in Section 6.
6 \Delta J. M. Calvin and M. K. Nakayama
3. BASIC IDEA
Our goal now is to create a new estimator for ff. We begin by giving a heuristic
explanation of how it is constructed by considering the simple example illustrated
in

Figure

1. For simplicity, we depict a continuous sample path on a continuous
state space S. The T 1 sequence corresponds to hits to the state v, and the T 2
sequence corresponds to hits to state w. The top graph shows the original sample
path generated having regenerative 1-cycles. For this path, there are
occurrences of stopping times from sequence T 2 . To make it easier to
see the individual 2-cycles, each is depicted using a different line style. Now we
can construct new sample paths from the original path by permuting the 2-cycles,
resulting in (M possible paths. The second graph shows one such
permuted path. Here, the original third 2-cycle is now first, the original first 2-cycle
is now second, and the original second 2-cycle is now third. The new 1-cycle times
are T 0
and the new 2-cycle times are T 0
3. The
third graph contains another permuted path, in which the original second 2-cycle
is now first, the original third 2-cycle is now second, and the original first 2-cycle
is now third. The 1-cycle times are now T 00
and the new 2-cycle
times are T 00
3. Note that for each new path, the number of 1-cycles
is the same as in the original path, but the paths of some of the 1-cycles have
changed. We show in Theorem 1 that all of the paths have the same distribution.
For each possible path, we can compute an estimator of ff based on the m 1 new
1-cycles by applying the performance function h j hm1 to it. Our new estimator
is then the average over all estimators constructed.
It turns out that we do not actually have to construct all permuted paths to
calculate the value of our new estimator. The basic reason for this is that we can
break up any sample path into a collection of segments of different types. After
any permutation, the path changes, but the collection of segments does not. To
calculate our new estimator for a given (original) path, we need to determine the
different ways the segments can be put together when 2-cycles are permuted. In
particular, since we form an estimator based on the 1-cycles for every permutation,
we want to understand how 1-cycles are formed from the segments.
Another key factor that will allow us to explicitly compute our new estimator
without actually computing all permutations is that for the performance measures
we consider, the contribution of each 1-cycle to the overall estimator can be expressed
as a function of the segments in the cycle. For instance, in Example 1 with
2, we can express the contribution from each 1-cycle as the square of the sum
of contributions of the segments in the cycle.
We now examine more closely the four different types of path segments that can
arise. We focus on the example in Figure 1.
(1) The first type is a 1-cycle that does not contain a hit to w. The segments of this
type in the original path of the figure are the first and third 1-cycles; i.e., the
segment from T 1 (0) to T 1 (1) and the segment from T 1 (2) to T 1 (3). Segments of
this type never change under permutation, although they may occur at different
times. For example, the third 1-cycle in the original path appears as the fourth
1-cycle in the second permuted path. This segment is the third 1-cycle in the
first permuted path, but it occurs at a different time. The first 1-cycle in our
Using Permutations in Regenerative Simulations to Reduce Variance \Delta 7
Original
Path
Another
Permuted
One
Permuted
Path
Fig. 1. A sample path and some corresponding permuted paths.
M. Calvin and M. K. Nakayama
example always appears in the same place in all permutations.
(2) Now consider any 2-cycle in which state v is not hit, such as the third 2-cycle
in the original path in the figure. After any permutation this 2-cycle will be in
the interior of some 1-cycle. For example, the third 2-cycle in the original path
is in the interior of the fifth 1-cycle in the original path, and in the interior of
the second (resp., third) 1-cycle in the first (resp., second) permuted path.
(3) The next type of segment goes from w to v before hitting w again. No matter
how the 2-cycles are permuted, this type of segment is always the end of some
1-cycle. For example, consider the path segment from T 2 (0) to T 1 (2) in the
original sample path. In this path, the segment is the end of the second 1-
cycle. In the first permuted path, this segment is again the end of the second
1-cycle, but this new second 1-cycle is different from that in the original path.
On the other hand, this segment in the second permuted path is the end of the
third 1-cycle. In general, any segment that goes from w to v before hitting w
again will be the end of some 1-cycle in any permuted path.
(4) The final type of segment goes from v to w before hitting v again. In any
permutation, this segment will be the beginning of a 1-cycle. For example,
consider the path segment from T 1 (3) to T 2 (1) in the original sample path. In
this path, the segment is the beginning of the fourth 1-cycle. In the second
permuted path this segment is the beginning of the fifth 1-cycle. In the first
permuted path, the segment is again the beginning of the fourth 1-cycle. In
general, any segment that goes from v to w before hitting v again will be the
beginning of some 1-cycle in any permuted path.
Note that the original sample path in Figure 1 consists of segments of types
appearing in the following order: 1, 4, 3, 1, 4, 3, 4, 2, 3. In any permutation,
the segments will appear in a different order, but the collection of segments never
changes.
Recall that for each permuted path, we compute an estimate of the performance
measure ff based on the m 1 1-cycles. So now we examine how 1-cycles can be
constructed from the permutations. For the original path, we divide up the 1-
cycles into those that hit state w and those that do not. The ones that do not hit
w are type-1 segments and are unaffected by permutations.
Now we examine how permutations affect 1-cycles that hit w. These cycles always
start with a type-4 segment, followed by some number (possibly zero) of segments
of type 2, and end with a type-3 segment. For example, the fifth 1-cycle in the
original path starts with the type-4 segment from T 1 (4) to T 2 (2), is followed by the
type-2 segment from T 2 (2) to T 3 (3), and concludes with the type-3 segment from
Also, the fourth 1-cycle in the original path begins with a type-
4 segment from T 1 (3) to T 2 (1) and terminates with a type-3 segment from T 2 (1)
to T 1 (4). This characterization of 1-cycles holds not only for the original sample
path, but also for any permuted path. Moreover, for any 2-cycle that hits v in the
original path, the type-3 segment and type-4 segment in it will always be in the
same 2-cycle in any permutation, and so these two segments can never be in the
same 1-cycle since the type-3 segment will always be the end of one 1-cycle and the
segment will always be the beginning of the following 1-cycle. For example,
the type-3 segment from T 2 (0) to T 1 (2) and the type-4 segment from T 1 (3) to T 2 (1)
Using Permutations in Regenerative Simulations to Reduce Variance \Delta 9
are always in the same 2-cycle in any permutation, and as such, they are always
in successive 1-cycles. Also, in our example the first type-4 segment from T 1 (1) to
and the last type-3 segment from T 2 (3) to T 1 (5) will never be in the same
1-cycle. Any other pair of type-3 segment and type-4 segment will be in the same
1-cycle in some permutation. Thus, to construct all 1-cycles that hit w that are
possible under permutations of 2-cycles, we have to consider all valid pairings of
the type-4 and type-3 segments, and allocate the type-2 segments among the pairs.
The proofs in Sections A.2 and A.3 basically use this reasoning.
4. FORMAL DEVELOPMENT OF GENERAL METHOD
Now we more formally show how to construct our new estimator. We begin with
some new notation. Let X be the space of paths S that are
right continuous with left limits on [0; i(x)). For define a new element
and
Thus, the new path x 1 obtained by concatenating x 2 on to the end of x 1 .
Given the original sample path ~
X, which consists of m 1 1-cycles, we begin by
constructing a new sample path ~
X 0 from ~
X such that ~
equality in distribution. This is done by first taking the original sample path ~
X and
determining the number of times M 2 that the stopping times from the sequence
occur during the m 1 1-cycles. Note that if M then the path ~
X has
no 2-cycles. If M 2, then there is only one 2-cycle. Assume now that M 2 - 3.
Then for the given path ~
X, we can now look at the (M 2-cycles in the path.
We generate a uniform random permutation of the within the
path ~
, and this gives us our new sample path ~
which also has m 1 1-cycles.
More specifically, define M
~
X . If M 2 - 3, then we break up the path ~
~
where
~
is the initial path segment until the first time a stopping time from sequence T 2
occurs,
~
is the final path segment from the last time a stopping time from sequence T 2 occurs
until the end of the path, and
~
is the kth 2-cycle of the original path ~
1)) be a uniform random permutation of 1. Then we define our new
M. Calvin and M. K. Nakayama
path ~
to be
~
which is the original path ~
X with the 2-cycles permuted. Note that ~
have
the same number m 1 of 1-cycles, and we prove in Section A.1 that ~
Now for the new sample path ~
, we can calculate
which is just the estimator obtained from the new sample path ~
and is based
on m 1 1-cycles (recall that h j hm1 ). The number of possible paths ~
0 we can
construct from ~
X is N( ~
which depends on ~
X and is therefore
random. We label these paths ~
each of which has the
same distribution as ~
X, and for each one we construct b ff( ~
finally define
our new estimator for ff to be
e
h( ~
Another way of looking at our new estimator is as follows. We first generate the
original path ~
X and use it to construct the N( ~
X) new paths ~
X (N) . We
then choose one of the new paths at random uniformly from ~
this be ~
X, we can think of b
standard estimator of ff since
it has the same distribution as b
X). Then we construct our new estimator e
to be the conditional expectation of b
respect to the uniform random
choice of ~
given the original path ~
X. That is, if E   denotes expectation with
respect to choosing ~
X 0 from the uniform distribution on ~
then we write
e
Assuming that E[-ff( ~
the new estimator has the same mean as the
original since
because ~
X. Moreover, decomposing the variance by conditioning on ~
us
implies that the variance of the new estimator
e
is no greater than that of the original estimator -
X). This
calculation, combined with the fact that ~
(which will be proved in Section
A.1), establishes the following theorem.
Theorem 1. Let T 1 and T 2 be two distinct sequences of stopping times, and
construct the estimator e
X) defined by (11). Assume that E[bff( ~
Using Permutations in Regenerative Simulations to Reduce Variance \Delta 11
X)], and
and so the mean squared error of our new estimator e
X) is no greater than
that of the original estimator b ff( ~
X). Strict inequality is obtained in (12) unless
In Theorem 1 we see that there is no variance reduction when for every possible
original sample path ~
X, the value of the function h in (1) is unaffected by permutations
of the 2-cycles. For example, this is the case in Example 1 with
since
h( ~
Z
Z
Z
Z
Z
and so e
X). Similarly, by choosing g(x) j 1, we see that permuting
2-cycles does not alter the estimator for E[-(1)]. Thus, our method has no effect
on the standard ratio estimator for steady-state performance measures ff that can
be expressed as
However, for p ? 1 in Example 1, we have in general that h( ~
so typically e
X). Also, we usually have that the standard time-average
variance estimator in Example 2 for a regenerative simulation will differ from the
new estimator defined by (11). Finally, applying the above idea separately to the
numerator and denominator in the ratio expression for the mean cumulative reward
until hitting some set of states F as in Example 3 will result in a new estimator.
5. ESTIMATING THE SECOND MOMENT OF CUMULATIVE CYCLE REWARD
For our new estimator e
X) to be computationally efficient, we need to calculate
explicitly the conditional expectation in (11) without having to construct all
possible permutations. We first do this for Example 1 with
and our standard estimator of ff is
where we have dropped the dependence of Y on g to simplify the notation. Our
new estimator of ff is then
e
M. Calvin and M. K. Nakayama
is the same as Y (k) except that it is for the sample path ~
than ~
X.
Now to explicitly calculate (15) in this particular setting, we will divide up
the original path into segments using the approach described in Sections 3 and
4. We need some new notation to do this. For our two sequences of stopping
times denote the set of indices of the
1-cycles in which a T 2 stopping time occurs, and define the complementary set
specifically, H(1;
jg. We analogously define the set H(2; 1) with the roles of
some lg, which is the first occurrence of a stopping time
from sequence T 2 after the 1)st stopping time from the sequence T 1 . Similarly
define e
some lg, which is the last occurrence
of a stopping time from sequence T 2 before the kth occurrence of the stopping-time
sequence T 1 . Then, for k 2 H(1; 2), we let
which is the contribution to Y (k) until a stopping time from sequence T 2 occurs,
and let
Y
Z
e
which is the contribution to Y (k) from the last occurrence of a stopping time from
sequence T 2 in the kth 1-cycle until the end of the cycle. Also, for l 2 J(2; 1), let
Y 22
Z
which is the integral of g(X(t)) over the lth 2-cycle in which there is no occurrence
of a stopping time from sequence T 1 . We now define B k ae J(2; 1) to be the set of
indices of those 2-cycles that do not contain any occurrences of the stopping times
from the sequence T 1 and that are between It then follows
that for k 2 H(1; 2),
Hence,
\Theta Y 12
A
Using Permutations in Regenerative Simulations to Reduce Variance \Delta 13
Y 22 (l)
Y 22 (l)
In the last expression, the first term does not change if we replace the original
sample path ~
X with the new sample path ~
the last term does change.
In Section A.2, we compute explicitly the conditional expectation of (16), when ~
is replaced with ~
with respect to a random permutation given the original path
~
X.
The expression for this involves some more notation. Define
and
Finally, define fi l to be the lth smallest element of the set H(1; 2) for
and define fi
is the index in H(1; 2) that occurs just before k if k is
not the first index and is the last element in H(1; 2) if k is the first element. The
following theorem is proved in Section A.2. (Pseudo-code for our estimator is given
in

Appendix

B.)
Theorem 2. Suppose we want to estimate ff defined in (13), and assume that
our new estimator is given by e
and otherwise by
e
\Theta Y 12

Y
Y 22
Y 22
l6=m
Y 22 (l)Y 22 (m)C A : (17)
The estimator satisfies E[eff( ~
X) is
the standard estimator of ff as defined in (14).
5.1 A New Estimator for the Time-Average Variance
We can use Theorem 2 to construct a new estimator for the variance parameter in a
regenerative simulation of the process X . We start by first giving a more complete
explanation of Example 2 in Section 2.
cost function. Define
Z tf(X(s)) ds:
14 \Delta J. M. Calvin and M. K. Nakayama
Since X is a regenerative process, there exists some constant r such that r t ! r
as Theorem 2.2 of [Shedler 1993]). Also, r satisfies
the ratio formula Assuming that E[Z(f ;
exists a finite positive constant oe such that
oe
as t !1. The constant oe 2 is called the time-average variance of X and is given in
(3). Given the central limit theorem described by (18), construction of confidence
intervals for r therefore effectively reduces to developing a consistent estimator for
oe 2 . The quality of the resulting confidence interval is largely dependent upon the
quality of the associated time-average variance estimator.
The standard consistent estimator of oe 2 is b
X) defined in (4). Note
that b oe 2 ( ~
X) can be expressed as
Now we define our new estimator e
X) to be the conditional expectation of b oe 2 ( ~
with respect to a random permutation of 2-cycles, given the original sample ~
X .
Hence, letting b r 0 , Y be the corresponding values of b r, Y (f \Gamma
k), and -(k) for the sample path ~
we get that
\Theta P m1

since
k=1 -(k) is independent of the permutation of
2-cycles. Also, observe that
Z
r
is independent of the permutation of 2-cycles, so
\Theta P m1

i.e., we can replace b r 0 with b r. The following is a direct consequence of Theorem 2.
Corollary 3. Suppose we want to estimate oe 2 defined in (3), and assume that
our new estimator e oe 2 ( ~
X) is given by (19), where the
numerator is as in (17) with the function
r. The estimator satisfies
X)] and Var[eoe 2 ( ~
X)].
5.2 Continuous-Time Markov Chains
We now consider the special case of an irreducible, positive-recurrent, continuous-time
Markov chain on a countable state space S having
generator matrix
be the embedded discrete-time Markov chain, and
be the sequence of random holding times of the continuous-time
Markov chain; i.e., Wn is the time between the nth and (n 1)st transitions
Using Permutations in Regenerative Simulations to Reduce Variance \Delta 15
of X . Define A
which is the time of
the nth transition. It is well known that conditional on V , the holding time in
state Vn is exponentially distributed with mean 1=-(Vn ) and that W i and W j are
(conditionally) independent for i 6= j.
Assume that the sequences of stopping times T 1 and T 2 correspond to hitting
times to fixed states v 2 S and w 2 S, respectively, with w 6= v, and assume that
specifically, define - 1
which is the sequence of hitting times to state v for
the discrete-time Markov chain. Similarly, define -
and
Suppose that we want to estimate ff as defined in (13), and our standard estimator
of ff is given in (14). Now note that
Z
Using discrete-time conversion [Hordijk et al. 1976; Fox and Glynn 1990] gives us
\Theta W 2
since are conditionally independent given V . For a function f
!, let
which is the cumulative reward over the kth cycle for the discrete-time chain V ,
and define the functions to be g 1
Therefore, we get
and
To create our new estimator of ff, we then compute the conditional expectation of
X) with respect to a random uniform permutation of 2-cycles given the original
M. Calvin and M. K. Nakayama
path ~
X. Define
The last term in (20) is independent of permutations of 2-cycles, and so we get the
following expression for e -
, which follows from Theorems 1 and 2.
Theorem 4. Suppose X is an irreducible, positive-recurrent, continuous-time
Markov chain on a countable state space S, and we want to estimate ff defined in
(13). Assume that T 1 and T 2 correspond to the hitting times to states v and w,
respectively, with w 6= v. Assume that E[Y (g our new estimator
is given by e -
X) is defined by (21), which can be computed from (17) with the function
1 . The estimator satisfies E[e -
is the standard estimator of ff as defined in (14).
If we had instead first converted to discrete time and then computed -
the discrete-time Markov chain and its conditional expectation with respect to the
permutation, we would have obtained e -
X) as our estimator for ff. However,
since E[ -
function g 6= 0, E[e -
X) is biased.
On the other hand, our estimator e -
X) is unbiased.
6. EXPECTED CUMULATIVE REWARD UNTIL HITTING A SET
Recall that we can express the expected cumulative reward until a hitting time
given in (5) as the ratio in (6), and the standard estimator of j is defined
in (7). Also, recall that the numerator - is estimated using the sample path ~
and the denominator fl is estimated from path ~
We will examine both the cases
when ~
are independent.
In the context of estimating the mean time to failure of highly reliable Markovian
systems, Goyal, Shahabuddin, Heidelberger, Nicola, and Glynn [1992] and
Shahabuddin [1994] estimate - and fl independently; i.e., ~
are indepen-
dent. This is useful because then different sampling techniques can be applied to
estimate the two quantities. In particular, fl is the probability of a rare event and
so it is estimated using importance sampling. On the other hand, we can efficiently
estimate - using naive simulation (i.e., no importance sampling). Below, we do not
apply importance sampling to estimate fl, but one can also derive a new estimator
of fl when using importance sampling.
Our new estimator of j is defined as
e j( ~
e
-( ~
e
Using Permutations in Regenerative Simulations to Reduce Variance \Delta 17
where
e
-( ~
-( ~
e fl( ~
Now to explicitly calculate the numerator and denominator, we will divide up the
original path into segments using the approach described in Sections 3 and 4. We
need some new notation to do this. For k 2 H(1; 2), let
I 12
I 22
I
with
Fg. Hence, I 12 (k) (resp., I 21 (k))
is the indicator of whether the set F is hit in the initial 1-2 segment (resp., final
2-1 segment) of the 1-cycle with index k 2 H(1; 2). Also, I 22 (l) is the indicator
whether the set F is hit in the 2-cycle with index l 2 J(2; 1).
We first consider the denominator fl. To derive the new estimator of fl from
permuting the 2-cycles, we first write
The first term on the right-hand side is independent of permutations of the 2-cycles.
For the second term we note that for k 2 H(1; 2),
I 12 (k); max l2Bk
I 22 (l); I 21 (k)
Thus,
e
I 12 (k); max
I 22 (l); I 21 (ae(k))
where ae(k) is the index of the I 21 variable that follows the I 12 (k) variable after
a permutation of the 2-cycles, and B 0
l is the same as B l except that B 0
l is after a
permutation. We work out in Section A.3 the conditional expectation appearing
above.
We now examine the estimation of -. Note that the standard estimator of -
satisfies
-( ~
The first term is not affected by permuting the 2-cycles, but the second term is.
M. Calvin and M. K. Nakayama
For
Y
j!l
Y
where
Z T (2)
F (l)-T2 (l)
Z T (1)
e
Hence, to compute the new estimator, we need to compute the conditional expectation
of the second term in (24), which we can do by using the representation for
given in (25); this is done in Section A.3.
To present what the new estimator actually works out to, we need some more
notation. Define
Also, define
I 22 (l);
r
D 22 (l)I 22 (l):
Then we have the following result, whose proof is given in Section A.3.
Theorem 5. Suppose we want to estimate j in (6), and assume that
1. Then, our new estimator is given by (22) where
Using Permutations in Regenerative Simulations to Reduce Variance \Delta 19
(i). e -( ~
-( ~
e
-( ~
r
(ii). e fl( ~
The estimators e
-( ~
-( ~
-( ~
-( ~
-( ~
are the standard
estimators of - and fl, respectively.
In Theorem 5 the variables used in part (i) are defined for the sample path ~
and the variables in part (ii) are for the sample path ~
. For example, h 12 in
part (i) is the cardinality of the set H(1; 2) for the path ~
in part (ii) it is
the same but instead for the path ~
Theorem 5 shows that our new estimator for j has unbiased and lower-variance
estimators for both the numerator and denominator, but the effect on the resulting
ratio estimator is more difficult to analyze rigorously. Instead, we now heuristically
examine the bias and variance of the ratio estimator.
To do this, we generically let -
- fl, and -
-fl be estimators of -, fl, and j,
respectively. Then using first- and second-order Taylor series expansions, we have
the following approximations for the bias and variance of - j:
Cov
and
see p. 181 of [Mood et al. 1974].
We now use these approximations to analyze the standard and new estimators
for j. First, consider the case when ~
are independent. Then b
-( ~
are independent, so Cov
-( ~
Similarly, e
-( ~
are independent, so Cov
-( ~
it follows from Theorem 5
and (26) and (27) that
J. M. Calvin and M. K. Nakayama
and
where we use the notation a - b to mean that a is approximately no greater than b.
Hence, the mean square error of e j( ~
approximately no greater than that
of b j( ~
In the case when ~
-( ~
Also, we have that
Cov
-( ~
-( ~
by Lemma 2.1.1 of [Bratley et al. 1987]. But since the variances of the new estimators
of the numerator and denominator are smaller than those for the original esti-
mators, we cannot compare the biases and variances of b j( ~
However, we examine this case empirically in Section 8 and find that there is a
variance reduction and smaller mean squared error.
7. STORAGE AND COMPUTATION COSTS
We now discuss the implementation issues associated with constructing our new
estimator e ff( ~
X) given in (17) for the case when ff is defined in (13). First note that
the first term in the second line of (17), excluding the factor 2=(h
Y
Y
Also, the last term in the last line of (17) satisfies
l6=m
Y 22 (l)Y 22 (m) =@ X
Y 22 (k)A\Gamma@ X
Y 22
Hence, to construct our estimator e
X), we need to calculate the following quantities

-the sum of the Y (k) 2 over the 1-cycles k 2 J(1; 2);
-the sums of the Y 12 (k), Y 21 (k), Y 12 over the 1-cycles k 2 H(1; 2);
-the sum of the Y 12 (k)Y 21 (/(k)) over the 1-cycles k 2 H(1; 2);
-the sums of the Y 22 (k) and Y 22 (k) 2 over the 2-cycles k 2 J(2; 1).
To compute these quantities in a simulation, we do not have to store the entire
sample path but rather only need to keep track of the various cumulative sums as
the simulation progresses. Also, the amount of storage required is fixed and does
not increase with the simulation run length. Therefore, compared to the standard
estimator, the new estimator can be constructed with little additional computational
effort and storage. (Pseudo-code for this estimator is given in Appendix B.)
Using Permutations in Regenerative Simulations to Reduce Variance \Delta 21
A similar situation holds when estimating j using the estimator defined in Theorem
5.
We conclude this section with a rough comparison of the work required for the
new estimator with that for the standard regenerative method when estimating ff
given in (13). Let W s be the (random) amount of work to generate a particular
sample path of m 1 1-cycles in a discrete-event simulation, where W s includes the
work for the random-variate generation, determining transitions, and appropriately
updating data structures needed in the sample-path generation. This quantity is
the same for our new method and the standard method.
We now study the work needed for the output analysis required by the standard
regenerative method. After every transition in the simulation, we need to update
the value of the current cycle-quantity Y (g; k); see (2). Let ' 1 denote this (deter-
ministic) amount of work, and if there are N total transitions in the sample path,
then the total work for updating the Y (g; during the entire simulation is N' 1 . At
the end of every 1-cycle, we have to square the current cycle quantity Y (g;
add it to its accumulator; see (14). Let ' 2 denote this (deterministic) amount of
work required at the end of each 1-cycle, and since there are m 1 1-cycles along the
path, the total work for accumulating the sum of the Y (g; Therefore,
the cumulative work (including sample-path generation and output analysis) for
the standard regenerative method is
Now we determine the amount of work needed for the output analysis of our new
permutation method. By examining the pseudo-code in Appendix B, we see that
after every transition, a single accumulator is updated. Every time a stopping time
from either sequence T 1 or T 2 occurs, we compute either a square or a product of
two terms and update at most three accumulators, and the amount of work for
this is essentially at most 3' 2 . Since the number of times this needs to be done is
the cumulative work for our permutation method is
Therefore, the ratio of cumulative work of our permutation method relative to the
standard regenerative method is
Typically in a regenerative simulation, the amount of time W s required to generate
the sample path is much greater than the time needed to perform the output
analysis, and so RW will usually be close to 1. Hence, the overhead of using our
method in a simulation will most likely be very small, and this is what we observed
in our experimental results in the next section.
8. EXPERIMENTAL RESULTS
Our example is based on the Ehrenfest urn model. The transition probabilities for
this discrete-time Markov chain are given by P
s
s:
22 \Delta J. M. Calvin and M. K. Nakayama
In our experiments we take 8. The stopping-time sequences T 1 and T 2 for
our regenerative simulation correspond to hitting times to the states v and w,
respectively, and so state v is the return state for the regenerative simulation.
We ran several simulations of this system to estimate two different performance
measures: oe 2 , which is the time-average variance constant from Section 5.1, and j,
which is the mean hitting time to a set F from Section 6. For each performance
measure, we ran our experiments with several different choices for v, and for each v,
we examined all possible choices for w. Choosing no effect on the resulting
estimator, so this corresponds the standard estimator. We ran 1,000 independent
replications for each choice of v and w. Tables 1-3 and 5-6 present the results from
estimating the two performance measures, giving the sample average and sample
variance of our new estimator over the 1,000 replications. The average cycle lengths
change with different choices of v; in order to make the results somewhat comparable
across the tables, we changed the number of simulated cycles for each case so that
the total expected number of simulated transitions remains the same. For Table 1,
corresponding to simulated 1,000 cycles, and a greater number for the
other tables. For example, the expected cycle length is 3:5 times as long for state
1 as for state 2, so in Table 2, we simulated 3,500 cycles. Since our new estimator
reduces the variance but at the cost of extra computational effort, we also compare
the efficiencies (inverse of the product of the variance and the time to generate the
estimator) of our new estimator and the standard one, as suggested by Hammersley
and Handscomb [1964] and Glynn and Whitt [1992].
8.1 Results from Estimating Variance
We first examine the results from estimating the time-average variance oe 2 with cost
function performed 3 experiments, corresponding to return states
and 4, and these results are given in Tables 1-3, respectively.
The transition probabilities are symmetric around state 4 (the mean of the binomial
stationary distribution), so our first choice of return state
is fairly far from the mean. Notice that the variability of the variance estimator is
smaller with w near the mean state 4, and that the variance reduction is greater
v. The reason for this is that the excursions from v that go below 1 have
little variability; because of the strong restoring force of the Ehrenfest model, such
excursions tend to be very brief. On the other hand, excursions that get as far
as the mean are likely to be quite long (and thus the contribution to the variance
estimator tends to have large variability). In the second table we ran the same
experiment with obtained similar results.
In

Table

3 we examine the same model, but now with our return state v chosen
to be the stationary mean, 4. The first thing to notice is that, compared with the
other choices of the return state, the variance reduction is relatively small. State 4
is the best return state in the sense of minimizing the variance of the regenerative-
variance estimator. Therefore, for this example, it appears that our estimator is a
significant improvement over the standard regenerative estimator if the standard
regenerative estimator is based on a relatively "bad" return state. However, if one
is able to choose a near-optimal return state to begin with, our estimator yields
a modest improvement. (Unfortunately, there are no reliable rules for choosing a
priori a good return state.) Comparing the three tables, we see that the minimum
Using Permutations in Regenerative Simulations to Reduce Variance \Delta 23

Table

1. Estimating variance with
w Avg. of eoe 2 Sample Var.

Table

2. Estimating variance with 2.
w Avg. of eoe 2 Sample Var.
variability does not change much across tables (0:17 for Table 3, and 0:18 for the
other tables). This example suggests that it may be possible to compensate for a
bad choice of v by an appropriate choice of w.
Finally, we illustrate the computational burden of our method. Table 4 shows
the work required for the results obtained in Table 1. The first column gives, for
each choice of w, the relative work (CPU time for generating the sample path and
output analysis) required for our new estimator, that is, the CPU time with our new
estimator divided by the CPU time for the standard regenerative method. (The
row with corresponds to the standard regenerative method, so all entries are
the other entries are normalized with respect to these.) The second column gives
the relative variance; that is, the sample variance of our new estimator divided by
the sample variance for the standard regenerative estimator. The last column gives

Table

3. Estimating variance with
w Avg. of eoe 2 Sample Var.
M. Calvin and M. K. Nakayama

Table

4. Comparison of efficiency,
w Relative Work Relative Var. Relative Efficiency
the relative efficiency; that is, the inverse of the product of the relative work from
column 1 and the relative variance from the second column.
Notice that in the cases where the variance reduction is small, the increase in
work is also small. More work is needed when there is a larger variance reduction,
but this is still no more than a few percent increase. Note that in the best case
the efficiency was improved by nearly a factor of eight. Because little
additional work is needed by our method, within each of the other tables, the run
times are approximately the same for the different values of w. Therefore, one can
roughly approximate the relative efficiencies for Tables 2 and 3 as the ratio of the
sample variances for states v and w.
It should also be noted that the Ehrenfest model considered here is very simple
compared with typical simulation models. The additional work required to compute
our new estimators is independent of the model, and so if the work to generate the
sample path is much larger than for the Ehrenfest model, the relative increase in
work would be correspondingly small.
8.2 Results from Estimating Hitting Times to a Set
We now consider estimating j, which is the mean hitting time to a set of states F
starting from a state v for our Ehrenfest model with reward function g(x) j 1.
We take which is hit infrequently. Tables 5 and 6 show our results from
generating 1,000 independent replications for In
each replication, we generated a sample path which we used to compute the new
estimators for both the numerator and denominator. Hence, using the terminology
of Section 6, we let ~
. For each path, we generated 1,000 cycles for
and 2,500 cycles for 4, so that the expected sample path length is the same
in each case. We calculated (i.e., not using simulation) the theoretical values to be
in addition to examining the
variance of our new estimator, we can also study the mean squared error. Note
that the theoretical value of j depends on the starting state v.
First observe that there is no change in the estimates of j and their variances
for certain choices of w. This is due to the fact that permuting w-cycles in these
cases has no effect on the estimators of either the numerator or denominator. For
the other values of w, the (relative) variance reduction is significantly greater when

Table

than when (Table 5). In addition, although the absolute
bias is the greatest for both choices of v, its magnitude is quite small,
Using Permutations in Regenerative Simulations to Reduce Variance \Delta 25

Table

5. Estimating expected hitting time to state 7 with
w Avg. of ej Sample Var. MSE

Table

6. Estimating expected hitting time to state 7 with 2.
w Avg. of ej Sample Var. MSE
and when examining the mean squared error, the variance reduction overwhelms
the effect of the bias.
We now explain why the choice of w that results in the most variance reduction
is In the original sample path without permuting the w-cycles, the set F
is hit a certain number of times. By permuting the w-cycles, we get a variance
reduction if there are some v-cycles that hit w but not F and if within a particular
v-cycle that hits F , there is more than one w-cycle that hits F and does not hit v.
Permuting the w-cycles then can distribute the hits to F to more of the v-cycles.
The amount of variance reduction in estimating fl is largely determined by the
difference between the maximum and minimum number of v-cycles that hit F from
permuting the w-cycles. Choosing results in no variance reduction
because we are working with a birth-death process and so the process always hits
F no later than it hits w within a v-cycle, and so permuting w-cycles has no effect.
Of the remaining choices for maximizes the number of w-cycles that hit
F , hence the largest variance reduction. Therefore, in general, we suggest that the
state w should be chosen so that w 62 F and it is as "close" as possible to the set
F to maximize the number of w-cycles that hit F .
9. DIRECTIONS FOR FUTURE RESEARCH
We are currently investigating how to construct confidence intervals based on our
new permuted estimators. This is a difficult problem because of the complexity of
the estimators. Another area on which we are currently working is determining how
to choose the two sequences of regenerative times T 1 and T 2 when there are more
than two possibilities. For example, this arises when simulating a Markov chain,
26 \Delta J. M. Calvin and M. K. Nakayama
since successive hits to any fixed state form a regenerative sequence. We explored
this to some degree experimentally in Section 8, but further study is needed.

ACKNOWLEDGMENTS

The authors would like to thank the Editor in Chief and the two anonymous referees
for their helpful comments on the paper.




A. PROOFS
A.1 Proof of Theorem 1
We need only prove that ~
that is, the paths have the same distribution when
2-cycles are permuted. Recall M 2 is the number of times that stopping times in
the sequence T 2 occur by time
as in (8)-(10). For the path ~
to be the number of times a
stopping time from sequence T 1 occurs in ~
to be the number of times that stopping times from the sequence T 1 occur outside
of the 2-cycles (we do not need to bother with the case since then we do
not change the path).
If f is a (measurable) function mapping sample paths to nonnegative real values,
then
i-m1;n
i-m1;n
for some (measurable) functions f n;i , and so it suffices to show that
~
is invariant under 2-cycle permutations, where the A i are (measurable) sets. Note
that
~
\Theta P (M
Given M 2 and L, the initial 1-2 segment (i.e., ~
final 2-1 segment (i.e., ~
are conditionally independent of the 2-cycles, so the last probability can be written
\Theta P
In examining the effect of permutations of the 2-cycles, we need consider only the
Using Permutations in Regenerative Simulations to Reduce Variance \Delta 27
last probability, which we rewrite
~
Since we are interested in the effect of permutations, we only look at the numerator
of the last expression:
and we are finally left with the task of showing that for any permutation oe,
P@ ~
But this follows from the fact that
~
~
Therefore, ~
, and the theorem is proved. Notice that we only used the
conditional exchangeability of the cycles, and not the full independence.
A.2 Proof of Theorem 2
Recall that in Section 5 we defined
for the original sample path ~
X. Using a permuted path ~
instead of the
original path ~
X in (16), we get
\Theta Y 0
A
22 (l) +@ X
22 (l)A1
28 \Delta J. M. Calvin and M. K. Nakayama
where the Y variables and sets are the same as the Y; H; J; B variables
and sets, respectively, with ~
replacing ~
X in (16). Recall that E   is the conditional
expectation operator corresponding to a random (uniform) permutation of 2-cycles
(as was done when constructing the path ~
X 0 from ~
X) given the original sample
path ~
X. Also, recall that we define our new estimator to be
e
which we will now show is equivalent to (17).
First note that by our construction of the path ~
X 0 from ~
X, the first term in (28)
does not change when replacing ~
X with ~
\Theta Y 0
\Theta Y 12
Now we compute the conditional expectation of the second term of (28). We can
assume that h 12 - 3. Let ae(k) be the index of the Y 21 segment that follows the
segment after a permutation of the 2-cycles. Note that ae(k) 6= /(k) since
Y 12 (k) and Y 21 (/(k)) are always in the same 2-cycle, no matter how the 2-cycles
are permuted. Any of the other h indices from H(1; 2) are equally likely,
however, so that
Y
For the second summand in the second term of (28), we have
Y 22 (l)7 5
Y 22 (l)7 5
Y 22 (l)7 5
Y 22 (l)7 5
Y
Y 22 (l)
Using Permutations in Regenerative Simulations to Reduce Variance \Delta 29
Y 22 (l)5
Y
Y 22 (l)5
Y 22 (l): (31)
For the third summand in the second term of (28), we have
Y 22 (l)
Y 22 (l) 1 fl2Bk gA3
l6=m
Now note that
Y 22 (l) 2
Y 22
Also,
l6=m
Y 22 (l)Y 22 (m) 1 fl2Bk;m2Bk g7 5
l6=m
Y 22 (l)Y 22 (m)
l6=m
Y 22 (l)Y 22 (m)
Now use the fact that
which follows from
Lemma 6. Suppose that p white balls, numbered 1 to p, are placed along with q
black balls into p+q boxes arranged in a line, with each box getting exactly one ball.
Apply a uniformly chosen random permutation to the balls. Then the probability
that ball 1 and ball 2 are not separated by a black ball is 2=(2
To apply this lemma to (32), we let be the number of 2-cycles that
include one of the T 1 stopping times, and p be the number of remaining 2-cycles.
Proof. Let D be the number of boxes in between the boxes containing ball 1
M. Calvin and M. K. Nakayama
and ball 2, and let L i be the box number containing ball i
Given that the probability that balls 1 and 2 are not separated by at least
black ball is the probability that all q black balls are chosen from the p+q
boxes that are not between ball 1 and ball 2, which is
Thus the desired probability is
Hence,
Y 22 (l)
l6=m
Y 22 (l)Y 22 (m): (33)
Using Permutations in Regenerative Simulations to Reduce Variance \Delta 31
Finally, putting together (29), (30), (31), and (33), we get that e
X) is as in (17).
The unbiasedness and variance reduction follow directly from Theorem 1.
A.3 Proof of Theorem 5
We first need the following result.
Lemma 7. Suppose that q black balls, r white balls, and p red balls are placed at
random in q +p+ r boxes arranged in a line (one ball per box). The probability that
there are no white balls in an interval formed by two particular black balls (or the
start or end of the boxes) is q=(q
Proof. Count boxes from the left until a non-red ball is encountered. The
desired probability is the probability that the first non-red ball is a black ball. Since
of the non-red balls q are black and r are white, this probability is q=(q
Now we prove Theorem 5. First recall our definitions of I 12 (k), I 22 (l), I 21 (k),
22 (l), and D 21 (k) given in Section 6. Also recall
(23). The first term in (23) is independent of the permutations, and so we now
consider the second term.
Note that I 12 (i), I 22 (j), and I 21 (k) are independent for any i; j; k. Then
I 12 (k); max
I 22 (l); I 21 (ae(k))
I 12 (k); max
I 22 (l); I 21 (ae(k))
I 12 (k); max
I 22 (l); I 21 (ae(k))
I 22
I 22
where we have used Lemma 7.
We now examine the estimation of -. Recall (24) and (25). The first term in (24)
is not affected by permuting the 2-cycles, but the second term is, and so we now
M. Calvin and M. K. Nakayama
examine the second term. First note that
since the D 12 (k) are unaffected by permutations of 2-cycles. Also,
Y
Y
Y
I 22
Finally, after the permutation of the 2-cycles, define R(k) to be the number of 2-
cycles in J(2; 1) that immediately follow the path segment corresponding to D 12 (k),
and let be the indices of those 2-cycles in J(2; 1) that
immediately follow the path segment corresponding to D 12 (k) in the order they
appear. Then
Y
j!l
D 22 (ffi l (k))
Y
j!l
Using Permutations in Regenerative Simulations to Reduce Variance \Delta 33
Lemma 8.
D 22 (ffi l
r
Proof. Suppose m balls are placed in n - m boxes in a line. Let Z be the
number of empty boxes on the left end. Then
for
and
(substitute
use the identity
To get the mean number of D 22 's that do not hit F , we use the above formula with
Finally, putting together Lemma 8 and (34), (35), and (36), we get our new
estimator for -. The unbiasedness and variance reduction of our two estimators
34 \Delta J. M. Calvin and M. K. Nakayama
follow directly from Theorem 1.
B. PSEUDO-CODE
B.1 Estimator in Theorem 2
Below is the pseudo-code for the estimator in Theorem 2. (The pseudo-code for estimator
in Theorem 5 is similar.) Note that it is specifically for a discrete-event simulation
of a continuous-time process. Discrete-event processes
can be handled by letting the inter-event time \Delta always be 1. The estimator is
denoted by alpha.
number of regenerative 1-cycles to simulate;
k / 0; // counter for number of regenerative 1-cycles
cardinality of H(1; 2)
sumy12 / 0; // sum of the Y 12 (k) over H(1; 2)
sum of the Y 21 (k) over H(1; 2)
sumy22 / 0; // sum of the Y 22 (k) over J(2; 1)
sumysq / 0; // sum of the Y
sumy12sq / 0; // sum of the Y 12
sumy21sq / 0; // sum of the Y 21
sumy22sq / 0; // sum of the Y 22
sum of the Y 12 (k)Y 21 (/(k)) over H(1; 2)
accum / 0; // accumulator for Y over the current segment
laststoptime / or 2) of last stopping time to occur
generate initial state x;
do while (k ! m)
generate inter-event time \Delta;
accum
generate next state x;
occurs at current time t) then
else
lasty21 / accum;
endif
accum / 0;
endif
occurs at current time t) then
Using Permutations in Regenerative Simulations to Reduce Variance \Delta
else
endif
else
sumy22
endif
accum / 0;
endif
alpha / sumysq/m;
else
endif



--R

A Guide to Simulation
A new variance-reduction technique for regenerative simulations of Markov chains
Practical Nonparametric Statistics
Simulating stable stochastic systems
methods: Another look at the jackknife.

A joint central limit theorem for the sample mean and regenerative variance estimator.
The asymptotic efficiency of simulation estimators.

Monte Carlo Methods.

Topics on Regenerative Processes.
Regenerative Phenomena.
Introduction to the Theory of Statistics
Some invariance principles relating to jackknifing and their role in sequential analysis.
Approximation Theorems of Mathematical Statistics.
Importance sampling for highly reliable Markovian systems.
Regenerative Stochastic Simulation.
--TR
A guide to simulation (2nd ed.)
Discrete-time conversion for simulating finite-horizon Markov processes
The asymptotic efficiency of simulation estimators
A Unified Framework for Simulating Markovian Models of Highly Dependable Systems
Importance sampling for the simulation of highly reliable Markovian systems

--CTR
James M. Calvin , Peter W. Glynn , Marvin K. Nakayama, On the small-sample optimality of multiple-regeneration estimators, Proceedings of the 31st conference on Winter simulation: Simulation---a bridge to the future, p.655-661, December 05-08, 1999, Phoenix, Arizona, United States
James M. Calvin , Marvin K. Nakayama, Exploiting multiple regeneration sequences in simulation output analysis, Proceedings of the 30th conference on Winter simulation, p.695-700, December 13-16, 1998, Washington, D.C., United States
James M. Calvin , Marvin K. Nakayama, Output analysis: a comparison of output-analysis methods for simulations of processes with multiple regeneration sequences, Proceedings of the 34th conference on Winter simulation: exploring new frontiers, December 08-11, 2002, San Diego, California
Michael A. Zazanis, Asymptotic Variance Of Passage Time Estimators In Markov Chains, Probability in the Engineering and Informational Sciences, v.21 n.2, p.217-234, April 2007
James M. Calvin , Marvin K. Nakayama, SIMULATION OF PROCESSES WITH MULTIPLE REGENERATION SEQUENCES, Probability in the Engineering and Informational Sciences, v.14 n.2, p.179-201, April 2000
Wanmo Kang , Perwez Shahabuddin , Ward Whitt, Exploiting regenerative structure to estimate finite time averages via simulation, ACM Transactions on Modeling and Computer Simulation (TOMACS), v.17 n.2, p.8-es, April 2007
James M. Calvin , Marvin K. Nakayama, Improving standardized time series methods by permuting path segments, Proceedings of the 33nd conference on Winter simulation, December 09-12, 2001, Arlington, Virginia
James M. Calvin , Peter W. Glynn , Marvin K. Nakayama, The semi-regenerative method of simulation output analysis, ACM Transactions on Modeling and Computer Simulation (TOMACS), v.16 n.3, p.280-315, July 2006
variance reduction;permutations;efficiency improvement;regenerative simulation