Proceedings of the 1996 Winter Simulation Conference

ed. J. M. Charnes, D. J. Morrice, D. T. Brunner, and J. J. S~vain

DISCRETE-EVENT SIMULATION INPUT PROCESS MODELING

Lawrence M. Leemis

Department of Mathematics
College of William & Mary
Williamsburg, VA 23187-8795, U.S.A.

ABSTRACT much broader treatll1ent of input 1110deling than pre-

sented here (e.g., Law and Kelton 1991). These texts
General guidelines for selecting probabilistic input include more specific inforll1ation on statistical tests
models as part of a discrete-event sin1ulation study for independence, graphical 111ethods for model selec-
are presented. Two short examples illustrating input tion, parameter estimation techniques, and goodness-
modeling decisions are also presented, as opposed to of-fit tests. Advanced input lnodeling is considered by
a complete treatment of the subject. Nelson et al. (1995).
An input model can be specified in a variety of
1 INTRODUCTION \vays, such as a cUll1ulative distribution function, haz-
ard function, intensity function, or a variate-genera-
Discrete-event sill1ulation models typically have stoch- tion algorithm. An input model characterizes each of
astic cOll1ponents that mimic the probabilistic nature the stochastic elell1ents of a discrete-event simulation.
of the system under consideration. Successful input Figure 1 contains a taxonoll1y whose purpose is to
modeling requires a close 111atch between the input illustrate the scope of potential input 1110dels that are
model and the true underlying probabilistic mecha- available to simulation analysts. There is certainly no
nism associated with the system. The general ques- uniqueness in the branching structure of the taxon-
tion considered here is how to model an element (e.g .. omy. The branches under stochastic processes, for ex-
arrival process, service times) in a discrete-event Si1l1- ample, could have been state followed by time, rather
ulation given a data set collected on the element of than tim.e followed by state, as presented.
interest. Exalnples of specific models that could be placed
Since time and space for this tutorial is limited, on the branches of the taxonomy appear at the far
the following sill1plifying assull1ptions have been ll1ade. right of the diagram. Mixed, univariate, time-in-
dependent input models have empirical/trace-driven
• The modeler has access to a reliable source of
given as an possible model. All of the branches in-
randoll1 numbers. Most introductory sill1ula-
clude this particular model. A trace-driven input
tion textbooks (e.g., Law and Kelton 1991) con- model simply generates a process that is identical
sider random number generation algorithn1s. to the collected data val ues so as not to rely on a
• An algorithm is available for converting these parametric model. A sill1ple example is a sequence
random nUll1bers to randoll1 variates associated of arrival tin1es collected over a 24-hour tillle period.
with the input model to drive the simulation The trace-driven input model for the arrival process is
(Devroye 1986). generated by having arrivals occur at the same times
as the observed values.
• Data is available on the aspect of the simulation The upper half of the taxonomy contains models
of interest. For examples of input modeling in that are independent of time. These lnodels could
the absence of data, see Schmeiser and Deutsch have been called !v! onte Carlo models. Models are
(1977) or Law, McComas, and Vincent (1994). classified by whether there is one or several variables
of interest, and whether the distribution of these ran-
With these assumptions lill1iting the scope of this dom variables is discrete. continuous, or contains both
tutorial, the focus turns to selecting the appropriate continuous and discrete elen1ents. Exall1ples of uni-
probabilistic models for the random cOll1ponents in a variate discrete models include the binoillial distribu-
simulation model. Many simulation textbooks have a tion and a degenerate distribution with all of its
40 Leenlis

Discrete Binomial( n, p)
Degenerate ( c )

Normal (Il, 0- 2 )
Univariate Continuous Exponential( A)
Bezier curve

Mixed Empirical / Trace- driven

Time-independent
models

Discrete
Independent bino mial( n, p)

Multi variate Continuous

Normal( Il, 2:)

Mixed
Bivariate exponent

Input Models

Stationary
Markov chain
Discrete-state

Nonstationary

Discrete- ti me

Stationary
ARMA(p, q)
Continuous-state

Nonstationary
ARIMA(p, d, q)
Stochastic Processes

Stationary Poisson process( A)

Renewal process
Discrete-state Semi-Markov chain

Nonstationary Nonhomogeneous Poisson

process
Continuous-time

Stationary
Mar kov process
Continuous-state

Nonstationary

Figure 1: .:\ Taxonol11Y for Input Models

 Discrete Event Sinlulation Input Process l\Iodeling .Jl

mass at one value. Examples of continuous distribu- ulation of a queuing systenl. The serVIce times In
tions include the normal distribution, an exponential seconds are
distribution with a random parameter A (see, for ex-
ample, Martz and Waller 1982), and Bezier curves 105.84 28.92 98.64 .55.56 128.04 45.60
(Flanigan-Wagner and Wilson 1993). Bezier curves 67.80 105.12 48.48 51.84 173.40 51.96
offer a unique combination of the paranletric and non- .54.12 68.64 93.12 68.88 84.12 68.64
parametric approaches. An initial distribution is fit- 41.52 127.92 42.12 17.88 33.00.
ted to the data set, then the modeler decides whether [Although these service tinles COIlle from the life test-
differences between the empirical and fitted nl0dels ing literature (Lieblein and Zelen 1956), the sanle
represent sampling variability (chance variation) or principles apply to both input nl0deling and survival
an aspect of the distribution that should be included analysis.]
in the input model. The first step is to assess \vhether the observations
Examples of k-variable multivariate input 1110d- are independent and identically distributed (iid). The
els (see Johnson 1987) include a sequence of k in- data must be given in the order collected for inde-
dependent binomial random variables, a nlultivari- pendence to be assessed. Situations where the iid
ate normal distribution with Inean J1 and variance- assumption \\Tould not be valid include:
covariance lllatrix L and a bivariate exponential dis-
tribution (Barlow and Proschan 1981). • A new teller has been hired at a bank and the 23
The lower half of the taxonomy contains stochas- service tinles represent a task that has a steep
tic process models. These models are often used to learning curve. The expected service time is
solve problems at the system level, in addition to likely to decrease as the new teller learns how
serving as input models for simulations with stochas- to perform the task more efficiently.
tic elements. Models are classified by how tilne is
measured (discrete/continuous), the state space (dis- • The service times represent 23 conlpletion times
crete/ continuous) and whether the 1110del is station- of a physically demanding task during an 8-hour
ary in time. For Markov models, the discrete-state/ shift. If fatigue is a significant factor, the ex-
continuous-state branch typically determines whether pected time to cOlllplete the task is likely to
the model will be called a "chain" or a "process" , and increase with time.
the stationary/nonstationary branch typically deter- If a sinlple linear regression of the observation num-
mines whether the model will be preceded with the bers regressed against the service times shows a signif-
term "homogeneous" or "nonhomogeneous". Exanl- icant nonzero slope, then the iid assumption is prob-
pIes of discrete-time stochastic processes include ho- ably not appropriate.
mogeneous, discrete-time Markov chains (Ross 1993) Assume that there is a suspicion that a learning
and ARIMA time series models (Box and Jenkins curve is present. An appropriate hypothesis test is
1976). Since point processes are counting processes,
they have been placed on the continuous-time, dis-
crete-space branch. Although the Poisson, renewal
and nonhomogeneous Poisson processes are all pure HI : PI < 0
birth processes, more general point processes, such
associated with the linear model (Neter, Wasserman,
as one to model the number of custonlers in a queue,
and Kutner 1989)
can be placed on one of the continuous time, discrete-
space branches.

2 EXAMPLES where X is the observation nunlber, }l" is the service

time, 130 is the intercept, /31 is the slope, and ( is an
Two simple examples illustrate the types of decisions error term. Figure 2 shows a plot of the (Xi, Yi) pairs
that often arise in input modeling. The first exaln- for i = 1, 2, ... , 23, along wi th the estimated regres-
pIe determines an input model for service times and sion line. The p-value associated with the hypothesis
the second example determines an input model for an test is 0.14, which is not enough evidence to conclude
arrival process. that there is a statistically significant learning curve
present.
2.1 Service Time Model There are a number of other graphical and statis-
tical methods for assessing independence. These in-
Consider a data set of n =
23 service times collected clude analysis of the sample autocorrelation function
to determine an input model in a discrete-event sim- associated with the observations and a scatterplot of
 42 Leemis

Service
Time

150

100

50
..

0 Observation
'-----~-----.----.-----y---- Number
10 15 20

Figure 2: Service Time Vs. Observation Number Figure 3: Histogram of Service Times

adjacent observations. For this particular example, 68.64 seconds, and the observation in the far right-
assume that we are satisfied that the observations hand tail of the distribution, 173.40 seconds, tend to
are truly iid in order to perform a classical statistical indicate that a parametric analysis is more appropri-
analysis. ate. Since the input model is for service times, the
The next step in the analysis of this data set in- accurate modeling of the right-hand tail of the dis-
cludes plotting a histogram and calculating the values tribution is critical. These long service times signifi-
of some sample statistics. A histogram of the obser- cantly impact queuing statistics. For this particular
vations is shown in Figure 3. Although the data set data set, a parametric approach is chosen.
is small, a skewed bell-shaped pattern is apparent. There are dozens of choices for a univariate para-
The largest observation lies in the far right-hand tail metric model for the service times. These include gen-
of the distribution, so care must be taken to assure eral families of scalar distributions, modified scalar
that it is representative of the population. The sam- distributions and commonly-used parametric distri-
ple mean, standard deviation, coefficient of variation, butions (see Schmeiser 1990). Since the data is drawn
and skewness are from a continuous population and the support of the
s distribution is positive, a time-independent, univari-
x = 72.22 s = 37.49 - = 0.52
x ate, continuous input model is chosen. The shape
of the histogram indicates that the gamma, inverse

;?=
I n ( Xi -
-s-
_)3 = 0.88.
X Gaussian, log logistic, log normal, and Weibull dis-
tributions (Lawless 1982) are good candidates. The
t=l
Weibull distribution is analyzed in detail here. Simi-
Examples of the interpretations of these sample statis- lar approaches apply to the other distributions.
tics are: Parameter estimates for the Weibull distribution
can be found by least squares, the method of mo-
• A coefficient of variation s / x close to 1, along
ments, and maximum likelihood. Due to desirable
with the appropriate histogram shape, indicates
statistical properties, maximum likelihood is empha-
that the exponential distribution is a potential
sized here. The Weibull distribution has probability
input model.
density function
• A sample skewness close to 0 indicates that a
symmetric distribution (e.g., a normal distribu-
tion) is a potential input model.
where A is a positive scale parameter and K is a posi-
The next decision that needs to be made is whether tive shape parameter. Let Xl, X2, .•. ,X n be the data
a parametric or nonparametric input model should be values. The likelihood function is
used. One simple nonparametric model would repeat-
edly select one of the service times with probability
1/23. The small size of the data set, the tied value,
 Discrete Event Sinlulation Input Process Alodeling

The log likelihood function is F(O

1.0
n

log L('x, K) == n log K + n:n log'x + (n: - 1) L log Xi

i=1 0.8
n
-,x~ ~ X~
LJ z' 0.6
i=1

The 2 x 1 score vector has elements

0.4

alog L('x, K)
o,x 0.2

and
0.0

50 100 150

Figure 4: Empirical and Fitted Cunlulative Distribu-

When these equations are equated to zero, the simul-
tion Functions for the Service Times
taneous equations have no closed-form solution for ~
and K:
evaluated at the maxinlum likelihood estimators is
logL('x,~) = -113.691. Figure 4 shows the empiri-
cal cUDlulative distribution function (a step function
n n
with each step height l/n) along with the Weibull fit
~ + nlog'x + Llogxi - L(Axi)" logAxi = o. to the data.
K i=l i=l
The observed inforl11ation matrix is
To reduce the problem to a single unknown, the first
equation can be solved for ,x in terms of K yielding O(~ h:~'.) = [681,875000 875 ]
1
10.4 '

revealing a positive correlation between the elements

of the score vector. lJsing the fact that the like-
Law and Kelton (1991, p. 334) give an initial esti- lihood ratio statistic, 2[log L(~,~) - log L('x, K)]' is
mate for K that can be used in Newton's lnethod to asymptotically ,,2
with 2 degrees of freedom and that
numerically solve for the maximum likelihood esti- \~,O.05 = 5.99, a 95% confidence region for the pa-
mators. Qiao and Tsokos (1994) consider nUlllerical rameters is all ,x and K satisfying
problems with Newton's method and give an alterna-
2[-113.691-log L('x, K)] < 5.99.
ti ve algorithm for calculating the nlaxilllum likelihood
estimators. The 9591cl confidence region is shown in Figure 5. The
The score vector has a mean of 0 and a variance- line K = 1 is not interior to the region, indicating
covariance matrix I('x, n:) given by the 2 x 2 Fisher that the exponential distribution is not an appropri-
information nlatrix ate model for this particular data set.
As further proof that K: is significantly different
from 1, the standard errors of the distribution of the
parameter estimators can be computed by using the
inverse of the observed inforDlation 111atrix
The observed information matrix
0- 1 (~ ~.) = [0.00000165 -0.000139]
,K: -0.000139 0.108 .

This is the asymptotic variance-covariance 111atrix for

the parameter estimators ~ and K. The standard er-
can be used to estimate I (,x, K:). rors of the parameter estimators are the square roots
For the 23 service times, the fitted Weibull dis- of the diagonal elements
tribution has maximunl likelihood estimators ~ =
0.0122 and K- = 2.10. The log likelihood function 0->. = 0.00128
 44 Leemis

/I

F
0.020
1.0

0.015 0.8

0.6
O.OlO

0.4

0.005
0.2

0.0 0.0
L-.,----.-------.--------.-----,----K F
o 0.0 0.2 0.4 0.6 0.8 1.0

Figure 5: 95% Confidence Region Based on the Like- Figure 6: A P-P Plot for the Service Times
lihood Ratio Statistic

package also performs a goodness-of-fit test such as

Thus an asymptotic 95% confidence interval for f\, is the Kolmogorov-Smirnovor chi-square test, the dis-
tribution that best fits the data set can quickly be
2.10 - (1.96)(0.329) < f\, < 2.10 + (1.96)(0.329) determined.
or p-p and Q-Q plots can also be used to assess
1.46 < Ii < 2.74, n10del adequacy. A P-P plot, for example, is a plot
of the fitted cumulative distribution function at the
since ZO.025 == 1.96. Since this confidence interval does ith order statistic X(i), i.e., F(X(i))' versus the ad-
not contain 1, the inclusion of the Weibull shape pa- justed empirical cumulative distribution function, i.e.
r ameter Ii is justified. F(X(i») == i-~.5, for i == 1,2, .. . ,n. A plot where the
At this point, model adequacy should be assessed. points fall close to a line indicates a good fit. For
Since the chi-square goodness-of-fit test suffers from the 23 service times, a P-P plot for the Weibull fit is
arbitrary interval limits and should not be applied to shown in Figure 6, along with a line connecting (0, 0)
small data sets, the I(011l10gorov-SIl1irnov, Cramer- and (1, 1). P-P plots should be constructed for all
von Mises, or Anderson-Darling goodness-of-fit tests competing models.
(Lawless 1982) are appropriate here. The I(oln10gor-
ov-Smirnov test statistic, for example, for this data
2.2 Arrival Process Model
set with a Weibull fit is 0.152, which 111easures the
maximum difference between the empirical and fitted Arrival times to a lunch wagon between 10:00 AM and
cumulative distribution functions. This test statistic 2:30 PM are collected on three days. The realizations
corresponds to a p-value of approximately 0.1,5 (Law were generated from a hypothetical arrival process
and Kelton 1991, page 391), so the Weibull distri- given by Klein and Roberts (1984). A total of n ==
bution provides a reasonable 1110del for these service 1.50 arrival times were observed, including n1 == 56,
times. The l\olmogorov-Smirnov test statistic values n2 == 42 and n3 == 52 on the k == 3 days. Defining
for several potential 1110dels are shown below. (0,4.5] be the time interval of interest (in hours) the
three realizations are
Model Test statistic
Exponential 0.301 0.2152 0.3494 0.3943 4.175 4.248,
Weibull 0.1,52
Gamll1a 0.123 0.3927 0.6211 0.7504 4.044 4.374,
Inverse Gaussian 0.099
Log norn1al 0.090 and

Many of the discrete-event sin1ulation packages 0.4499 0.5495 0.6921 3.643 4.357.
exhibited at the ~Vinter Simulation Conference have
the capability of determining 111axin1UITI likelihood es- One preliminary statistical issue concerning this
tin1ators for several paran1etric distributions. If the data is whether the three days represent processes
 Discrete Event Simulation Input Process l\Jodeling -J5

drawn from the same population. External factors A(I)

60
such as the weather, day of the week, advertisement,
,.'
and workload should be kept fixed. For this partic-
50
ular example, these factors have been fixed and the I
I

three processes are representative of the population

40 r-
of arrival processes to the lunch wagon. I "
I

The input model for the process comes from the I

I

lower branch (stochastic processes) of the taxonomy 30 ,

in Figure 1. Furthermore, the arrival times consti- I

tute realizations of a continuous-time, discrete-state 20 ,

I
I

stochastic process, so the remaining question con- "

I
cerns whether or not the process is stationary. 10 ,- -.,

If the process proves to be stationary, the tech- ,,' ,--'

niques from the previous example, such as drawing
a histogram, and choosing a parametric or nonpara-
o
metric model for the interarrival times are appropri-
ate. This results in a Poisson or renewal process.
On the other hand, if the process is nonstationary, a Figure 7: Point and 95% Confidence Interval Estinla-
nonhomogeneous Poisson process might be an input tors for the Cumulative Intensity Function
appropriate model. A nonhomogeneous Poisson pro-
cess is governed by an intensity function A(t) which set since the intensity function can only increase, de-
gives an arrival rate [e.g., A(2) == 10 means that the crease or remain constant, and can not model an in-
arrival rate is 10 customers per hour at time 2] that tensity function that increases, then decreases. Since
can vary with time. the intensity function is analogous to the hazard func-
Figure 7 contains a plot of the empirical cumula- tion for tinle-independent nlodels, an appropriate 2-
tive intensity function estimator suggested by Leemis parameter distribution to consider would be one with
(1991) for the three realizations. The solid line de- a hazard function that increases initially, then de-
notes the point estimator for the cumulative inten- creases. A log-logistic process, for example, with in-
sity function A(t) == J~ A(r)dr and the dashed lines tensity function
denote 95% confidence intervals. The cumulative in-
tensity function estimator at time 4.5 is 150/3 == 50, AK( At r~-1
A(t) - - - - t > 0,
the point estimator for the expected number of arriv- - 1 + (At)~
ing customers per day. If A( t) is linear, a stationary
for A > 0 and K > 0, would certainly be an im-
model is appropriate. Since people are more likely
proved choice. A more general EPTF (exponential-
to arrive to the lunch wagon between 12:00 (t == 2)
polynomial-trigonometric function) model is given by
and 1 :00 (t == 3) than at other times and the cumu-
Lee, Wilson and Crawford (1991) with intensity func-
lative intensity function estimator has an S-shape, a
tion
nonstationary model is indicated. More specifically,
a nonhomogeneous Poisson process will be used to
nlodel the arrival process. t > O.
The next question to be determined is whether a
parametric or nonparametric model should be chosen The trigonometric function is capable of modeling the
for the process. Figure 7 indicates that the inten- intensity function that increases, then decreases.
sity function increases initially, remains fairly con- In all of the parametric nlodels, the likelihood
stant during the noon hour, then decreases. This function for the vector of unknown parameters B ==
may be difficult to model parametri~ally, so a non- (B 1 , ()2, ... , Bp ) from a single realization on (0, c] is
parametric approach, possibly using A(t) in Figure 7
might be appropriate.
There are many potential paranletric models for
nonstationary arrival processes. The power law, or
L(()) = [g'\(d] t ex p [- 1 c ,\ ( t )dt] .

Weibull process has intensity function Maximum likelihood estimators can be determined
by maximizing L( B) or its logarithm with respect to
t > 0, all unknown parameters. Confidence intervals for the
where A and K are positive parameters. This pop- unknown parameters can be found in a similar man-
ular model would not be appropriate for this data ner to the service time exanlple.
46 LeeIllis

ACKNOWLEDGMENTS 199,5. Input modeling when simple models fail. In

Proceedings of the 1995 Winter Simulation Con-
The author thanks Steve Tretheway for his help in ference, ed. C. Alexopoulos, K. Kang, W. R. Lileg-
developing Figure 1 while at The lTniversity of Ok- don, and D. Goldsman, 93-100. Institute of Elec-
lahoma. Thanks also goes to Steve Park for reading trical and Electronics Engineers, Piscataway, New
a draft of this tutorial. The figures have been done Jersey.
using the Splus and xfig programs. Neter, J., W. Wasserman, and M. H. Kutner. 1989.
.4pplied Linear Regression Models. 2d ed. Boston:
REFERENCES Ir\vin.
Ross, S. 11. 1993. Introduction to probability models.
Barlow, R. E., and F. Proschan. 1981. Statistical the- 5th ed. Boston: Academic Press.
ory of reliability and life testing: probability mod- Schmeiser, B. 1990. Simulation experiments. In Hand-
els. Silver Springs, Maryland: To begin with. books in OR fj MS, ed. D. P. Heyman and M. J.
Box, G., and G. Jenkins. 1976. Time series analy- Sobel, 296-330. New York: Elsevier Science Pub-
sis: forecasting and control. Oakland, C;alifornia: lishers.
Holden-Day. Schmeiser, B., and S. J. Deutsch. 1977. A versatile
Devroye, L. 1986. Non-uniform random variate gen- four-parameter family of probability distributions,
eration. New York: Springer-Verlag. suitable for simulation. AIlE Transactions 2:176-
Flanigan-Wagner, M., and J. R. Wilson. 1993. lTs- 182.
ing univariate Bezier distributions to 1110del sin1u- Qiao, H., and C. P. Tsokos. 1994. Parameter es-
lation input processes. In Proceedings of the 1993 timation of the Weibull probability distribution.
Winter Simulation Conference, ed. G. W. Evans, !vIathematics and Computers in Simulation 37:47-
M. Mollaghasemi, E. C. Russell, and vV. E. Biles, 55.
365-373. Institute of Electrical and Electronics
Engineers, Piscataway, Ne\v Jersey. AUTHOR BIOGRAPHY
Johnson, M. E. 1987. Multivariate statistical sznl,ula-
tion. New York: John Wiley & Sons. LAWRENCE M. LEEMIS is a professor in the
Klein, R. W., and S. D. Roberts. 1984. A time- Mathe111atics Depart111ent at the College of William
varying Poisson arrival process generator. S'im,u- & Mary. He received his BS and MS degrees in Math-
lation 43:193-195. ematics and his PhD in Industrial Engineering from
Law, A. M., and W. D. Kelton. 1991. Sim'ulation Purdue lTniversity. He has also taught courses at
modeling and analysis. 2d ed. Ne\\' \'"ork: i\lc- Baylor University, The lTniversity of Oklahoma, and
Graw-Hill. Purdue University. His consulting, short course, and
Law, A. M., M. G. McComas, and S. G. Vincent. research contract work includes contracts with AT&T,
1994. The crucial role of input 1110deling in suc- NASA/Langley Research Center, Delco Electronics,
cessful simulation studies. Industrial Engineering Departlnent of Defense, Air Logistic COlnmand, leASE,
26:55-59. !(omag, Federal Aviation Administration, Tinker Air
Lawless, J. F. 1982. Statistical models 6' m,ethods for Force Base, Magnetic Peripherals, and Argonne Na-
lifetime data. New York: John Wiley &: Sons. tional Laboratory. His research and teaching interests
Lee, S., J. R. Wilson, and M. M. Crawford. 1991. are in reliability and simulation. He is a member of
Modeling and simulation of a nonho1110geneous ASA, liE, and INFORMS.
Poisson process having cyclic behavior. Com.mu-
nications in Statistics - Simulation and Compu-
tation 20:777-809.
Leelnis, L. M. 1991. Nonparametric esti111ation of the
intensity function for a nonho1110geneous Poisson
process. Management Science 37:886-900.
Lieblein, J., and M. Zelen. 1956. Statistical investi-
gation of the fatigue life of deep-groove ball bear-
ings. Journal of Research of the National Bureau
of Standards 57:273-316.
Martz, H. F., and R. A. Waller. 1982. Bayesian reli-
ability analysis. New yTork: John Wiley & Sons.
Nelson, B. L., P. Ware, M. C. Cario, C. A. Harris,
S. A. Jamison, J. O. Miller, J. Steinbugl, J. \Tang .