\documentclass[acmtomacs]{acmtrans2m}
\usepackage{graphicx}
\usepackage{color}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{conjecture}[theorem]{Conjecture}
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{lemma}[theorem]{Lemma}
\newdef{definition}[theorem]{Definition}
\newdef{remark}[theorem]{Remark}
\newcommand{\Exp}{\mbox{Exp}}
\newcommand{\Unif}{\mbox{Unif}}
\newcommand{\myrule}{\rule[1ex]{\textwidth}{0.15mm}}
\newenvironment{algorithm}{\begin{tabbing}\myrule\\%
zzzz\=zzzz\=zzzz\=zzzz\=zzzz\=zzzz\=zzzzzzzzzzzz\=\kill}{\myrule\end{tabbing}}
\newcommand{\X}{{\cal X}}
\renewcommand{\baselinestretch}{2}
\bibliographystyle{acmtrans}
\markboth{D. J. Murdoch and G. Takahara}{Perfect Sampling for Queues and Network Models}
\title{Perfect Sampling for Queues and Network Models}
\author{DUNCAN J. MURDOCH \\ University of Western Ontario \and GLEN TAKAHARA \\ Queen's University}
\begin{abstract}
We review Propp and Wilson's [1996] CFTP algorithm and
Wilson's [2000] ROCFTP algorithm. We then use these to construct
perfect samplers for several queueing and network models: Poisson
arrivals and exponential service times, several types of customers,
and a trunk reservation protocol for accepting
new customers; a similar protocol on a network switching model; a queue with
a general arrival process; and a queue with both general arrivals and
service times. Our samplers give effective ways to generate random samples
from the steady-state distributions of these queues.
\end{abstract}
\category{I.6.8}{Simulation and Modelling}{Types of Simulation}[Monte Carlo]
\category{I.6.3}{Simulation and Modelling}{Applications}
\terms{Algorithms}
\keywords{Coupling from the past, perfect simulation, queues, networks}
\begin{document}
\begin{bottomstuff}
Author's addresses: D. J. Murdoch, Dept. of Stat. and Act. Sci., University of Western
Ontario, London, Ontario, Canada N6A 2C8; G. Takahara, Dept. of Math and Stats, Queen's University, Kingston, Ontario, Canada K7L 3N6.
\end{bottomstuff}
\maketitle
\section{Introduction}
The aim of {\em perfect sampling} is to take approximate simulation
methods (e.g. simple Markov chain Monte Carlo) and to simulate from the
limiting distribution, usually using some sort of coupling method. In
this paper we demonstrate how to apply these methods to queues and
network models.
There is a huge literature on simulation of queues and network models in
general. In particular, Andrad\'{o}ttir and Ott \shortcite{Andradottir95} and
Andrad\'{o}ttir and Hosseini-Nasab \shortcite{Andradottir03} described
methods of simulation of Markovian models on parallel processors that
involved coupling and reduction of initialization bias. Their work
did not use perfect sampling as in Propp and Wilson
\shortcite{Propp96} (described below), and did not yield samples drawn
perfectly from the steady-state model, but the coupling they used is the same as
the coupling used in Section \ref{sec:MM} below.
\subsection{Algorithms}
The {\em coupling from the past} (CFTP) algorithm was introduced by
Propp and Wilson \shortcite{Propp96}. It allows a simulation of a
Markov chain to be used to produce a simulated value whose
distribution is the limiting distribution of the chain. Since Markov
chain Monte Carlo algorithms are often used to produce approximate
inference about the steady state distribution, CFTP and related
algorithms have come to be known as {\em exact} or {\em perfect}
simulation.
The simplest description
of CFTP applies to an aperiodic discrete time Markov chain, though it can be
applied much more generally; indeed, our applications below are more
general. We suppose that the Markov chain $X_t$ may be simulated using
a recursive formulation
\begin{equation}
\label{eq:srs}
X_{t+1} = \phi(X_t, U_{t+1})
\end{equation}
where $U_{t+1}$ are i.i.d. from some known distribution.
The algorithm proceeds as follows:
\begin{enumerate}
\item Conceptually, generate the sequence $U_{0}, U_{-1}, \ldots$.
In
actual practice these values are generated as needed.
\item Search for a time $-T < 0$ such that paths started at every
possible starting point at that time, updated using (\ref{eq:srs})
repeatedly, all produce the same simulated value for $X_0$ (i.e. they
{\em coalesce}).
\item Output the common value $X_0$.
\end{enumerate}
Typically the search for $-T$ proceeds by trying some value $-T_0$,
then
trying $-2T_0$ if that fails, and continually doubling until success
is
achieved. The art of creating CFTP algorithms lies in constructing
$\phi(\cdot, \cdot)$ so that the search terminates. It is easy to
see
that if the probability of termination is greater than zero for some
finite $-T$, then the search terminates with probability 1: each new
segment further back in time is independent of the already-examined
tail
segment. If the state space is very large (e.g. a continuum), then
tricks such as monotonicity [\citeNP{Propp96}; see discussion below]
and gamma couplers
\cite{Murdoch98} are used to reduce the amount of computation to
manageable levels.
For continuous time jump processes, we need a more complicated
formalism, but the same general ideas apply. We need an easily
simulated point process $U_t$ (e.g., a constant rate Poisson process) to
drive $X_t$. For CFTP, we need to be able to simulate $U_t$ in a time
interval $[-T, 0]$, and then, conditional on these values, simulate it
in an earlier interval $[-2T, -T]$ (and so on recursively). The value
of $X_t$ at time 0 must be completely determined by the value of $X_{-T}
$ and the realization of $U_t$ on $[-T, 0]$. The general CFTP algorithm
then simulates $U_t$ on $[-T, 0]$ and doubles $T$ until every
possible value of $X_{-T}$ is updated to the same output $X_0$
(coalescence). In this general setting CFTP will not necessarily
terminate. For the cases where CFTP terminates with
probability 1, the value of $X_0$ is
determined by the realization of the $U_t$ process on some finite
interval $[-T, 0]$;
the same value would be drawn when started at any value at any earlier time.
Thus its distribution is the same as the limiting distribution of the chain.
Proving that CFTP terminates is typically straightforward in our settings. If realizations of
$U_t$ are independent in non-overlapping intervals (as with a Poisson
process), then CFTP will terminate provided the probability of
coalescence is positive. For more general conditions and a rigorous
proof of coupling in finite time see Theorem 8.2 of Chapter 10
in Thorisson \shortcite{Thorisson00}. We consider a more general case in Section \ref{sec:pareto}
below.
At this point, most readers newly exposed to CFTP will be wondering about
the complication of going backwards in time. Why not just run coupled chains
forwards until coalescence as in Andrad\'{o}ttir and Ott \shortcite{Andradottir95}?
The answer is that the state at the time of coalescence is typically not
a draw from the steady-state of the chain. CFTP
gives the same output no matter how large $T$ is chosen to be; thus
the output must be a draw from the limiting distribution.
Forward coupled chains converge to the limiting distribution, but the distribution
at any fixed time point, or fixed time after a coalescence, will generally
not be in steady-state: the event of coalescence may bias the distribution.
Nevertheless, we would like to avoid one of the difficulties with CFTP:
it requires the $U_t$
values
to be used repeatedly. When an attempt for coalescence fails and $T$
is
increased, it is important that the $U_t$ values already generated (i.e.
those that fall in $[-T, 0]$)
are
re-used, as this guarantees that the output value is the same
regardless
of the search strategy. In practice, this is usually done in one of
two ways.
Either the $U_t$ values can be stored as generated, or they can be
regenerated as needed, by storing random number seeds and re-using
them.
Both require an unbounded amount of storage availability, with the
needs of the former
growing exponentially faster than those of the latter.
The latter uses pseudo-random
number generators (PRNGs) in ways for which they were not likely
designed or tested:
whereas many PRNGs are extensively tested and found to be good
approximations to i.i.d. samples in the order generated, it is not
always clear that they will be as good when the values are
re-ordered. If hardware based random number generators are employed,
the storage strategy is the only one available.
Read-once CFTP (ROCFTP) addresses these concerns \cite{Wilson00}.
This
variation on the algorithm proceeds as follows. Again the
description is for aperiodic discrete time Markov chains.
\begin{enumerate}
\item Choose a block size $T$.
\item Simulate blocks of updates of $T$ steps forward in time until a block is found
in
which coalescence occurs, say block $n$, consisting of the $T$ steps between the
$T+1$ states at times $(n-1)T$ through $nT$. Set $X_{nT}$ (the ``special
path'') to the value at the end of this block.
\item Continue to simulate blocks forward in time, following paths
from
all states at the start of each block, but paying particular
attention
to the special path, until another coalescent block is
detected, say block $m>n$.
\item Output $X_{(m-1)T}$, the value of the special path at the {\em
start} of the next coalescent block.
\end{enumerate}
This algorithm requires that the $U_t$ values within non-overlapping
time blocks are
mutually independent, so the events that block $i$ coalesces,
which we denote by $C_i$, $i=1, \ldots$,
are i.i.d. events. ROCFTP can then be seen as equivalent to CFTP
with a relabelling of the time axis, where time $(m-1)T$ is
relabelled as time $0$. The output value depends only on $U_t$ values in
blocks $n, \ldots, m-1$, which are independent of the values
in blocks $m$ or greater; by the time-invariance of the structure, the relabelling
of the time axis has no effect on the distribution of the output value.
This is a discrete-time version of the PASTA property \cite{Wolff82}.
The distribution of a stationary process at the fixed time 0 is the
same as its distribution just before a randomly sampled time, provided
the sampling time is independent of the history of the process. The
reason ROCFTP requires two coalescence events instead of the one required
by CFTP is that the $C_i$ events play two roles: the second one determines
the sampling time, and the first guarantees that the value of the
special path will be known then.
Because the blocks are independent, ROCFTP will always terminate provided
$P(C_i)$ is positive.
Wilson \shortcite{Wilson00} shows that ROCFTP is most efficient when $T$
is chosen near the median forward coalescence time so that
$P(C_i)$ is around 1/2; in that case it is typically about
1/3 less time-efficient than CFTP, though it can be more efficient in
certain special cases. As discussed above, it generally makes more
efficient use of memory.
As we will see below, implementations of CFTP in which the $U_t$
values are not independent between blocks cannot easily be redone as
ROCFTP. However, in cases where independent blocks are
available, we would generally recommend ROCFTP over CFTP. It is simpler
and/or less resource-intensive to program.
\subsection{What Distribution to Simulate?}
When we say we are producing a draw from the steady-state
distribution of a queue or network model, we need to ask: which
steady-state distribution? As is well-known, the time-average
steady-state (i.e. the limit of the
distribution of the state at time $t$, as $t \rightarrow \infty$) may
be different from the distribution of the state at the $n^{\rm th}$
arrival, as $n \rightarrow \infty$. The descriptions of CFTP and ROCFTP
apply to time-average steady-states, as time 0 is a fixed time in a long
running (from the past) Markov chain.
In the examples in Sections \ref{sec:MM}, \ref{sec:queue} and
\ref{sec:network} below there is no distinction, because of the
Poisson arrival assumption. In Section \ref{sec:pareto} we show
how to simulate both steady states.
\section{Examples}
\subsection{Example: M/M/1/C Queue}
\label{sec:MM}
Consider a queue with a single server and a finite buffer, with total
capacity $C+1$. Customers
who arrive when the buffer is full are lost. Arrivals occur
in a Poisson process at rate $\lambda$; service is
at rate $\mu$ (i.e. there are exponential service times with mean $1/\mu$.)
The state of this queue is entirely represented by $X$, the number of customers
currently in the system. Transitions to a higher state occur
at rate $\lambda$, provided $X < C+1$; transitions
to a lower state occur at rate $\mu$ provided $X > 0$. Steady-state analysis
of this queue is straightforward; for example, see Gross and Harris
[\citeyearNP{Gross98} \S 2.4]. However, we will construct a perfect sampler
for it as a ``warm-up''.
Both CFTP or ROCFTP are straightforward to apply to this queue. There is only a finite
set of states; one could define random updates for each possible state,
and follow them all forward through time until they coalesced. However,
if $C$ is large, this will be impractical; there will be too many states
to follow, and it may take too long for them to coalesce.
One approach often used to reduce the amount of computation is to apply
a (partial) order to the state space, and to choose an update function
$\phi$ that preserves this order, i.e. $x \preceq y$ implies
$\phi(x,u) \preceq \phi(y,u)$
for all $u$, where ``$\preceq$'' represents the partial order relation.
This property is known as {\em
monotonicity} of the update function. Under monotonicity only maximal
and minimal elements need to be followed; all others are sandwiched between
them. The monotonicity approach works well here if we apply the same potential
updates to every state, because the natural integer ordering ``$\le$'' on
states is preserved under these updates. Thus we can determine coalescence
by following states $X=0$ and $X=C+1$ until those
two states coalesce.
One minor complication is that the $U_t$ process (here the arrival and
service events) must be generated independently of the current state of the
system, but when the queue is in state 0 no service events occur. However,
this is dealt with simply by generating service events uniformly across
time, even when the queue is empty, and having the queue
reject them if it is in state 0, just as it rejects arrivals when it is
in state $C+1$.
\begin{figure}
\center
\begin{minipage}{4in}
\begin{algorithm}
MMrocftp(T):\\
\>repeat \>\>\>\>\>\>\em Find first coalescence\\
\>\>$L \leftarrow 0$ \>\>\>\>\>\em Initialize\\
\>\>$X \leftarrow$ undefined \\
\>\>$R \leftarrow C+1$ \\
\>\>\mbox{MMblockupdate}(T,L,X,R)\\
\>until $L=R$\\
\>$X \leftarrow R$\\
\>repeat \>\>\>\>\>\>\em Find next coalescence\\
\>\>$L \leftarrow 0$\\
\>\>$R \leftarrow C+1$ \\
\>\>$Y \leftarrow X$ \>\>\>\>\>\em Save start value\\
\>\>\mbox{MMblockupdate}(T,L,X,R)\\
\>until $L=R$\\
\>return $Y$\\
\end{algorithm}
\end{minipage}
\caption{ROCFTP algorithm for M/M/1/C Queue. (Continued on next page) \label{fig:MMrocftp}}
\end{figure}
\addtocounter{figure}{-1}
\begin{figure}
\center
\begin{minipage}{4in}
\begin{algorithm}
MMblockupdate(T, L, X, R): \>\>\>\>\>\>\>\em Block update function\\
\>$ t \leftarrow 0$ \>\>\>\>\>\>\em time within block\\
\>repeat\\
\>\>$t \leftarrow t + \Exp(\lambda+\mu)$ \>\>\>\>\>\em Generate time of next event\\
\>\>if $t > T$ then\\
\>\>\>return modified $(L, X, R)$\\
\>\>else\\
\>\>\>$U \leftarrow \Unif(0,1)$ \>\>\>\>\em Used to determine event type\\
\>\>\>$L \leftarrow \mbox{MMupdate}(L,U)$\\
\>\>\>$X \leftarrow \mbox{MMupdate}(X,U)$\\
\>\>\>$R \leftarrow \mbox{MMupdate}(R,U)$\\
\> until return\\
\myrule\\
MMupdate(x,u): \>\>\>\>\>\>\>\em State update function\\
\>if $u < \lambda/(\lambda+\mu)$ then \\
\>\>$x \leftarrow \min(x+1, C+1)$ \>\>\>\>\>\em Arrival event \\
\>else\\
\>\>$x \leftarrow \max(x-1, 0)$ \>\>\>\>\>\em Service event \\
\>Return $x$\\
\end{algorithm}
\end{minipage}
\caption{(Continued from previous page) ROCFTP algorithm for M/M/1/C Queue.}
\end{figure}
With these choices, we implement ROCFTP with block
size $T$ as shown in Fig.~\ref{fig:MMrocftp}, where we use $L \le
X \le R$ to represent the state at time $t$ of the lower limit, steady-state
and upper limit chains respectively, $\Exp(\lambda + \mu)$ and $\Unif(0,1)$
are randomly generated values from the Exponential distribution with
rate $\lambda+\mu$ and the uniform distribution on $[0,1]$ respectively.
The arrival and service events are generated by creating unlabelled
events at the sum of the respective rates, and then randomly labelling
them as arrival or service events. The last event generated in {\em MMblockupdate}
will generally fall outside the block; by the memorylessness of Poisson arrivals,
it is safely discarded, and a new event
generated the next time this subroutine is called.
An implementation of CFTP would also be straightforward. There the main
implementation issue is how the events are generated and saved. A
simple pseudocode implementation that assumes they are stored and re-used is
given in Murdoch and Green \shortcite{Murdoch98}.
\subsection{Example: M/M/C/C Queue with Trunk Reservations}
\label{sec:queue}
In this section we generalize the M/M/1/C queue of Section~\ref{sec:MM}
to a priority queue with capacity $C$, and $p$ customer types.
Customer type $i=1,\ldots,p$ arrives in a Poisson process with rate
$\lambda_i$ and uses $b_i$ units of capacity when accepted for service.
However, there is a trunk reservation of $r_i$ units: an arriving customer
of type $i$ will only be accepted if there are currently at least $r_i+b_i$
units of spare capacity on the server. Customers who are not accepted are lost.
Service times are exponential with service rate $\mu_i$ for type $i$
customers. Because of the exponential
interarrival and service times, we can represent the state as the vector
$X=(X^{(1)},\ldots,X^{(p)})$, $X^{(i)}\in Z^{+}$, $\sum_iX^{(i)}b_i\le C$,
giving the counts of customers of each type currently being serviced.
With all $r_i=0$, an explicit analysis of this queue is straightforward
[\citeNP{Ross95}, Ch. 2]. However, positive trunk reservations complicate
the analysis significantly. For example, if $p=2$ and $r_10$.
The parameter $c_i>1$ is chosen to control the ``burstiness'' of the
arrival
process, with lower values more bursty; our simulations use $c_i=1.1$
for
all $i$.
The parameter $a_i$ is chosen to match the mean interarrival time
$a_i/(c_i-1)$
to the corresponding value $1/\lambda_i$ in the Section
\ref{sec:queue}
simulation.
We will label the average arrival rate as $\lambda_i$.
At steady-state, the time since the last type $i$ arrival will follow
the so-called
``current life'' distribution, with density function
$[1-F_i(x)]\lambda_i$.
For our Pareto arrivals, this distribution is
Pareto$(a_i,\ c_i-1)$.
We still assume exponential service times. Then in order that we
have a
Markov
process, it is sufficient to supplement the counts $X^{(i)}$ of
each type of customer currently in the system with the times
since the last arrival of each type of customer.
Since arrivals are independent of the service process, we can easily
simulate the arrival process at steady-state, and use this to drive
the
CFTP algorithm.
\begin{enumerate}
\item \label{step:initialization} Simulate $Y_i$, $i=1,\ldots,p$ from
the density
$(1-F_i(x))\lambda_i$. Set the most recent arrival time of customers
of type $i$ to $A^{(i)}_1 = -Y_i$ for $i=1,\ldots,p$.
\item Choose $-T < 0$.
\item \label{step:arrivals} Simulate arrival times $A^{(i)}_j$,
$j=2,\ldots$ of each customer type $i=1,\ldots,p$
backwards until $A^{(i)}_{j_i} < -T$ for each $i$. These are
simulated recursively by drawing the interarrival times from $F_i$.
Save these arrival times, and the current value of the random number
seed.
\item Generate a labelled service time following time $-T$ as in
Section \ref{sec:queue}.
\item Set the current hyper-rectangle of states to cover all possible
states.
\item Repeat the following loop forward until time 0 is reached.
\begin{enumerate}
\item Choose whichever of the $p$ arrival times or the labelled
service
time comes first. With Pareto interarrival times and exponential
service times ties will not occur, but in general
ties could arise. In that case ties should be broken as they would be
in the system being modelled, e.g. higher priority customers could be
treated as arriving first.
\item Process that event as in Section \ref{sec:queue} to update the
hyper-rectangle.
\item Replace an arrival time by choosing the next saved
$A^{(i)}_{j_i}$ value. If none remain, use an arbitrary positive
value instead. Replace a service time by generating a new labelled
service time using the appropriate exponential inter-service time
distribution.
\end{enumerate}
\item Check for coalescence of the hyper-rectangle to a single state.
If it has occurred, output the state at time 0. Otherwise, double
$T$, and repeat from step \ref{step:arrivals}.
\end{enumerate}
If we are interested in simulating the state at the time of a type
$i$ arrival, we simply replace the distribution used in step
\ref{step:initialization} for that type with $F_i$ (i.e. we assume
that there was a type $i$ arrival just after time 0).
This algorithm requires that large amounts of storage be devoted to
the $A^{(i)}_j$ values. An alternative implementation with fixed
storage would regenerate them as needed, by restoring the random number generator seed
and cycling through the values. However, this could be very
time-consuming.
Unfortunately, implementing ROCFTP is not straightforward here. The
arrivals within blocks of time of a fixed length are not independent
of other blocks. If there were only one customer type, the blocks
could be defined in terms of the number of arrivals of that type of
customer; then independence would be achieved. However, this is not
possible with $p>1$.
It might be possible to implement ROCFTP using the gamma coupling
techniques of Murdoch and Green \shortcite{Murdoch98}. These would
require expanding the state space to include the times since the last
arrival at the start of each block, and coupling all possible
conditional updates. However, these are difficult to implement, and
not always successful, in that an infinite number of paths might
result. It seems simpler to use the CFTP algorithm in
this case. Termination is guaranteed because there is a non-zero
probability of coalescence in any interval (e.g. realized
service times fast enough to empty even a full queue after the last
arrival).
Fig.~\ref{fig:pareto} shows the output from 3 runs of the CFTP
algorithm for this queue, illustrating both qualitative and quantitative
differences between this queue and the M/M/C/C queue illustrated in
Fig.~\ref{fig:mmccstationary}. On the left, the time-average steady-state
distribution has been sampled. The
queue spends about forty percent of the time completely empty, as shown
by the large concentration of points at $(0,0)$. Fifty percent of the
time there is only one type of customer in the queue. The centre and
right plots show that arriving customers see very different behaviour.
Arrivals almost never (only 0.3\% of the simulations) see an empty
system, and thirty percent of the time see customers of both types.
Refusal rates can also be calculated from the simulations by looking at
the proportion of simulations on the solid line for type 1 (i.e. 66\%),
and on or above the dashed line for type 2 (i.e. 72\%).
\begin{figure}
\begin{center}
\includegraphics[width=5in]{pareto.ps}
\end{center}
\caption{\label{fig:pareto} Simulations of the G/M/C/C queue with trunk reservations.
The plots show the three different steady-state distributions under bursty Pareto
arrivals.}
\end{figure}
\subsection{Example: G/G/C/C Queue}
When both arrivals and service times come from general distributions,
perfect sampling is more difficult. These chains are Markov only if the
arrival times of each individual in the system are included as part of the state
space; that makes the state space very large, and in general it appears
that gamma coupling techniques \cite{Murdoch98}
may be necessary. However, in the
special case where the service times are bounded, an approach similar to
the one above may be successful.
One simple perfect sampler for the case of bounded service times would
be to use CFTP as follows. Simulate arrival times backwards from time 0
until a gap in arrivals is observed that is longer than the upper bound
on service times. By the end of that gap, the system will be known to
be empty. Start a path at that known state and simulate forward in time
to time zero. The
difficulty with this approach is that such regeneration points may not
occur frequently enough; the search would have to go too far into the
past to be practical. Here we describe a CFTP scheme which detects
coalescence without requiring regeneration.
For simplicity, we assume in this section that there is only one type of
customer (i.e. $p=1$). The methods may be extended to more types.
We start by assuming there are $C$ servers and no buffer. Customers
arrive in a process with interarrival distribution $F(\cdot)$; those
who arrive when a server is free occupy it for a service time with
distribution function $G(\cdot)$, while those who arrive when all
servers are busy are lost. We make the same assumptions on $F$ and $G$
as in the previous section (i.e., finite
variance and non-arithmetic).
For the general discussion we
will assume without loss of generality that the service time
distribution is bounded above by 1; for concreteness we will
assume $G(\cdot)$ is Uniform$(0,1)$.
The boundedness of the service time distribution is useful because of
the following Markov-like property: the state at time $t$ is determined
entirely by events (arrivals, their acceptance or rejection, and their
service) in the time interval going back to time $t-1$. To implement
CFTP, we can use the
same technique as in the previous section to simulate the arrival times
backwards from time 0. Whether a particular arrival is accepted or not
depends on the number of customers in the system at the time of the
arrival, which we do not know: some unknown number of customers (the
``old customers") will be left from arrival times which have not been
simulated.
Furthermore, service completion times
for the old customers are not known, other than
the fact that they will occur sometime within the first time unit in our
simulation.
To handle this, we store the simulation in two parts. First, we store
the values of all arrivals in an ``arrivals list'' $0 > A_1 > A_2 >
\cdots$. As per usual in CFTP, we only simulate these values as needed.
For each simulated
arrival $A_j$, we store a virtual ``service completion time'' $s_j >
A_j$, i.e. the
time at which that individual would be serviced if they were accepted
into the system at time $A_j$. When doing CFTP, we generate enough
$A_j$'s to know all arrival events on the interval $[-T, 0]$.
The second part of the simulation is the storage of the possible current
states of
the system at an instant $t \in [-T, 0]$. Possible states are of the form $X_t =
(n,
\{j_i\}_{i=1}^m)$, where $\{0,\ldots,n\}$ is the range of possible
numbers of old
customers still in the
system, $m$ is the number of new customers currently in the system, and
$j_1, \ldots, j_m$ are the indices of the customers' arrival times from
the
the $A_j$ array. A state of the queue
changes
at the unknown times of service for old customers (when $n$ decreases by
one, reaching 0 by $-T+1$),
at the times of accepted arrivals (when $m$ increases by one, and the
new
service time is added to the service completion list), and at the times
$s_{j_i}$ in the list,
when $m$ decreases by one and the corresponding $j_i$ is dropped. However,
because the service times for old customers are unobserved, we cannot
update $n$ until $t=-T+1$, at which point it becomes 0.
The
algorithm is then described informally
as follows:
\begin{enumerate}
\item Simulate $Y$ from
the density
$(1-F(x))\lambda$. Set the most recent arrival time to $A_1 = -Y$.
\item Choose $-T < 0$.
\item Simulate arrival times $A_j$,
$j=2,\ldots$ backwards until $A_{j} < -T$. These are
simulated recursively by drawing the interarrival times from $F$.
Save these arrival times, and the current value of the random number
seed.
\item Simulate the corresponding service durations from $G$; add these
to the $A_j$ values to give $s_j$.
\item Set $t=-T$.
\item Set $n$ to $C$, set $m$
to
0, and initialize the list of $j_i$ to be empty.
\item Repeat the following loop forward until time 0 is reached.
\begin{enumerate}
\item Set the new $t$ to whichever of the arrival times $A_j$ or
service times
$s_{j_i}$
comes first.
\item If $t > -T+1$, set $n$ to \{0\}.
\item If we have a service event first and $j_i$ is in the list of
active arrivals, delete
$j_i$, and decrement $m$ by 1.
\item If the arrival comes first, then our action depends on $n$.
\begin{enumerate}
\item If $n+m \ge C$, then
it is possible that the arrival will be rejected, so keep the current
state as a possible state of the process. Otherwise, delete it.
\item If $m < C$, then
create a state which accepts the arrival (i.e. increments $m$, adds $j$
to the service time list, and, if necessary, reduces the range of $n$
so that $n+m \le C$ is guaranteed. If not, do nothing.
\end{enumerate}
Note that either one or two states may be output from the tests above.
\item Remove any duplicate states, or other redundant states. For
example, if two states have the same set of $s_j$ values but different
values of $n$, then it is only necessary to keep the larger
value.
\end{enumerate}
\item Check for coalescence at time 0 by simply counting how many unique
states are being followed. If only one, output it; if not double
$T$, and repeat from step \ref{step:arrivals}.
\end{enumerate}
\section{Concluding Remarks}
In the models that we considered in this paper, monotonicity was not always achievable,
nor were all of the models Markov. Nevertheless, CFTP or ROCFTP allowed
relatively straightforward simulation from the steady-state distribution
of the chain. What was needed was a simulation that was driven by
random inputs which could be simulated backward or forward in time.
Once we had such a simulation we used it to calculate steady-state properties.
For example, the simulations illustrated in Figs.~\ref{fig:mmccstationary}
and \ref{fig:pareto} confirm the expected
qualitative differences between Poisson arrivals and bursty
arrivals, but perfect simulation also allows us to obtain
quantitative measures. In practical applications, this information
could be used for efficient planning of telephone networks or
other queues.
We have provided a set of examples that give a flavor of how
perfect sampling may be applied to queues. A further useful
technique may be state space aggregation (called multistage
backward coupling in Meng \shortcite{Meng2000}). In conventional
simulation scenarios, one is often interested in relatively few
performance statistics, such as blocking probabilities. The
probability of being in a given subset of the state space could be
efficiently simulated by mapping the state vector to a Bernoulli
random variable and declaring ``coalescence" based on the values of
the Bernoulli sequences being followed. In carrying this out,
it is essential that the law of the original process be used
to carry out the updates; efficiency comes from the fact that
coalescence of the original process to a point is not needed, only
coalescence of the Bernoulli sequence. This technique should be
efficient even for high dimensional state spaces. A benefit of
perfect sampling is that one knows the exact distribution of the
target probability estimators.
One might ask whether the extra effort to produce a perfect sample is worthwhile.
We believe that the answer is in the affirmative. First, stochastic
simulation is an effective way to calculate properties of models that are
algebraically intractable, and even for those models where exact algebraic solutions
exist, is often more easily implemented and checked than the exact solution.
Perfect sampling improves
on other stochastic simulation methods by obviating questions of convergence:
the ``burn-in'' stage is not required, since results are guaranteed to be distributed
according to the steady-state distribution. It is true that
each random value generated takes a lot
of computation time, but Murdoch and Rosenthal \shortcite{Murdoch99} describe ways to
combine perfect simulation methods with standard stochastic simulation in order to
achieve the unbiasedness of the former and the efficiency of the latter. We believe
algorithms in this class should always be employed when they are available.
\section*{Acknowledgments}
This work was supported in part by NSERC grants to both authors. The authors are
grateful to the anonymous referees and editors for some very helpful comments.
\bibliography{coupling,reflist,queues}
\begin{received}
...
\end{received}
\end{document}