9
Input Modeling
Input data provide the driving force for a simulation model. In the simulation of
@ queueing system, typical input data are the distributions of time between ar-
tivals and service times, For an inventory system simulation, input data include
the distributions of demand and lead time. For the simulation of a reliability
system, the distribution of time-to-failure of a component is an example of
input data.
In the examples and exercises in Chapters 2 and 3, the appropriate distri-
butions were specified for you. In real-world simulation applications, however,
determining appropriate distributions for input data is a major task from the
standpoint of time and resource requirements. Regardless of the sophistication
of the analyst, faulty models of the inputs will lead to outputs whose interpre-
tation may give rise to misleading recommendations.
‘There are four steps in the development of a useful model of input data:
1, Collect data from the real system of interest. This often requires asubstan-
tial time and resource commitment. Unfortunately, in some situations it is
not possible to collect data (for example, when time is extremely limited,
when the input process does not yet exist, or when laws or rules prohibit
the collection of data). When data are not available, expert opinion and
knowledge of the process must be used to make educated guesses.
. Identify a probability distribution to represent the input process. When
data are available, this step typically begins by developing a frequency
distribution, or histogram; of the data. Based on the frequency distri-
bution and structural knowledge of the process, a family of distributions
is chosen. Fortunately, as described in Chapter 5, several well-known
distributions often provide good approximations in practice.
ny
323324° Chap.9 Input Modeling
3. Choose parameters that determine a Specific instance of the distribution
family. When data are available, these parameters may be estimated from
the data.
. Evaluate the chosen distribution afd the associated parameters for good-
ness-of-fit. Goodness-of-fit may be evaluated informally via graphical
methods, or formally via statistical tests, The chi-square and the Kolmo-
gorov-Smirnov tests are standard goodness-of-fit tests. If not satisfied
that the chosen distribution is a good approximation of the data, then the
analyst returns to the second step, chooses a different family of distribu-
tions, and repeats the procedure. If several iterations of this procedure
fail to yield a fit between an assumed distributional form and the collected
data, the empirical form of the distribution may be used as described in
Section 8.1.5 of the previous chapter.
-
Each of these steps is discussed in this chapter. Although software is now widely
available to accomplish st "ps 2, 3 and 4 — including standalone programs such
as ExpertFit® and Stat::Fit™, and integrated programs such as Arena’s Input
Processor and @Risk’s BestFit — it is still important to understand what the
software does so that it can be used appropriately. Unfortunately, software is
not as readily available for input modeling when there is a relationship between
two or more variables of interest, or when no data are available. These two
topics are discussed toward the end of the chapter.
9.1 Data Collection
Problems are found at the end of each chapter, as exercises for the reader,
in mathematics, physics, chemistry, and other technical subject texts. Years
and years of working these problems may give the reader the impression that
data are readily available. Nothing could be further from the truth. Data
collection is one of the biggest tasks in solving a real problem. It is one of the
most important and difficult problems in simulation. And even when data are
available, they have rarely been recorded in a form that is directly useful for
simulation input modeling.
“GIGO,” or “garbage-in, garbage-out,” is a basic concept in computer
science and it applies equally in the area of discrete system simulation. Many
are fooled by a pile of computer output or a sophisticated animation, as if these
were the absolute truth. Even if the model structure is valid, if the input data
are inaccurately collected, inappropriately analyzed, or not representative of
the environment, the simulation output data will be misleading and possibly
damaging or costly when used for policy or decision making,
EXAMPLE 9.1 (The Laundromat)
As budding simulation students, the first two authors had assignments to sim-
ulate the operation of an ongoing system. One of these systems, which seemed
to be a rather simple operation, was a self-service laundromat with 10 washing
machines and six dryers.Sec. 9.1 Data Collection 325
However, the data-collection aspect of the problem rapidly became rather
enormous. The interarrival-time distribution was not homogeneous; that is, the
distribution changed by time of day and by day of week. The laundromat was
open 7 days a week for 16 hours per day, or 112 hours per week. It would have
been impossible to cover the operation of the laundromat with the limited
resources available (two students who were also taking four other courses) and
Service-time distributions also presented a difficult problem from many
perspectives, The proportion of customers demanding the various service com-
desiring one washer followed by one dryer. However, a customer might choose
wo washing machines followed by one dryer, one dryer only, and so on. Since
the customers used numbered machines, it was possible to follow them us-
ing that reference, rather than remembering them by personal characteristics.
Because of the dependence between washer demand and dryer demand for
an individual customer, it would have been inappropriate to treat the service
times for washers and dryers separately as independent variables,
Some customers waited patiently for their clothes to complete the washing
or drying cycle, and then they removed their clothes Promptly. Others left the
Premises and returned after their clothes had finished their cycle on the machine
being used. Ina very busy period, the manager would remove a customer’s
clothes after the cycle and set them aside in a basket. It was decided that
service termination would be measured as the point in time when the machine
was emptied of its contents.
Also, machines would break down from time to time. The length of
the breakdown varied from a few moments, when the manager repaired the
machine, to several days (a breakdown on Friday night, requiring a part not
in the laundromat storeroom, would not be fixed until the following Monday).
The short-term repair times were tecorded by the student team. The long-term
repair completion times were estimated by the manager. Breakdowns then
became part of the simulation, <
Many lessons can be learned from an actual experience in data collection.
The first five exercises at the end of this chapter suggest some situationsin which
the student can gain such experience.
The following suggestions may enhance and facilitate data collection, al-
though they are not all-inclusive.
1. A useful expenditure of time is in planning. This could begin by a practice
Or preobserving session. Try to collect data while preobserving. Devise
forms for this purpose. It is very likely that these forms will have to be326
»
Chap.9 Input Modeling
modified several times before the actual data collection begins. Watch
for unusual circumstances and consider how they will be handled. When
Possible, videotape the system and extract the data later by viewing the
tape. Planning is important even if data will be collected automatically
(e.g., via computer data collection) to insure that the appropriate data are
available. When data have already been collected by someone else, be
sure to allow plenty of time for converting the data into a usable format,
‘Try to analyze the data as they are being collected. Determine if the data
being collected are adequate to provide the distributions needed as input
to the simulation. Determine if any data being collected are useless to
the simulation. There is no need to collect superfluous data.
. Try to combine homogeneous data sets. Check data for homogeneity in
Successive time periods and during the same time period on successive
days. For example, check for homogeneity of data from 2:00 pM. to 3:00
PM. and 3:00 PM. to 4:00 pM. and check to see if the data are homogeneous
for 2:00 pa. to 3:00 eM, on Thursday and Friday. When checking for
homogeneity, an initial test is to see if the means of the distributions (the
average interarrival times, for example) are the same. The two-sample t
test can be used for this purpose. A more thorough analysis would require
a determination of the equivalence of the distributions using, perhaps, a
quantile-quantile plot (described later).
Be aware ofthe Possibility of data censoring, in which aquantity ofinterest
is not observed in its entirety. This problem most often occurs when the
analyst is interested in the time required to complete some Process (for
example, produce a part, treata patient, or have acomponent fail), but the
Process begins prior to, or finishes after the completion of, the observation
Period. Censoring can result in especially long process times being left
Out of the data sample.
To determine whether there is a relationship between two variables, build
a scatter diagram. Sometimes an eyeball scan of the scatter diagram
will indicate if there is a relationship between two variables of interest.
Section 9.6 describes models for statistically dependent input data.
Consider the possibility that a sequence of observations which appear to
be independent may possess autocorrelation. Autocorrelation may exist
in successive time periods or for successive customers. For example, the
service time for the ith customer may be related to the service time for the
(+ n)th customer. A brief introduction to autocorrelation was provided
in Section 7.4.3, and some input models that account for autocorrelation
are presented in Section 9.6.
Keep in mind the difference between input data and output or perfor-
mance data, and be sure to collect input data. Input data typically rep-
resent the uncertain quantities that are largely beyond the control of the
system and will not be altered by changes made to improve the system.
Output data, on the other hand, represent the performance of the sys-
tem when subjected to the inputs, performance that we may be trying toSec. 9.2 Identifying the Distribution with Data 327
improve. Ina queueing simulation, the customer arrival times are usu-
ally inputs, while the customer delay is an output. Performance data are
useful for model validation, however (see Chapter 10).
Again, these are just a few suggestions, As a rule, data collection and analysis
must be approached with great care,
9.2 Identifying the Distribution with Data
In this section we discuss methods for selecting families of input distributions
when data are available. The specific distribution within a family is specified
by estimating its parameters, as described in Section 9.3. Section 9.5 takes up
the case when no data are available.
9.2.1 Histograms
A frequency distribution or histogram is useful in identifying the shape of a
distribution. A histogram is constructed as follows:
1. Divide the range of the data into intervals (intervals are usually of equal
width; however, unequal widths may be used if the heights of the frequen-
cies are adjusted).
2, Label the horizontal axis to contorm to the intervals selected.
3. Determine the frequency of occurrences within each interval.
4. Label the vertical axis so that the total occurrences can be plotted for
each interval.
5. Plot the frequencies on the vertical axis.
‘The number of class intervals depends on the number of observations and the
amount of scatter or dispersion in the data. Hines and Montgomery [1990]
state that choosing the number of class intervals approximately equal to the
Square root of the sample size often works well in practice. If the intervals are
too wide, the histogram will be coarse, or blocky, and its shape and other details
will not show well. If the intervals are too narrow, the histogram will be ragged
and will not smooth the data. Examples of a ragged, coarse, and appropriate
histogram using the same data are shown in Figure 9.1. Modern data-analysis
software often allows the interval sizes to be changed easily and interactively
until a good choice is found.
The histogram for continuous data corresponds to the probability density
function of a theoretical distribution. If continuous, a line drawn through the
center point of each class interval frequency should result in a shape like that
ofa pdf.
Histograms for discrete data, where there are a large number of data
points, should have a cell for each value in the range of the data. However,
if there are few data points, it may be necessary to combine adjacent cells to
eliminate the ragged appearance of the histogram. Ifthe histogram is associated
with discrete data, it should look like a probability mass function.328 Chap. 9 Input Modeling
Frequency
024 6 8 012 14 16 18 2022 ©
@
Pieb
a
gab
a
8
ab
=
oF Hs Te :
o
4
nb
10
Bost
a
Boob
E
AG
2
20
02 35 68 9-11 12-14 15-17 18-20 21-24 x
©)
Figure 9.1. Ragged, coarse, and appropriate histograms:
(a) original data — too ragged; (b) combining adjacent cells
— too coarse; (c) combining adjacent cells — appropriate.
EXAmpLe 9.2 (Discrete Data)
‘The number of vehicles arriving at the northwest corner of an intersection in a5-
minute period between 7:00 a.m. and 7:05 A.M. was monitored for five workdays
over a 20-week period. Table 9.1 shows the resulting data. The first entry in the
table indicates that there were 12 5-minute periods during which zero vehicles
arrived, 10 periods during which one vehicle arrived, and so on.
Since the number of automobiles is a discrete variable, and since there
are ample data, the histogram can have a cell for each possible value in the
range of the data. The resulting histogram is shown in Figure 9.2. <Sec. 9.2. Identifying the Distribution with Data
Table 9.1, Number of Arrivals in a 5-
Minute Period
Arrivals Arrivals,
per Period Frequency _per Period _ Frequency
0 2 6 7
1 10 7 3
2 19 8 $
3 a 9 3
4 10 10 3
5 8 u 1
20
18
16
pi
3
212
E
10
8
6
4
2
o4123 45678 9 0
Figure 9.2. Histogram of number of arrivals per period.
‘Number of arrivals per period
Example 9,3 (Continuous Data)
Life tests were performed on a random samp!
329
le of electronic chips at 1.5 times
the nominal voltage, and their lifetime (or time to failure) in days was recorded:
79.919
3.027
6.769
18.387,
144.695
0.941
0.624
0.590
7.004
3.217
3.081
6.505
59.899
0.141
2.663
0.878
5.380
1.928
31.764
14.382
0.062
0.021
1,192
43.565
17.967
3.371
3,148
0.300
1.005
1.008
1.961
0.013
34.760
24.420
0.091
2.157
7.078
0.002
1.147
2.336
5.845
0.123
5.009
0.433
9.003
7579
23.960
0.543,
0.219
4.562330 Chap.9 Input Modeling
Lifetime, usually considered a continuous variable, is recorded here to three-
decimal-place accuracy. The histogram is prepared by placing the data in class
intervals. The range of the data is rather large, from 0.002 day to 144.695 days.
However, most of the values (30 of 50) are in the zero-to-S-day range. Using
intervals of width three results in Table 9.2. The data of Table 9.2 are then used
to prepare the histogram shown in Figure 9.3. <
‘Table 9.2. Electronic Chip Data
Chip Life
(Days) Frequency
Osx <3 23
3
oo. When using Newton’s method, B is approached
through increments of size f (B)-1)/f'(Bj-1)- Equation (9.12) is used to com-
pute f@;-1) and Equation (9. 15) is used to compute f'(B;-1) as follows:
a ™ xP enx;
p@ <2, 22 xPCnxi)? 4 Aten Xr en Xs) “ x; (Sta xPenxs)”
= (9.15)
B Die XP i xy342 Chap.9 Input Modeling
Equation (9.15) can be derived from Equation (9.12) by differentiating f(B)
with respect to 8. The iterative process continues until f @) = = 0, for example,
until fl < < 0.001.
Consider the data given in Example 9.: "3. These data concern the failure
of electronic components and may come from an exponential distribution. In
Example 9.11, the parameter 2 was estimated on the hypothesis that the data
were from an exponential distribution. If the hypothesis that the data came
from an exponential distribution is rejected, an alternative hypothesis is that the
data come from a Weibull distribution. The Weibull distribution is suspected,
since the data pertain to electronic component failures which occur suddenly.
Equation (9.14) is used to determine Bo. For the data in Example 9.3,
n = 50, ¥ = 11.894, X? = 141.467, and 72, X? = 37,575.850, so that S? is
found by Equation (9.2) to be
37,578.850 — 50(141.467)
2 = = 622.650
s a 622.6
and $ = 24,953. Thus,
> _ 11.894
= 0.477
bo = s4953 = °
To compute B; using Equation (9.13) requires the determination of f Bo) and
f'(Bo) using Equations (9.12) and (9.15). The following sacinonel values
are needed: 303, X# = 115.125, 2%, enX; = 38.294, 120, XP tnx, =
i=l
292,629, and 732, Theat 1057.781. Thus,
50(292.629)
f@o) = aa + 38.294 — Tee = 16.024
and
A =50 50(1057.781) , 50(292.629)?
"By = 250 _ 50(1057.781) , 50(292.629)" _
FO) = Care 1is125* ~(@15.125)2 356.110
‘Then, by Equation (9.13),
= 16.024
Ay = 0477 - Te = 0522
After four iterations, |f(B3)| < 0.001, at which point B = By = 0.525
is the approximate solution to Equation (9.10). Table 9.4 contains the values
needed to complete each iteration.
Now, & can be determined using Equation (9.11) with 6 = 8 = 0.525 as
follows:
= 6.227
30 oe
If By is sufficiently close to B, the procedure converges quickly, usually
in four to five iterations. However, if the procedure appears to be diverging,
[sey 525Sec. 9.4 Goodness-of-fit Tests 343
‘Table 9.4. Iterative Estimation of Parameters of the Weibull Distribution
ar rr) @ A ~
JB ExP Bxbenx, Exexye 58) £8) Bras
0 0477 115.125 292.629 1057.781 16.024 ~356.110 0.522
1 0.522 129,489 344.713 1254.11 1.008 -313.540 0.525
2 0.525 130.603 348.769 1269.547 0.004 -310.853 0.525
30.525 130.608 348.786 1269.614 0.000 -310.841 0.525
try other initial guesses for By — for example, one-half the initial estimate or
twice the initial estimate.
The difficult task of determining parameters for the Weibull distribution
by hand emphasizes the value of having software support for input model-
ing. <
9.4 Goodness-of-Fit Tests
Hypothesis testing was discussed in Section 7.4 with respect to testing random
numbers. In Section 7.4.1 the Kolmogorov-Smirnov test and the chi-square
test were introduced. These two tests are applied in this section to hypotheses
about distributional forms of input data.
Goodness-of-fit tests provide helpful guidance for evaluating the suit-
ability of a potential input model. However, since there is no single correct
distribution in a real application, you should not be a slave to the verdict of
such tests. It is especially important to understand the effect of sample size. If
very little data are available, then a goodness-of-fit test is unlikely to reject any
candidate distribution; but if a lot of data are available, then a goodness-of-fit
test will likely reject all candidate distributions. Therefore, failing to reject a
candidate distribution should be taken as one piece of evidence in favor of that
choice, while rejecting an input model is only one piece of evidence against the
choice.
9.4.1 Chi-Square Test
One procedure for testing the hypothesis that a random sample of size n of
the random variable X follows a specific distributional form is the chi-square
goodness-of-fit test. This test formalizes the intuitive idea of comparing the
histogram of the data to the shape of the candidate density or mass function.
The test is valid for large sample sizes, for both discrete and continuous distri-
butional assumptions, when parameters are estimated by maximum likelihood.
The test procedure begins by arranging the n observations into a set of k class
intervals or cells. The test statistic is given by
k
Oj; — E;)*
xe yp ey (9.16)
i
i=l344° Chap.9__ Input Modeling
where 0; is the observed frequency in the ith class interval and E; is the
expected frequency in that class interval. The expected frequency for each
class intervalis computed as E; = np; , where p; is the theoretical, hypothesized
probability associated with the ith class inferval.
Itcan be shown that x} approximately follows the chi-square distribution
with k — s ~ 1 degrees of freedom, where s represents the number of param-
eters of the hypothesized distribution estimated by the sample statistics. The
hypotheses are:
‘Ao: the random variable, X, conforms to the distributional assumption with
the parameter(s) given by the parameter estimate(s)
#4: the random variable X does not conform
‘The critical value eee is found in Table A.6. The null hypothesis, Ho, is
rejected if x2 > epee
In applying the test, if expected frequencies are too small, X@ will reflect
not only the departure of the observed from the expected frequency but the
smallness of the expected frequency as well. Although there is no general
agreement regarding the minimum size of E;, values of 3, 4, and 5 have been
widely used. In Section 7.4.1, when the chi-square test was discussed, the
minimum expected frequency of five was suggested. If an E; value is too small,
it can be combined with expected frequencies in adjacent class intervals, The
corresponding O; values should also be combined, and k should be reduced by
one for each cell that is combined.
If the distribution being tested is discrete, each value of the random vari-
able should be a class interval, unless it is necessary to combine adjacent class
intervals to meet the minimum expected cell-frequency requirement. For the
discrete case, if combining adjacent cells is not required,
Pi = pi) = P(X = x)
Otherwise, p; is determined by summing the probabilities of appropriate ad-
jacent cells,
If the distribution being tested is continuous, the class intervals are given
by [ai-1, a)), where a, and a; are the endpoints of the ith class interval. For
the continuous case with assumed pdf f(x), or assumed cdf F(x), p; can be
computed by:
a
Pi -[ fQ@)dx = Fla) — F@-)
ain
For the discrete case, the number of class intervals is determined by the
number of cells resulting after combining adjacent cells as necessary. How-
ever, for the continuous case the number of class intervals must be specified.
Although there are no general rules to be followed, the recommendations in
Table 9.5 are made to aid in determining the number of class intervals for con-
tinuous data.Sec. 9.4 Goodness-of-Fit Tests 345
Table 9.5, Recommendations for
Number of Class Intervals
for Continuous Data
‘Sample Size, Number of Class intervals,
n k
20 Do not use the chi-square test
50 5to10
100 10 to 20
>100 sfiton/s
EXAMPLE 9.13 (Chi-Square Test Applied to Poisson Assumption)
In Example 9.7, the vehicle-arrival data presented in Example 9.2 were ana-
lyzed. Since the histogram of the data, shown in Figure 9.2, appeared to follow
a Poisson distribution, the parameter, @ = 3.64, was determined. Thus, the
following hypotheses are formed:
Hp: the random variable is Poisson distributed
Hy: the random variable is not Poisson distributed
The pmf for the Poisson distribution was given in Equation (5.18) as follows:
Hey
>) = 0,1,2,... (0.17)
0, otherwise
For a = 3.64, the probabilities associated with various values of x are obtained
using Equation (9.17) with the following results:
(0) = 0.026, p(6) = 0.085
PCT) = 0.044
(8) = 0.020
(9) = 0.008
(10) = 0.003,
pG) = 0.140, p(11) = 0.001
With thisinformation, Table 9.6is constructed. The value of Ey is given by npy =
100(0.026) = 2.6. In a similar manner, the remaining E; values are determined.
Since E, = 2.6 < 5, E; and E2 are combined. In that case O; and O2 are
also combined and k is reduced by one. The last five class intervals are also
combined for the same reason, and k is further reduced by four. .
The calculated x? is 27.68. The degrees of freedom for the tabulated
value of x? is k—s -1=7~-1-1=5. Here, s = 1, since one parameter,
@, was estimated from the data.. At the 0.05 level of significance, the critical
value XGos,5 is 11.1. Thus, Ho would be rejected at level of significance 0.05.
The analyst may therefore want to search for a better-fitting model or use the
empirical distribution of the data. <346 ~~ Chap. 9" Input Modeling
Table 9.6. Chi-Square Goodness-of-Fit Test for Example 9.2
Xi Observed Frequency, 0; Expected Frequency, E; oer
0 2 26
: to | n oe | 122 | 787
2 19 174 oas
. 3 7 24 0.80
4 10 192 441
5 8 14.0 257
6 7 as 0.26
1 5 44
8 5 20
9 shir 08} 76 Jr
10 3 03
u ef 01 —
1 100.0 TB
9.42 Chi-Square Test with Equal Probabilities
Tf a continuous distributional assumption iy being tested, class intervals that
are equal in probability rathér than equal in width of interval should be used.
‘This has been recommended by a number of authors (Mann and Wald [1942];
Gumbel [1943]; Law and Kelton [2000]; Stuart, Ord, and Arnold [1998]). It
should be noted that the procedure is not applicable to data collected in class
intervals, where the raw data have been discarded ot lost.
Unfortunately, there is as yet no method for determining the probability
associated with each interval that maximizes the power for a test of a given size.
(The power of a test is defined as the probability of rejecting a false hypothesis.)
However, if using equal probabilities, then p; = 1/k. Since we recommend
E, = np, > 5
substituting for p; yields
"4s
z2
and solving for k yields
n
kez .18)
35 (9.18)
Equation (9.18) was used in determining the recommended maximum number
of class intervals in Table 9.5.
Ifthe assumed distribution is normal, exponential, or Weibull, the method
described in this section is straightforward. Example 9.14 indicates how the
Procedure is accomplished for the exponential distribution. If the assumedSec.9.4 Goodness-of-Fit Tests 347
distribution is gamma (but not Erlang), or certain other distributions, then
the computation of endpoints for class intervals is complex and may require
numerical integration of the density function. Statistical-analysis software is
very helpful in such cases.
EXAMPLE 9.14 (Chi-Square Test for Exponential Djstribution)
In Example 9.11, the failure data presented in Example 9.3 were analyzed.
Since the histogram of the data, shown in Figure 9.3, appeared to follow an
exponential distribution, the parameter X = 1/X =.0.084 was determined.
‘Thus, the following hypotheses are formed:
Ho: the random variable is exponentially distributed
Hy: the random variable is not exponentially distributed
In order to perform the chi-square test with intervals of equal probability, the
endpoints of the class intervals must be determined. Equation (9.18) indicates
that the number of intervals should be less than or equal to n/5. Here, n = 50, so
that k < 10. In Table 9:5, it is recommended that 7 to 10 class intervals be used.
Let k =8, then each interval will have probability p = 0.125. The endpoints for
each interval are computed from the cdf for the exponential distribution, given
in Equation (5.27), as follows:
F(a) = 1-6 (9.19)
where a; represents the endpoint of the ith interyal, i = 1,2,...,k. Since
F (ai) is the cumulative area from zero to a;, F(a:) = ip, so Equation (9.19)
can be written as .
ip =1-e%%
or
ei = 1—ip
Taking the logarithm of both sides and solving for a; gives a general result
for the endpoints of k equiprobable intervals for the exponential distribution,
namely
a = ~ dena ~ip, i =OL.0k (9.20)
Regardless of the value of A, Equation (9.20) will always result in ay = 0 and
= 00. With 4 = 0.084 and k = 8, a; is determined from Equation (9.20) as
1,
Toa" 0.125) = 1,590
Continued application of Equation (9.20) fori = 2,3, ..., 7resultsin az, ...,a7
as 3.425, 5.595, 8.252, 11.677, 16.503, and 24.755. Since k = 8, ag = oo. The
first interval is [0, 1.590), the second interval is [1.590, 3.425), and so on. The
expectation is that 0.125 of the observations will fall in each interval. The OA
servations, expectations, and the contributions to pe calculated value of x
are shown in Table 9.7. The calculated value of x? is 39.6. The degrees of
aq=-348 Chap.9_ Input Modeling
Table 9.7. Chi-Square Goodness-of-Fit Test for Example 9.14
Class Interval Observed Frequency, 0; Expected Frequency, E; eae
{0, 1.590) 19 6.25 26.01
{[1.590, 3.425) 10 6.25 225
(3.425, 5.595) 3 625 O81
15,595, 8.252) 6 625 001
(8.252, 11 677) 1 6.25 441
(11.677, 16.503) 1 6.25
(16.503, 24.755) 4 6.25
(24.755, 00)- 6 6.25
50 50
freedom are given by ks ~1=8~1-1=6. Ata =0.05, the tabulated
value of 956 i8 12.6. Since x3 > x2os ,, the null hypothesis is rejected. (The
value of X01, i8 168, o the mull hypothesis would also be rejected at level of
significance a = 0.01.) <
9.4.3 Kolmogorov-Smirnov Goodness-of-Fit Test
The chi-square goodness-of-fit test can accommodate the estimation of param-
eters from the data with a resultant decrease in the degrees of freedom (one for
cach parameter estimated). The chi-square test requires that the data be placed
in class intervals, and in the case ofa continuous distributional assumption, this
srouping is arbitrary. Changing the number of classes and the interval width
affects the value of the calculated and tabulated chi-square. A hypothesis may
be accepted when the data are grouped one way but rejected when grouped
another way. Also, the distribution of the chi-square test statistic is known only
approximately, and the power of the test is sometinies rather low. As a result
of these considerations, goodness-of-fit tests, other than the chi-square, are de-
sired. The Kolmogorov-Smirnov test formalizes the idea behind examining a
@-q plot.
The Kolmogoroy-Smirnov test was presented in Section 7.4.1 to test for
the uniformity of numbers and again in Section 7.4.4 to perform the gap test.
Both of these uses fall into the category of testing for goodness-of-fit. Any
Continuous distributional assumption can be tested for goodness-of-fit using
the method of Section 7.4.1, while discrete distributional assumptions can be
tested using the method of Section 7.
The Kolmogorov-Smirnov test is particularly useful when sample sizes
are small and when no parameters have been estimated from the data. When
Parameter estimates have been made, the critical values in Table A.8 are biased;
in particular, they are too conservative. In this context “conservative” means
that the critical values will be too large, resulting in smaller ‘Type I (a) errorsSec. 9.4 Goodness-of Fit Tests 349
than those specified. The exact value of & can be determined in some instances
as discussed at the end of this section.
The Kolmogorov-Smirnov test does not take any special tables when an
exponential distribution is assumed. The following example indicates how the
test is applied in this instance. (Notice that it is not necessary to estimate the
Parameter of the distribution in this example, permitting the use of Table A.8.)
EXAMPLE 9.15 (Kolmogoroy-Smirnoy Test for
Exponential Distribution)
Suppose that 50 interarrival times (in minutes) are collected over the following
10-minute interval (arranged in order of occurrence):
oad” 053° 204 2.74 2.00 030 254 0.52 2.02 if 153 021
280 0.04 135 832 234 1.95 0.10 142 046 0.07 1.09 0.76
555 3.93 1.07 2.26 288 0.67 1.12 026 457 537 012 3.19
163 146 1.08 2.06 085 083 244 2:11 3.15 290 658 0.64
‘The null hypothesis and its alternate are formed as follows:
Ho: the interarrival times are exponentially distributed
Aj: the interarrival times are not exponentially distributed
The data were collected over the interval 0 to T = 100 minutes. It can be
shown that if the underlying distribution of interarrival times {T;, Ts, ..) is
exponential, the arrival times are uniformly distributed on the interval (0, 7).
‘The arrival times 71, T, + Tp, Ty + T + Ts,..., Ty + +--+ Too are obtained
by adding interarrival times. The arrival times are then normalized to a (0, 1)
interval so that the Kolmogorov-Smirnov test, as presented in Section 7.4.1, can
be applied. Ona (0, 1) interval, the points will be [71/T,, (T,+75)/T, ..., (Ti+
++: +T50)/T}. The resulting 50 data points are as follows:
0.0044 0.0097 0.0301 0.0575 0.0775 0.0805 0.1059 0.1111 0.1313 01502
0.1655 0.1676 0.1956 0.1960 0.2095 0.2927 03161 0.3356 0.3366 03508
0.3553 0.3561 0.3670 0.3746 0.4300 0.4694 0.4796 0.5027 0.5315 0.5382
0.5494 0.5520 0.5977 0.6514 0.6526 0.6845 0.7008 0.7154 0.7262 0.7468
0.7553 0.7636 0.7880 0.7982 0.8206 0.8417 08732 0.9022 0.9680 0.9744
Following the procedure in Example 7.6 yields a D+ of 0.1054 and a D~ of
0.0080. Therefore, the Kolmogorov-Smirnov statistic is D = max(0.1054,
0.0080) = 0.1054. The critical value of D obtained from Table A.8 for a
level of significance of a = 0.05 and n = 50 is Doos = 1.36//n = 0.1923.
Since D = 0.1054, the hypothesis that the interarrival times are exponentially
ributed cannot be rejected. <Y
350 Chap.9 Input Modeling
‘The Kolmogorov-Smirnov test has been modified so that it can be used in
several situations where the parameters are estimated from the data. The com-
putation of the test statistic is the same, but different tables of critical values
are used. Different tables of critical values-are required for different distribu-
tional assumptions. Lilliefors [1967] developed a test for normality. The null
hypothesis states that the population is one of the family of normal distributions
without specifying the parameters of the distribution. The interested reader
may wish to study Lilliefors’ original work, as he describes how simulation was
used to develop the critical values.
Lilliefors [1969] also modified the critical values of the Kolmogorov-
Smirnov test for the exponential distribution. Lilliefors again used random
sampling to obtain approximate critical values, but Durbin [1975] subsequently
obtained the exact distribution. Conover [1980] gives examples of Kolmogo-
rov-Smirnov tests for the normal and exponential distributions. He also refers
to several other Kolmogorov-Smirnov-type tests which may be of interest to
the reader. ,
A test that is similar in spirit to the Kolmogorov-Smirnov test is the
Anderson-Darling test. Like the Kolmogorov-Smirnov test, the Anderson.
Darling test is based on the difference between the empirical cdf and the fitted
cdf; unlike the Kolmogorov-Smirnov test, the Anderson-Darling test is based
on a more comprehensive measure of difference (not just the maximum dif-
ference) and is more sensitive to discrepancies in the tails of the distributions.
The critical values for the Anderson-Darling test also depend on the candidate
distribution and on whether or not parameters have been estimated. Fortu-
nately, this test and the Kolmogorov-Smirnoy test have been implemented in a
number of software packages that support simulation input modeling.
9.4.4 p-Values and “Best Fits”
To apply a goodness-of-fit test a significance level must be chosen. Recall
that the significance level is the probability of falsely rejecting Hp: the random
variable conforms to the distributional assumption. The traditional ‘Significance
levels are 0.1, 0.05, and 0.01. Prior to the availability of high-speed computing,
having a small set of standard values made it possible to produce tables of
useful critical values. Now most statistical software computes critical values as
needed, rather than storing them in tables. Thus, if the analyst prefers a level
of significance of, say, 0.07, then he or she can choose it.
However, rather than require a prespecified significance level, many soft-
ware packages compute a p-value for the test statistic. The p-value is the
significance level at which one would just reject Ho for the given value of the
test statistic. Therefore, a large p-value tends to indicate a good fit (we would
have to accept a large chance of error in order to reject), while a small p-value
Suggests a poor fit (to accept we would have to insist on almost no risk).
Recall Example 9.13, in which a chi-square test was used to check the
Poisson assumption for the vehicle arrival data. The value of the test statisticSec. 9.5. Selecting Input Models without Data 351
was x = 27.68 with 5 degrees of freedom. The p-value for this test statistic is
0.00004, meaning that we would reject the hypothesis that the data are Poisson
at the 0.00004 significance level (recall that we rejected the hypothesis at the
0.05 level; now we know that we would also Teject it at even lower levels).
The p-value can be viewed as a measure of fit, with larger values being
better. This suggests that we could fit every distribution at our disposal, com.
pute a test statistic for each fit, and then choose the distribution that yields
the largest p-value. While we know of no input modeling software that im-
plements this specific algorithm, many such packages do include a “best-fit”
option in which the software recommends an input model to the user based on
evaluating all feasible models. While the software may also take into account
other factors — such as whether the data are discrete or continuous, bounded
or unbounded — in the end some summary measure of fit, ike the p-value, is
used to rank the distributions. There is nothing wrong with this, but there are
several things to keep in mind:
1. The software may know nothing about the physical basis of the data, and
that information can suggest distribution families that are appropriate
(see the list in Section 9.2.2). Remember that the goal of input modeling
is often to fill in gaps or smooth the data, rather than find an input model
that conforms as closely as possible to the given sample.
2. Recall that both the Erlang and the exponential distributions are special
Cases of the gamma, while the exponential is also a special case of the more
flexible Weibull. Automated best-fit procedures tend to choose the more
flexible distributions (gamma and Weibull over Erlang and exponential)
because the extra flexibility allows closer conformance to the data anda
better summary measure of fit. But again, close conformance to the data
may not always lead to the most appropriate input model.
3. A summary statistic, like the p-value, is just that, a summary measure. It
says little or nothing about where the lack of fit occurs (in the body of the
distribution, in the right tail or in the left tail). A human, using graphical
tools, can see where the lack of fit occurs and decide whether or not it is
important for the application at hand.
Our recommendation is that automated distribution selection be used as one of
several ways to suggest candidate distributions. Always inspect the automatic
Selection using graphical methods, and remember that the final choice is yours.
9.5 Selecting Input Models without Data
Unfortunately, it is often necessary in practice to develop a simulation model
—perhaps for demonstration pirposes or a preliminary study—before any
Process data are available. In this case the modeler must be resourceful in
choosing input models and must carefully check the sensitivity of results to the
chosen models.352 Chap.9_ Input Modeling
‘There are a number of ways to obtain information about a process even
if data are not available:
Engineering data Often a product or process has performance ratings pro-
vided by the manufacturer (for example; the mean time to failure of a disk
drive is $000 hours; a laser printer can produce 4 pages/minute; the cutting
speed of a tool is 1 inch/second, etc.). Company rules may specify time
or. production standards. These values provide a starting point for input
modeling by fixing a central value.
Expert option Talk to people who are experienced with the process or similar
processes, Often they can provide optimistic, pessimistic, and most likely
times. They may also be able to say if the process is nearly constant or highly
variable, and they may be able to define the source of variability.
Physical or conventional limitations Most real processes have physical limits
on performance (for example, computer data entry cannot be faster than a
person can type). Because of company policies, there may be upper limits
‘on how long a process may take. Do not ignore obvious limits or bounds
that narrow the range of the input process.
The nature of the process The description of the distributions in Section 9.2.2
can be used to justify a particular choice even when no data are available.
‘When data are not available, the uniform, triangular, and beta distributions are
often used as input models. The uniform can be a poor choice, because the
upper and lower bounds are rarely just as likely as the central values in real
processes. If, in addition to upper and lower bounds, a most-likely value can be
given, then the triangular distribution can be used. The triangular distribution
places much of its probability near the most-likely value and much less near
the extremes (see Section 5.4). If a beta distribution is used, then be sure to
plot the density function of the selected distribution, since the beta can take
unusual shapes.
A useful refinement is obtained when a minimum, maximum, and one or
more “breakpoints” can be given. A breakpoint is an intermediate value and
a probability of being less than or equal to that value. The following example
illustrates how breakpoints are used.
EXamPLe 9.16
For a production planning simulation, the sales volume of various products is
required. The salesperson responsible for product XYZ-123 says that no fewer
than 1000 units will be sold because of existing contracts, and no more than 5000
units will be sold because that is the entire market for the product. Based on
her experience, she believes that there is a 90% chance of selling more than
2000 units, a 25% chance of selling more than 3500 units, and only a 1% chance
of selling more than 4500 units.
Table 9.8 summarizes this information. Notice that the chances of exceed-
ing certain sales goals have been translated into the cumulative probability of
being less than or equal to those goals. With the information in this form the
method of Section 8.1.5 can be employed to generate simulation inputdata.Sec. 9.6 Multivariate and Time-Series Input Models 353
Table 9.8. Summary of Sales Information
Interval ‘Cumulative
i (Hours) Frequency, Ci
1 1000