8.1 Types of Statistical Variables
A simulation experiment occurs when the modeler sets the input parameters to the model and executes the simulation. This causes events to occur and the simulation model to evolve over time. During the execution of the simulation, the behavior of the system is observed and various statistical quantities computed. When the simulation reaches its termination point, the statistical quantities are summarized in the form of output reports.
A simulation experiment may be for a single replication of the model or may have multiple replications. A replication is the generation of one sample path which represents the evolution of the system from its initial conditions to its ending conditions. If you have multiple replications within an experiment, each replication represents a different sample path, starting from the same initial conditions and being driven by the same input parameter settings. Because the randomness within the simulation can be controlled, the underlying random numbers used within each replication of the simulation can be made to be independent. Thus, as the name implies, each replication is an independently generated “repeat” of the simulation.
Within a single sample path (replication), the statistical behavior of the model can be observed.
Within replication statistics are collected based on the observation of the sample path and include observations on entities, state changes, etc. that occur during a sample path execution. The observations used to form within replication statistics are not likely to be independent and identically distributed. Since across replication statistics are formed from the final values of within replication statistics, one observation per replication is available. Since each replication is considered independent, the observations that form the sample for across replication statistics are likely to be independent and identically distributed. The statistical properties of within and across replication statistics are inherently different and require different methods of analysis. Of the two, within replication statistics are the more challenging from a statistical standpoint.
For within replication statistical collection there are two primary types of observations: tally and time-persistent. Tally data represent a sequence of equally weighted data values that do not persist over time. This type of data is associated with the duration or interval of time that an object is in a particular state or how often the object is in a particular state. As such it is observed by marking (tallying) the time that the object enters the state and the time that the object exits the state. Once the state change takes place, the observation is over (it is gone, it does not persist, etc.). If we did not observe the state change, then we would have missed the observation. The time spent in queue, the count of the number of customers served, whether or not a particular customer waited longer than 10 minutes are all examples of tally data.
Time-persistent observations represent a sequence of values that persist over some specified amount of time with that value being weighted by the amount of time over which the value persists. These observations are directly associated with the values of the state variables within the model. The value of a time-persistent observation persists in time. For example, the number of customers in the system is a common state variable. If we want to collect the average number of customers in the system over time, then this will be a time-persistent statistic. While the value of the number of customers in the system changes at discrete points in time, it holds (or persists with) that value over a duration of time. This is why it is called a time-persistent variable.
Figure 8.1 illustrates a single sample path for the number of customers in a queue over a period of time. From this sample path, events and subsequent statistical quantities can be observed.
Figure 8.1: Sample Path for Tally and Time-Persistent Data
Let \(A_i \; i = 1 \ldots n\) represent the time that the \(i^{th}\) customer enters the queue
Let \(D_i \; i = 1 \ldots n\) represent the time that the \(i^{th}\) customer exits the queue
Let \(W_i = D_i - A_i \; i = 1 \ldots n\) represent the time that the \(i^{th}\) customer spends in the queue
Thus, \(W_i \; i = 1 \ldots n\) represents the sequence of wait times for the queue, each of which can be individually observed and tallied. This is tally type data because the customer enters a state (the queued state) at time \(A_i\) and exits the state at time \(D_i\). When the customer exits the queue at time \(D_i\), the waiting time in queue, \(W_i = D_i - A_i\) can be observed or tallied. \(W_i\) is only observable at the instant \(D_i\). This makes \(W_i\) tally based data and, once observed, its value never changes again with respect to time. Tally data is most often associated with an entity that is moving through states that are implied by the simulation model. An observation becomes available each time the entity enters and subsequently exits the state.
With tally data it is natural to compute the sample average as a measure of the central tendency of the data. Assume that you can observe \(n\) customers entering and existing the queue, then the average waiting time across the \(n\) customers is given by:
\[\bar{W}(n) = \dfrac{1}{n} \sum_{i=1}^{n} W_{i}\]
Many other statistical quantities, such as the minimum, maximum, and sample variance, etc. can also be computed from these observations. Unfortunately, within replication data is often (if not always) correlated with respect to time. In other words, within replication observations like, \(W_i \, i = 1 \ldots n\), are not statistically independent. In fact, they are likely to also not be identically distributed. Both of these issues will be discussed when the analysis of infinite horizon or steady state simulation models is presented.
The other type of statistical variable encountered within a replication is based on time-persistent observations. Let \(q(t), t_0 < t \leq t_n\) be the number of customers in the queue at time \(t\). Note that \(q(t) \in \lbrace 0,1,2,\ldots\rbrace\). As illustrated in Figure 8.1, \(q(t)\) is a function of time (a step function in this particular case). That is, for a given (realized) sample path, \(q(t)\) is a function that returns the number of customers in the queue at time \(t\).
The mean value theorem of calculus for integrals states that given a function, \(f(\cdot)\), continuous on an interval \([a, b]\), there exists a constant, \(c\), such that
\[\int_a^b f(x) \mathrm{d}x = f(c)(b-a)\]
The value, \(f(c)\), is called the mean value of the function. A similar function can be defined for \(q(t)\). In simulation, this function is called the time-average:
\[\bar{L}_q(n) = \frac{1}{t_n - t_0} \int_{t_0}^{t_n} q(t) \mathrm{d}t\]
This function represents the average with respect to time of the given state variable. This type of statistical variable is called time-persistent because \(q(t)\) is a function of time (i.e. it persists over time).
In the particular case where \(q(t)\) represents the number of customers in the queue, \(q(t)\) will take on constant values during intervals of time corresponding to when the queue has a certain number of customers. Let \(q(t)\) = \(q_k\) for \(t_{k-1} \leq t \leq t_k\) and define \(v_k = t_k - t_{k-1}\) then the time-average number in queue can be rewritten as follows:
\[ \begin{aligned} \bar{L}_q(n) &= \frac{1}{t_n - t_0} \int_{t_0}^{t_n} q(t) \mathrm{d}t\\ & = \frac{1}{t_n - t_0} \sum_{k=1}^{n} q_k (t_k - t_{k-1})\\ & = \frac{\sum_{k=1}^{n} q_k v_k}{t_n - t_0} \\ & = \frac{\sum_{k=1}^{n} q_k v_k}{\sum_{k=1}^{n} v_k} \end{aligned} \]
Note that \(q_k (t_k - t_{k-1})\) is the area under \(q(t)\) over the interval \(t_{k-1} \leq t \leq t_k\) and
\[t_n - t_0 = \sum_{k=1}^{n} v_k = (t_1 - t_0) + (t_2 - t_1) + \cdots (t_{n-1} - t_{n-2}) + (t_n - t_{n-1})\]
is the total time over which the variable is observed. Thus, the time average is simply the area under the curve divided by the amount of time over which the curve is observed. From this equation, it should be noted that each value of \(q_k\) is weighted by the length of time that the variable has the value. This is why the time average is often called the time-weighted average. If \(v_k = 1\), then the time average is the same as the sample average.
With time-persistent data, you often want to estimate the percentage of time that the variable takes on a particular value. Let \(T_i\) denote the total time during \(t_0 < t \leq t_n\) that the queue had \(q(t) = i\) customers. To compute \(T_i\), you sum all the rectangles corresponding to \(q(t) = i\), in the sample path. Because \(q(t) \in \lbrace 0,1,2,\ldots\rbrace\) there are an infinite number of possible value for \(q(t)\) in this example; however, within a finite sample path you can only observe a finite number of the possible values. The ratio of \(T_i\) to \(T = t_n - t_0\) can be used to estimate the percentage of time the queue had \(i\) customers. That is, define \(\hat{p}_i = T_i/T\) as an estimate of the proportion of time that the queue had \(i\) customers during the interval of observation.
Let’s look at an example. Consider Figure 8.1, which shows the change in queue length over a simulated period of 25 time units.
Compute the time average number in queue over the interval of time from 0 to 25.
Compute the percentage of time that the queue had \(\lbrace 0, 1, 2, 3, 4\rbrace\) customers
Since the queue length is a time-persistent variable, the time average queue length can be computed as:
\[ \begin{aligned} \bar{L}_q & = \frac{0(2-0) + 1(7-2) + 2(10-7) + 3(13-10) + 2(15-13) + 1(16-15)}{25}\\ & + \frac{4(17-16) + 3(18-17) + 2(21-18) + 1(22-21) + 0(25-22)}{25} \\ & = \frac{39}{25} = 1.56 \end{aligned} \]
To estimate the percentage of time that the queue had \(\lbrace 0, 1, 2, 3, 4\rbrace\) customers, the values of \(v_k = t_k - t_{k-1}\) need to be summed for whenever \(q(t) \in \lbrace 0, 1, 2, 3, 4 \rbrace\). This results in the following:
\[\hat{p}_0 = \dfrac{T_0}{T} = \dfrac{(2-0) + (25-22)}{25} = \dfrac{2 + 3}{25} = \dfrac{5}{25} = 0.2\]
\[\hat{p}_1 = \dfrac{T_1}{T} = \dfrac{(7-2) + (16-15) + (22-21)}{25} = \dfrac{5 + 1 + 1}{25} = \dfrac{7}{25} = 0.28\]
\[\hat{p}_2 = \dfrac{T_2}{T} = \dfrac{3 + 2 + 3)}{25} = \dfrac{8}{25} = 0.28\]
\[\hat{p}_3 = \dfrac{T_3}{T} = \dfrac{3 + 1}{25} = \dfrac{4}{25} = 0.16\]
\[\hat{p}_4 = \dfrac{T_4}{T} = \dfrac{1}{25} = \dfrac{1}{25} = 0.04\]
Notice that the sum of the \(\hat{p}_i\) adds to one. To compute the average waiting time in the queue, use the supplied values for each waiting time.
\[\bar{W}(8) = \frac{\sum_{i=1}^n W_i}{n} =\dfrac{0 + 11 + 8 + 7 + 2 + 5 + 6 + 0}{8} = \frac{39}{8} = 4.875\]
Notice that there were two customers, one at time 1.0 and another at time 23.0 that had waiting times of zero. The state graph did not move up or down at those times. Each unit increment in the queue length is equivalent to a new customer entering (and staying in) the queue. On the other hand, each unit decrement of the queue length signifies a departure of a customer from the queue. If you assume a first-in, first-out (FIFO) queue discipline, the waiting times of the six customers that entered the queue (and had to wait) are shown in the figure.
Now that we understand the type of data that occurs within a replication, we need to develop an understanding for the types of simulation situations that require specialized statistical analysis. The next section introduces this important topic.