Proceedings of the 2015 Winter Simulation Conference

L. Yilmaz, W. K. V. Chan, I. Moon, T. M. K. Roeder, C. Macal, and M. D. Rossetti, eds.

STATISTICAL ANALYSIS OF SIMULATION OUTPUT DATA:

THE PRACTICAL STATE OF THE ART

Averill M. Law

Averill M. Law & Associates, Inc.

4729 East Sunrise Drive, #462
Tucson, AZ 85718, USA

ABSTRACT

One of the most important but neglected aspects of a simulation study is the proper design and analysis of
simulation experiments. In this tutorial we give a state-of-the-art presentation of what the practitioner re-
ally needs to know to be successful. We will discuss how to choose the simulation run length, the
warmup-period duration (if any), and the required number of model replications (each using different
random numbers). The talk concludes with a discussion of three critical pitfalls in simulation output-data
analysis.

1 INTRODUCTION
In many “simulation studies” a great amount of time and money is spent on model development and
“programming,” but little effort is made to analyze the simulation output data appropriately. As a matter
of fact, a common mode of operation is to make a single simulation run of somewhat arbitrary length and
then to treat the resulting simulation estimates as the “true” model characteristics. Since random samples
from probability distributions are typically used to drive a simulation model through time, these estimates
are just particular realizations of random variables that may have large variances. As a result, these esti-
mates could, in a particular simulation run, differ greatly from the corresponding true characteristics for
the model. The net effect is, of course, that there could be a significant probability of making erroneous
inferences about the system under study.
We now describe more precisely the random nature of simulation output. Let Y1 , Y2 , be an output
stochastic process [see, for example, section 4.3 in Law (2015)] from a single simulation run. For exam-
ple, Yi might be the delay in queue for the ith job to arrive at a single-server queueing system. Alterna-
tively, Yi might be the total cost of operating an inventory system in the ith month. The Yi ' s are random
variables that will not, in general, be independent or identically distributed (IID). Thus, many of the for-
mulas from classical statistics (see Section 2) will not be directly applicable to the analysis of simulation
output data.
Example 1. For the queueing system mentioned above, the delays in queue will not be independent,
since a large delay for one customer waiting in queue will tend to be followed by a large delay for the
next customer waiting in queue. Suppose that the simulation is started at time zero with no customers
in the system, as is usually the case. Then the delays in queue at the beginning of the simulation will
tend to be smaller than later delays and, thus, the delays are not identically distributed.
Let y11 , y12 ,..., y1m be a realization of the random variables Y1 , Y2 ,, Ym resulting from running the
simulation with a particular set of random numbers u11 , u12 ,... . If we run the simulation with a different

978-1-4673-9743-8/15/$31.00 ©2015 IEEE 1810

 Law

set of random numbers u21 , u22 ,... , then we will obtain a different realization y21 , y22 ,..., y2 m of the ran-
dom variables Y1 , Y2 ,, Ym . (The two realizations are not the same since the different random numbers
used in the two runs produce different samples from the input probability distributions.) In general, sup-
pose that we make n independent replications (runs) of the simulation (i.e., different random numbers are
used for each replication, each replication uses the same initial conditions, and the statistical counters for
the simulation are reset at the beginning of each replication) each of length m, resulting in the observa-
tions:

y11 ,..., y1i ,..., y1m

y 21 ,..., y 2 i ,..., y 2 m

  
y n1 ,..., y ni ,..., y nm

The observations from a particular replication (row) are clearly not IID. However, note that y1i , y 2 i ,..., y ni
(from the ith column) are IID observations of the random variable Yi , for i  1,2,..., m. More generally,
each entire replication is independent of any other replication, and each replication’s observations have
the same (joint) distribution. This independence across runs is the key to relatively simple output-data
analysis that is discussed in later sections of this paper. Then, roughly speaking, the goal of output-data
analysis is to use the observations y ji (i = 1, 2, …, m; j = 1, 2, …, n) to draw inferences about character-
istics of the random variables Y1 , Y2 ,, Ym .

Example 2. Consider a bank with five tellers and one queue, which opens its doors at 9 A.M., closes
its doors at 5 P.M., but stays open until all customers in the bank at 5 P.M. have been served. As-
sume that customers arrive with IID exponential inter-arrival times with mean 1 minute, that service
times are IID exponential random variables with mean 4 minutes, and that customers are served in a
first-in, first-out (FIFO) manner. Table 1 shows two typical output statistics from 5 independent rep-
lications of the bank, assuming that no customers are present initially. Note that results from differ-
ent replications can be quite different. Thus, one run clearly does not produce the “answers.”

Table 1: Results for 5 independent replications of the Bank Model.

Replication Average delay Average number
in queue in queue
1 1.53 1.52
2 1.66 1.62
3 1.24 1.23
4 2.34 2.34
5 2.86 2.83

Our goal in this paper is to discuss methods for statistical analysis of simulation output data and to
present the material with a practical focus. Section 2 of this paper reviews formulas from classical statis-
tics based on IID data, which we will find useful later in this paper. In Section 3, we discuss the two main
types of simulations with regard to output-data analysis, namely, terminating and non-terminating. Statis-
tical methods for analyzing each type are given in Sections 4 and 5, respectively. Finally, we give a
summary of this tutorial and three fundamental pitfalls in output-data analysis in Section 6.
Portions of this paper are based on chapters 4 and 9 of Law (2015). Other references on output-data
analysis are Alexopoulos (2007), Banks et al. (2010), Currie and Cheng (2013), and Nakayama (2008).

1811
 Law

2 REVIEW OF CLASSICAL STATISTICS

Suppose that X 1 , X 2 ,, X n are IID random variables with population mean and variance  and  2 , re-
spectively. Then unbiased point estimators for  and  2 are given by

X
n

X (n)  i 1
i
(1)
n

and

[ X  X (n)]2
n

S (n)  i 1
i

n 1
2 (2)

proximate 100(1   ) percent (0    1) confidence interval for  is given by

(An estimator is unbiased if its mean is equal to the target population characteristic.) Furthermore, an ap-

X (n)  tn1,1 /2 S 2 (n) / n (3)

where tn 1,1 / 2 is the upper 1   / 2 critical point for a t distribution with n  1 degrees of freedom. If

coverage probability arbitrarily close to 1   . Alternatively, if the X i ' s are normally distributed, then
the sample size n is “sufficiently large,” then the confidence interval given by Expression (3) will have a

the coverage probability will be exactly 1   . In practice, if the distribution of the X i ' s is reasonably
symmetric, then the coverage probability will be close to 1   [see Law (2015, pp. 233-237)]. If we in-
crease the sample size from n to 4n, then the half-length of the confidence interval, tn1,1 /2 S 2 (n) / n ,
will decrease by a factor of approximately 2, since there is an n in the denominator under the square-root
sign.
As stated above, the Yi ' s from one simulation run are not IID and, thus, Expressions (1), (2), and (3)
are not directly applicable to their analysis. However, if we take comparable output statistics from differ-
ent independent replications of a simulation model, then these observations are IID and the three expres-
sions are applicable.

Example 3. For the bank simulation of Example 2, the five average delays in queue from column 2 of
Table 1 are IID and, thus, Expressions (1), (2), and (3) could legitimately be used for their analysis.

3 TYPES OF SIMULATIONS WITH REGARD TO OUTPUT ANALYSIS

The options available for designing and analyzing simulation experiments depend on whether the simula-
tion of interest is terminating or non-terminating, which depends on whether there is an obvious way for
determining the simulation run length.
A terminating simulation is one for which there is a “natural” event E that specifies the length of each
run (replication). Since different runs use independent random numbers and the same initialization rule,
this implies that comparable random variables are IID. The event E often occurs at a time point that has
one of the following properties:

1812
 Law



The system is “cleaned out”


Beyond which no useful information is obtained
Specified by management.

The event E is specified before any runs are made, and the time of occurrence of E for a particular run
may be a random variable. Since the initial conditions for a terminating simulation generally affect the
desired measures of performance, these conditions should be representative of those for the actual system.

Example 4. A retail/commercial establishment (e.g., a bank) closes each evening. If the establish-
ment is open from 9 to 5, the objective of a simulation might be to estimate some measure of the qual-
ity of customer service over the period beginning at 9 A.M. and ending when the last customer who
entered before the doors closed at 5 P.M. has been served. In this case, E = {8 hours of simulated
time have elapsed and the system is empty}, and the initial conditions for the simulation should be
representative of those for the bank at 9 A.M.

Example 5. Consider a military ground confrontation between a blue force and a red force. Relative
to some initial force strengths, the goal of a simulation might be to determine the (final) force
strengths when the battle ends. In this case, E = {either the blue force or the red force has “won” the
battle}. An example of a condition that would end the battle is one side losing 30 percent of its force,
since this side would no longer be considered viable. The choice of initial conditions for the simula-
tion, e.g., the number of troops and tanks for each force, is generally not a problem here, since they
are specified by the military scenario under consideration.

A non-terminating simulation is one for which there is no natural event E to specify the length of a
run. This often occurs when we are designing a new system or modifying an existing system, and we are
interested in the behavior of the system in the long run when it is operating ”normally.” Unfortunately,

Consider the output stochastic process Y1 , Y2 ,... for a simulation model. Let Fi ( y  I )  P(Yi  y  I ) for
“in the long run” doesn’t naturally translate into a terminating event E.

i  1,2,... , where y is a real number and I represents the initial conditions used to start the simulation at
time 0. [The conditional probability P(Yi  y  I ) is the probability that the event {Yi  y} occurs given the

each machine is busy or idle, at time 0. We call Fi ( y  I ) the transient distribution of the output process
initial conditions I.] For a manufacturing system, I might specify the number of jobs present, and whether

at (discrete) time i for initial conditions I. Note that Fi ( y  I ) will, in general, be different for each value
of i and each set of initial conditions I. For fixed y and I, the probabilities F1 ( y  I ) , F2 ( y  I ) , … are just
a sequence of numbers. If Fi ( y  I )  F ( y ) as i   for all y and any initial conditions I, then F ( y) is
called the steady-state distribution of the output process Y1 , Y2 ,... . Note that the steady-state distribution
F ( y) does not depend on the initial conditions I.
A measure of performance for a non-terminating simulation is said to be a steady-state parameter if it
is a characteristic of the steady-state distribution of some output stochastic process Y1 , Y2 ,.... If the ran-

state mean   E (Y ) .
dom variable Y has the steady-state distribution, then we are typically interested in estimating the steady-

Example 6. Consider a company that is going to modify its manufacturing system and would like to
determine the long-run (steady-state) mean hourly throughput of their system after it has been running
long enough for workers to know their jobs and for mechanical difficulties to have been worked out.

1813
 Law

The system operates continuously 24 hours a day for 7 days a week. Let N i be the number of parts
manufactured in the ith hour. If the stochastic process N1 , N 2 ,... has a steady-state distribution with

  E ( N ).
corresponding random variable N, then they are interested in estimating the steady-state mean

Example 7. Suppose that a military organization is going to employ a new inventory system during (a
long) peacetime and would like to determine the long-run mean monthly cost of operating their sys-
tem.
Let Ci be the cost of operating the inventory system in the ith month. If the output process C1, C2, ...

mating the mean   E (C ).

has a steady-state distribution with corresponding random variable C, then they are interested in esti-

4 STATISTICAL ANALYSIS FOR TERMINATING SIMULATIONS

event E. Let X j be an output random variable defined over the jth replication, for j  1, 2,..., n ; it is as-
Suppose that we make n independent replications of a terminating simulation each terminated by the

sumed that the X j ' s are comparable for different replications. Then the X j ' s are IID random variables.

 D / N over a day from the jth replication,

N
For the bank of Example 4, X j might be the average delay
i 1
i

where N (a random variable) is the number of customers served in a day and Di is the delay in queue of
the ith arriving customer. For the combat model of Example 5, X j might be the number of red tanks de-
stroyed on the jth replication.

  E ( X ), where X is a random variable defined on a replication as described above. Make n independ-

Suppose that we would like to obtain a point estimate and confidence interval for the mean

ent replications of the simulation and let X 1 , X 2 ,, X n be the resulting IID random variables. Then, by
substituting the X j ' s into Expressions (1), (2), and (3), we get that X ( n) is an unbiased point estimator
for  , and an approximate 100(1   ) percent confidence interval for  is given by

X (n)  tn1,1 /2 S 2 (n) / n

Example 8. A small factory consists of a machine and an inspector, as shown in Figure 1. Unfin-
ished parts arrive to the factory with exponential inter-arrival times having a mean of 1 minute. Pro-
cessing times at the machine are uniformly distributed on the interval [0.65, 0.70] minute, and subse-
quent inspection times at the inspector are uniformly distributed on the interval [0.75, 0.80] minute.
(The assumption of uniformity is for ease of exposition, and is not likely to be valid in a real-world
application.) Ninety percent of inspected parts are “good” and leave the system immediately; 10 per-
cent of the parts are “bad” and are sent back to the machine for rework. (Both queues are assumed to
be of infinite capacity.) The machine is subject to randomly occurring breakdowns. In particular, a
new (or freshly repaired) machine will break down after an exponential amount of calendar time with
a mean of 6 hours. Repair times are uniform on the interval [8, 12] minutes. If a part is being pro-
cessed when the machine breaks down, then the machine continues where it left off upon the comple-
tion of repair. Assume that the factory is initially empty and idle.

1814
 Law

0.9 Good
Machine Inspector

0.1 Bad

Figure 1: Small factory.

The factory gets an order to produce 2000 parts and, thus, a simulation of this system can be con-
sidered to be terminating with E = {2000 parts have been completed}. Let T be the time required to

confidence interval for the mean   E (T ).

complete the required 2000 parts. Then the company would like a point estimate and a 95 percent

We made 10 independent replications of the simulation and obtained the following observed val-
ues for T (in hours):

T1  32.62 , T2  32.57 , T3  33.51 , T4  33.29 ,

T5  32.10 , T6  34.24 , T7  32.70 , T8  33.49 ,
T9  33.36 , T10  34.61

Substituting the T j ' s into Expressions (1), (2), and (3), gives the following results:

T (10)  33.25, S 2 (10)  0.606

and an (approximate) 95 percent confidence interval for   E (T ) is given by

33.25  0.56 or [32.69, 33.81]

Thus, we are approximately 95 percent confident that  is between 32.69 and 33.81 hours. (If

100 confidence intervals to contain the true  .) Note also that the interval is quite precise, with the
100 people performed this experiment independently, then we would expect that about 95 out of the

half-length of the confidence interval being less than 2 percent of the point estimate.
We discussed above how to estimate the mean of the random variable T. However, in some cases
we might also be interested in getting a measure of how large a value T can take on. For example, we
may want to estimate the 95th percentile, x0.95 , of the distribution of the random variable T, i.e., T is
larger than x0.95 only 5 percent of the time.

5 STATISTICAL ANALYSIS FOR NONTERMINATING SIMULATIONS

that we want to estimate the steady-state mean   E (Y ) , which is also defined by

Let Y1 , Y2 ,... be an output stochastic process from a single run of a non-terminating simulation. Suppose

  lim E (Yi )
i 

However, E (Yi )   for “small” i because we generally don’t know how to choose the initial conditions I
where E (Yi ) is the transient mean at time i. Thus, the transient means converge to the steady-state mean.

1815
 Law

tor of  for all finite values of m. The problem that we have just described is called the problem of the
to be representative of “steady-state behavior.” This causes the sample mean Y ( m) to be a biased estima-

initial transient in the simulation literature.

The technique most often suggested for dealing with this problem is called warming up the model.

ing observations to estimate  . In particular, given the observations Y1 , Y2 ,, Ym , we would use
The idea is to delete some number of observations from the beginning of a run and to use only the remain-

Y
m

Y (m, l )  i l 1
i

ml

(1  l  m  1) rather than Y (m) as an estimator of  . In general, one would expect Y ( m, l ) to be less bi-
ased than Y ( m) , since the observations near the “beginning” of the simulation may not be very repre-
sentative of steady-state behavior due to the choice of initial conditions.

would like to pick l (and m) such that E[Y (m, l )]   . If l and m are chosen too small, then E[Y ( m, l )]
The question naturally arises as to how to choose the warmup period (or deletion amount) l. We

may be significantly different than  . On the other hand, if l is chosen larger than necessary, then
Y ( m, l ) will probably have an unnecessarily large variance.

(1983) [see also Law (2015, pp. 513-520)]. Its specific goal is to determine l such that E (Yi )   for
The simplest and most general technique for determining l is a graphical technique due to Welch

i  l , where l is the warmup period. This is equivalent to determining when the transient-mean curve
E (Yi ) “flattens out” at level  . In general, it is difficult to determine l from a single replication due to the
inherent variability of the process Y1 , Y2 ,... . As a result, Welch’s procedure is based on making multiple
replications of the simulation in a pilot study.

5.1 The Replication/Deletion Approach

In this section, we discuss how to construct a point estimate and confidence interval for  . Suppose that
the warmup period has been determined by Welch’s procedure or by using “engineering judgment.”
Make n independent replications of the output process Y1 , Y2 ,... each of length m, where m should be

nating simulations.) Let Y ji be the ith observation from the jth replication, for j  1,2,..., n and
much larger than l. (There is no definitive way of picking the run length m here, as there was for termi-

i  1,2,..., m . Let

Xj  Y /(m - l ) for j  1, 2,..., n

m

i  l 1
ji

Note that i  l  1 is where we think that “steady state” begins. Then the X j ' s are IID random variables.
Furthermore, E ( X j )   since Y j ,l 1 , Y j ,l  2 ,..., Y j , m each have approximate mean  . Then, by substituting
the X j ' s into Expressions (1), (2), and (3), we get that X ( n) is an (approximately) unbiased point esti-
mator for  , and an approximate 100(1   ) percent confidence interval for  is given by

X (n)  tn1,1 /2 S 2 (n) / n

1816
 Law

We call the above method for constructing a point estimate and confidence interval for  the replica-
tion/deletion method. One criticism that was levied against this method historically is that l observations
must be discarded from each of the n replications. However, given the availability of PCs with fast multi-
core processors, this is no longer an issue for many, if not most, steady-state analyses.

Example 9. Consider a manufacturing system with a receiving/shipping station and five workstations
(see Figure 2), as described in Law (2015, chapter 14). Assume that there are 4, 2, 5, 3, and 2 ma-
chines in stations 1through 5, respectively. The machines in a particular station are identical, but ma-
chines in different stations are dissimilar. Jobs arrive to the system with exponential inter-arrival
times having a mean of 1/15th an hour. Thus, 15 jobs arrive in a typical hour. There are three types
of jobs, and jobs are of types 1, 2, and 3, with respective probabilities 0.3, 0.5, and 0.2. Job types 1,
2, and 3 require 4, 3, and 5 operations to be done, respectively, and each operation must be done at a
specified workstation in a prescribed order. Each job begins at the receiving/shipping station, travels
to the work stations on its routing, and then leaves the system at the receiving/shipping station. For
example, the routing for a type 1 job is 3, 1, 2, 5.
A job must be moved from one station to another by a forklift truck, which moves at a speed of 5
feet per second, and two forklifts are available. When a forklift becomes available, it processes re-
quests by jobs using a shortest-distance-first dispatching rule. If more than one forklift is idle when a
job requests transport, then the closest forklift is used. When a forklift finishes moving a job to a
work station, it remains at that station if there are no pending job requests. Each station has a single
FIFO

Workstation 2 Workstation 3 Workstation 4

Forklift

In Out

Workstation 1 Receiving/Shipping Workstation 5

Figure 2: Factory with five workstations.

queue. The time to perform an operation at a particular machine is a gamma random variable with a
shape parameter of 2, whose mean depends on the job type and the station to which the machine be-
longs. For example, the mean service time for a type 1 job at station 3 (the first station on its routing)
is 0.25 hour. When a machine finishes processing a job, the job blocks that machine (i.e., the ma-
chine cannot process another job) until the job is removed by a forklift.
The factory is open 8 hours a day, and thus the arrival rate is 120 jobs per day. The system con-
figuration described here is called system design 3 in Law (2015).
Let N1 , N 2 ,... be the output stochastic process corresponding to daily throughputs. Then we are

daily throughput   E ( N ). Since the simulation starts out with no jobs present at time zero, the
interested in obtaining a point estimate and 90 percent confidence interval for the steady-state mean

throughput will tend to be “small” during the early part of a run and a warmup period is needed. Us-

cess is l  15 days [see Law (2015, chapter 14)].

ing Welch’s graphical procedure, we determined that a reasonable warmup period for this output pro-

1817
 Law

We made n  10 (production) replications of length m  115 days, and used a warmup period of
l  15 days. Let

N
115

Xj  i 16
ji

100

where N ji is the throughput in the ith day of the jth replication.

Substituting the X j ' s into Expressions (1), (2), and (3), we get the following point estimate and
approximate 90 percent confidence interval for   E ( N ) :

ˆ  X (10)  120.29

and

120.29  0.63 or [119.66,120.92]

Thus, subject to the correct interpretation, we are approximately 90 percent confident that the steady-
state mean daily throughput is between 119.66 and 120.92 jobs per day. Note that this confidence in-

and forklifts. (In a real application,  would not, of course, be known.)

terval contains 120, which should be the mean daily throughput if the system has enough machines

Note also that the confidence interval is quite precise, with the half-length being less than 1 per-
cent of the point estimate. Also, since X j is the average of 100 N ji ' s , it should be approximately
normally distributed by a central-limit-theorem type effect. This suggests that the coverage of the
confidence interval should be close to the desired coverage probability of 0.9. Finally, if, for exam-
ple, we wanted to decrease the half-length by a factor of 3, then a total of approximately 90 replica-
tions would be required.

6 SUMMARY AND PITFALLS IN OUTPUT-DATA ANALYSIS

We have seen that both terminating and non-terminating analyses can be performed easily by making in-
dependent replications of the simulation model and using Expressions (1), (2), and (3), which come from
a first undergraduate course in statistics. In the case of steady-state parameters, we also have to determine
a warmup period, but this can be reliably addressed using Welch’s graphical approach. The method of
replication can also be easily applied to comparing alternative system configurations [see, for example,
Law (2015, chapters 10 and 11)] and to estimating multiple measures of performance. Moreover, multi-
ple replications can be made simultaneously on computers having multiple-core processors, on computers
connected by a local-area network, or by using cloud computing!

The following are three major pitfalls in output-data analysis:

 Analyzing simulation output data from one run using formulas [e.g., Expression (2)] that assume
independence, which might result in a gross underestimation of variances and standard devia-
tions. This problem is exacerbated by the use of these formulas by some simulation-software


packages.
Failure to have a warmup period for steady-state analyses

1818
 Law

 Failure to determine the statistical precision of simulation output statistics by the use of a confi-
dence interval, which can be accomplished easily using the replication approach.

REFERENCES

Alexopoulos, C. 2007. “Statistical Analysis of Simulation Output: State of the Art.” In Proceedings of the
2007 Winter Simulation Conference, edited by S. G. Henderson, B. Biller, M.-H. Hsieh, J. Shortle, J.
D. Tew, and R. R. Barton, 150-161. Piscataway, New Jersey: Institute of Electrical and Electronic
Engineers, Inc.
Banks, J., J. S. Carson, B. L. Nelson, and D. M. Nicol. 2010. Discrete-Event System Simulation, 5th ed.
Upper Saddle River, New Jersey: Prentice-Hall.
Currie, C. S. M., and R. C. H. Cheng. “A Practical Introduction to Analysis of Simulation Output Data.”
In Proceedings of the 2013 Winter Simulation Conference, edited by R. Pasupathy, S.-H. Kim, A.
Tolk, R. Hill, and M. E. Kuhl, 328-341. Piscataway, New Jersey: Institute of Electrical and Electronic
Engineers, Inc.
Law, A. M. 2015. Simulation Modeling and Analysis, 5th ed. New York: McGraw-Hill.
Nakayama, M. K. 2008. “Statistical Analysis of Simulation Output.” In Proceedings of the 2008 Winter
Simulation Conference, edited by S. J. Mason, R. R. Hill, L. Mönch, O. Rose, T. Jefferson, J. W.
Fowler, 62-72. Piscataway, New Jersey: Institute of Electrical and Electronic Engineers, Inc.
Welch, P. D. 1983. “The Statistical Analysis of Simulation Results.” In The Computer Performance Mod-
eling Handbook, edited by S. S. Lavenberg. New York: Academic Press.

AUTHOR BIOGRAPHY

AVERILL M. LAW is President of Averill M. Law & Associates, a company specializing in simulation
seminars, simulation consulting, and software. He has presented more than 555 simulation and statistics
short courses in 20 countries, including onsite seminars for AT&T, Boeing, Caterpillar, Coca-Cola, GE,
GM, IBM, Intel, Lockheed Martin, Los Alamos National Lab, NASA, NATO (Netherlands), NSA, Sasol
Technology (South Africa), 3M, UPS, U.S. Air Force, U.S. Army, U.S. Navy, and Verizon. Dr. Law has
been a simulation consultant to more than 50 organizations including Booz Allen & Hamilton, Cono-
coPhillips, Defense Modeling and Simulation Office, Kimberly-Clark, M&M/Mars, Oak Ridge National
Lab, U.S. Air Force, U.S. Army, U.S. Marine Corps, and U.S. Navy. He has written or coauthored nu-
merous papers and books on simulation, operations research, statistics, manufacturing, and communica-
tions networks, including the book Simulation Modeling and Analysis that has more than 157,000 copies
in print and 14,800 citations. He developed the ExpertFit® distribution-fitting software and also several
videotapes on simulation modeling. He was awarded the INFORMS Simulation Society Lifetime Profes-
sional Achievement Award in 2009. Dr. Law wrote a regular column on simulation for Industrial Engi-
neering magazine. He has been a tenured faculty member at the University of Wisconsin-Madison and the
University of Arizona. He has a Ph.D. in industrial engineering and operations research from the Univer-
sity of California at Berkeley. His e-mail address is <averill@[Link]> and his website is
<[Link]>.

1819