9 Input Modeling

Input models provide the driving force for a simulation. In the simulation ofa queueing system, typica input models are the distributions of time between arrivals and of service times. For a supply-chain simulation, input models include the distributions of demand and of lead time. For the simulation of a reliability problem, the distribution of time to failure of a component is an example of an input model.

In real-world simulation applications, choosing 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 interpretation could 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 a substantial 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.

2. Identify a probability distribution to represent the input process. When data are available, this step typically begins with the development of a histogram of the data. Given the histogram and a structural knowledge of the process, a family of distributions is chosen. Fortunately,

several well-known distributions often provide good approximations in practice.

3. Choose parameters that determine a specific instance of the distribution family. When data are available, these parameters may be estimated from the data.

4. Evaluate the chosen distribution and the associated parameters for goodness of fit. Goodness of fit may be evaluated informally, via graphical methods, or formally, via statistical tests. The chi-square and the Kolmogorov–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 distributions, and repeats the procedure. If several iterations of this procedure fail to yield a fit between an assumed distributiona form and the collected data, the empirical form of the distribution may be used.

Each of these steps is discussed in this chapter. Although software is now widely available to accomplish Steps 2, 3, and 4—including such stand-alone programs as ExpertFit® and $\mathrm{S t a t} {: :} \mathrm{F i t} ^ {\circledR}$ and such integrated programs 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.

1 Data Collection

Problems are found at the end of each chapter, as exercises for the reader, in textbooks about mathematics, physics, chemistry, and other technical subjects. Years and years of working these problems could 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 and it is one of the most important and difficult problems in simulation.

The focus of input modeling is on the statistical aspects of fitting probability distributions to data, distributions that will later provide the driving inputs to the simulation. However, before “fitting” can occur there needs to be accurate and relevant input data available (or collectible), and the properties o the data need to be well understood. Too often data are handed to input modeling software, which wil always fit some probability model to the data, whether it is justified or not. The following examples illustrate common pitfalls of this approach:

Example 1: Stale Data

A simulation study will investigate ways to reduce the time patients spend in a cancer screening clinic. One part of the screening protocol involves a nurse obtaining the patient’s history. In a 2002 study of nursing activities, observations of nurses taking patient histories were collected by a consulting firm. The hospital still has this data. Should it be used for input modeling in the current project?

In this case the answer is probably no, because deeper investigation by the engineers on the current project discovered that in 2008 the nurses stopped recording screening notes on a paper char and began entering them directly into a laptop computer. Therefore, the 2002 data are of a different history-taking process and may significantly misrepresent the current process. If current data canno

be collected, either due to the time required or the expense, then it would be better to use estimates provided by the nurses themselves rather than rely on the out-of-date study. Section 6 describes methods for converting expert opinion into probability distributions. The engineers might then check the nurses’ estimates by volunteering to have their own histories taken and timing the process to get a few actual data points.

Example 2: Unexpected Data

For security purposes people who enter a government building are required to pass through a security screening that includes a metal detector. Before new security protocols are implemented a simulation study will be conducted to assess their impact. To facilitate the study, 1000 observations of the time, in minutes, required to pass through the metal detector station were collected. The average value o these data was about 1/2 minute, with standard deviation almost the same. This suggests that a good input model for the time to pass through the metal detector is the exponential distribution, since its mean and standard deviation are always the same and exponential distributions are often used to model service times in queueing systems.

However, the analyst is uncomfortable with this input model because in her experience it only takes a few seconds to pass through a metal detector. She creates a histogram of the data, shown in Figure 1, which makes it immediately apparent that the time to pass through metal detection consists of two distinct distributions. Thinking carefully, she realizes that some people pass through the detector without incident, but others trigger the alarm and have to remove remaining meta

Figure 1 Histogram of 1000 observations of time through the metal detector measured in minutes.

objects and try again, in some cases even requiring a personal inspection. Therefore, three distinc input models are required for this data: one for the time to pass through without setting off the detector (the 764 left-most observations in the histogram); one for the time to pass through when the alarm does go off (the 236 right-most observations in the histogram); and a model for the chance tha a person sets off the alarm. While an exponential distribution could certainly match the summary statistics from the data (mean and standard deviation). the fitted exponential distribution would have missed a key feature of security screening, almost certainly leading to an inaccurate simulation.

Unfortunately, inspecting the data carefully the analyst spots an additional complication: The histogram in Figure 1 shows negative values in the sample (nine total). A negative time to pass through metal detection is clearly impossible, so what should the analyst do? It may be that the observations were recorded correctly but a mistake was made in calculating the times; if that is the case then the analyst may be able to reconstruct the correct values. If not, then the nine negative values should be discarded before fitting the distribution. But should they also be discarded when estimating the probability that someone sets off the alarm? That is, should the estimate be 236/1000 or 236/991? If the analyst can convince herself that the 9 negative values were actually small times that were misrecorded, then she might want to retain them for estimating this probability. If she is unsure, then the simulation should be run with both estimates to see if there is a substantial impact.

Example 3: Time-varying Data

A computer manufacturer provides a customer support center where owners of the product may cal with questions or concerns. The call center is staffed from 8A.M. to 8 P.M. Eastern Time.A simulation study is planned to help set staffing levels, since staff salaries and benefits are the primary costs. The call center has extensive records on daily call volume showing that, on average, they receive 12,277 calls per day. The number of random events in a fixed amount of time or space is often well modeled by a Poisson distributions, so a Poisson distribution is fit to the data (using methods described later in the chapter) and appears to fit well. When the number of arrivals is Poisson, then the time between arrivals—which is what we need to simulate calls—must be exponentially distributed. The natura estimate of the arrival rate for this exponential distribution is 12,277 arrivals per 12-hour day, bu since a “day” is a pretty large time unit, it is converted to 17 calls per minute.

The analysis in the previous paragraph seems to be sound, and the choice ofa Poisson distribution (implying exponential time between arrivals) was even supported by data. However, the level of data aggregation, a 12-hour day, masked a very important feature of call centers: The arrival rate of calls can, and typically does, vary substantially by time of day, day of week, and even time of the year (e.g., right after the holiday gift-giving season). We often assume that the process of interest is stationary, meaning that its basic properties do not change over time. However, if trends exist then accounting for them can be far more important than the choice of probability distribution to represent the remaining uncertainty in the system. In Section 5 we describe a method for fitting a nonstationary Poisson process, which is likely more appropriate for this call center.

Example 4: Dependent Data

Monthly demands for a particular brand of sports drink by a large retailer have been extracted from an orders database to use in a supply chain simulation. A simple input model for this data is that monthly demands are independent and identically distributed (i.i.d.), and this chapter describes methods for fitting a distribution to such data and evaluating the fit. But it is important to evaluate the “i.i.d.” assumption also.

Figure 2 Time-lagged demand data.

One informal way to check that these data are identically distributed is to plot the demands in time order to look for trends (e.g., demands are decreasing over time) or seasonality (there is a regular cycle of demand, for instance, high demand each summer followed by low demand each winter). Data with trends or seasonality should typically be “de-trended” first (usually by fitting a model to the trend or seasonal pattern) and then a distribution can be fit to any remaining variation.

Independence can be checked via time-lagged scatter plots. Figure 2 is a scatter plot of each month’s demand for the sports drink plotted against the following month’s demand (e.g., January’s demand vs. February’s demand). Clearly there is a strong negative correlation, which means large demands in one month tend to be followed by low demands in the next month, and vice versa. This feature might be the result of the retailer overstocking in one month, and then ordering less the following month so as to work off the excess inventory. Time lags greater than one month can also be examined. By blindly making the standard assumption of i.i.d. inputs this important phenomenon would be missed, leading to less effective supply chain management. In Section 7 we present methods for modeling dependent input data.

As these examples illustrate, simply having data is not enough to support effective input modeling. Data can be out of date or “dirty” (containing errors). Sometimes the effort or cost to transform data into a usable form, or “clean” it of errors, can be as significant as that required to obtain it. Inpu modeling nearly always requires the analyst to use their judgment as well as to apply appropriate statistical tools. Understanding which input models are the most and least reliable is important for judging the reliability of the conclusions drawn from the simulation study. And since uncertainty about the correctness of the input models can never be entirely eliminated, it is sensible to run the simulation with several plausible input models to see if the conclusions are robust or highly sensitive to the choices.

Many lessons can be learned from an actual experience at data collection. The first five exercises at the end of this chapter suggest some situations in which the student can gain such experience.

The following suggestions might enhance and facilitate data collection, although they are no 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 be 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 ensure 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.

2. Try to analyze the data as they are being collected. Figure out whether the data being collected are adequate to provide the distributions needed as input to the simulation. Find out whether any data being collected are useless to the simulation. There is no need to collect superfluous data.

3. 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 P.M. to 3:00 P.M. and 3:00 P.M. to 4:00 P.M., and check to see whether the data are homogeneous for 2:00 P.M. to 3:00 P.M. on Thursday and Friday. When checking for homogeneity, an initial test is to see whether 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 test of the equivalence of the distributions, perhaps via a quantile-quantile plot (described later).

4. Be aware of the possibility of data censoring, in which a quantity of interest 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, treat a patient, or have a component 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.

5. To discover whether there is a relationship between two variables, build a scatter diagram. Sometimes a look at the scatter diagram will indicate whether there is a relationship between the variables of interest.

6. Consider the possibility that a sequence of observations that appear to be independent actually has autocorrelation. Autocorrelation can exist in successive time periods or for successive customers.

7. Keep in mind the difference between input data and output or performance data, and be sure to collect input data. Input data typically represent 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 system when

subjected to the inputs, performance that we might be trying to improve. In a queueing simulation, the customer arrival times are usually inputs, whereas the customer delay is an output. Performance data are useful for model validation.

Again, these arejust a few suggestions.As a rule, data collection and analysis must be approached with great care.

2 Identifying the Distribution with Data

In this section, we discuss methods for selecting families ofinput distributions when data are available and are believed to be independent and identically distributed. The specific distribution within a family is specified by estimating its parameters, as described in Section 3. Section 6 takes up the case in which data are unavailable.

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 frequencies are adjusted.)

2. Label the horizontal axis to conform to the intervals selected.

3. Find 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 on the amount o scatter or dispersion in the data. Hines et al. [2002] 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 ragged, coarse, and appropriate histograms of the same data are shown in Figure 3. 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 theoretica distribution. If continuous, a line drawn through the center point of each class interval frequency should result in a shape like that of a 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 could be necessary to combine adjacent cells to eliminate the ragged appearance of the histogram. If the histogram is associated with discrete data, it should look like a probability mass function.

Figure 3 Ragged, coarse, and appropriate histograms: (a) original data—too ragged; (b) combining adjacent cells—too coarse; (c) combining adjacent cells—appropriate.

Example 5: Discrete Data

The number of vehicles arriving at the northwest corner of an intersection in a 5-minute period between 7:00 A.M. and 7:05 A.M. was monitored for five workdays over a 20-week period. Table 1 shows the resulting data. The first entry in the table indicates that there were 125-minute periods during which zero vehicles arrived, 10 periods during which one vehicle arrived, and so on.

Table 1 Number of Arrivals in a 5-Minute Period

Figure 4 Histogram of number of arrivals per period.

The number of automobiles is a discrete variable, and there are ample data, so the histogram may have a cell for each possible value in the range of the data. The resulting histogram is shown in Figure 4.

Example 6: Continuous Data

Life tests were performed on a random sample of electronic components at 1.5 times the nomina voltage, and their lifetime (or time to failure), in days, was recorded:

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-5-day range. Using intervals of width three results in Table 2. The data of Table 2 are then used to prepare the histogram shown in Figure 5.

Figure 5 Histogram of component life.

2.2 Selecting the Family of Distributions

The purpose of preparing a histogram is to infer a known pdf or pmf. A family of distributions is selected on the basis of what might arise in the context being investigated along with the shape of the histogram. Thus, if interarrival-time data have been collected, and the histogram has a shape similar to an exponential distribution. Similarly, if measurements of the weights of pallets of freigh are being made, and the histogram appears symmetric about the mean, the assumption of a norma distribution would be warranted.

Table 2 Electronic Component Data

The exponential, normal, and Poisson distributions are frequently encountered and are not difficult to analyze from a computational standpoint. Although more difficult to analyze, the beta, gamma, and Weibull distributions provide a wide array of shapes and should not be overlooked during modeling of an underlying probabilistic process. Perhaps an exponential distribution was assumed, but it was found not to fit the data. The next step would be to examine where the lack of fit occurred. If the lack of fit was in one of the tails of the distribution, perhaps a gamma or Weibull distribution would fit the data more adequately.

There are literally hundreds of probability distributions that have been created; many were created with some specific physical process in mind. One aid to selecting distributions is to use the physical basis of the distributions as a guide. Here are some examples:

Binomial: Models the number of successes in n trials, when the trials are independent with common success probability p; for example, the number of defective computer chips found in a lot of n chips.

Negative Binomial (includes the geometric distribution): Models the number of trials required to achieve k successes; for example, the number of computer chips that we must inspect to find 4 defective chips.

Poisson: Models the number of independent events that occur in a fixed amount of time or space; for example, the number of customers that arrive to a store during 1 hour, or the number of defects found in 30 square meters of sheet metal.

Normal: Models the distribution of a process that can be thought of as the sum of a number of component processes; for example, a time to assemble a product that is the sum of the times required for each assembly operation. Notice that the normal distribution admits negative values, which could be impossible for process times.

Lognormal: Models the distribution of a process that can be thought of as the product of (meaning to multiply together) a number of component processes; for example, the rate of return on an investment, when interest is compounded, is the product of the returns for a number of periods.

Exponential: Models the time between independent events, or a process time that is memoryless (knowing how much time has passed gives no information about how much additional time wil pass before the process is complete); for example, the times between the arrivals from a large population of potential customers who act independently of one another. The exponential is a highly variable distribution; it is sometimes overused, because it often leads to mathematically tractable models. Recall that, if the time between events is exponentially distributed, then the number of events in a fixed period of time is Poisson.

Gamma: An extremely flexible distribution used to model nonnegative random variables. The gamma can be shifted away from 0 by adding a constant.

Beta: An extremely flexible distribution used to model bounded (fixed upper and lower limits) random variables. The beta can be shifted away from 0 by adding a constant and can be given a range larger than [0, 1] by multiplying by a constant.

Erlang: Models processes that can be viewed as the sum of several exponentially distributed processes; for example, a computer network fails when a computer and two backup computers fail, and each has a time to failure that is exponentially distributed. The Erlang is a special case of the gamma.

Weibull: Models the time to failure for components; for example, the time to failure for a disk drive. The exponential is a special case of the Weibull.

Discrete or Continuous Uniform: Models complete uncertainty: All outcomes are equally likely. This distribution often is used inappropriately when there are no data.

Triangular: Models a process for which only the minimum, most likely, and maximum values of the distribution are known; for example, the minimum, most likely, and maximum time required to test a product. This model is often a marked improvement over a uniform distribution.

Empirical: Resamples from the actual data collected; often used when no theoretical distribution seems appropriate.

Do not ignore physical characteristics of the process when selecting distributions. Is the process naturally discrete or continuous valued? Is it bounded or is there no natural bound? This knowledge, which does not depend on data, can help narrow the family of distributions from which to choose. And keep in mind that there is no “true” distribution for any stochastic input process. An input mode is an approximation of reality, so the goal is to obtain an approximation that yields useful results from the simulation experiment.

The reader is encouraged to complete Exercises 6 through 11 to learn more about the shapes of the distributions mentioned in this section. Examining the variations in shape as the parameters change is very instructive.

2.3 Quantile-Quantile Plots

The construction of histograms, as discussed in Section 2.1, and the recognition of a distributiona shape, as discussed in Section 2.2, are necessary ingredients for selecting a family of distributions to represent a sample of data. However, a histogram is not as useful for evaluating the $\mathscr {f} t$ of the chosen distribution. When there is a small number of data points, say, 30 or fewer, a histogram can be rather ragged. Further, our perception of the fit depends on the widths of the histogram intervals. But, even if the intervals are well chosen, grouping data into cells makes it difficult to compare a histogram to a continuous probability density function. A quantile-quantile $( q \ - q )$ plot is a useful tool for evaluating distribution fit, one that does not suffer from these problems.

If X is a random variable with cdf F, then the q-quantile of X is that value $\gamma$ such that $F ( \gamma ) =$$P ( X \leq \gamma ) = q$, for $0 < q < 1$. When F has an inverse, we write $\gamma = F^{- 1} ( q )$ Now let $\{x_{i} , i = 1 , 2 , \dots , n \}$ be a sample of data from X. Order the observations from the smallest to the largest, and denote these as $\{y_{j} , j = 1 , 2 , . . . n \}$, where $y_{1} \leq y_{2} \leq \cdot \cdot \cdot \leq y_{n}$. Let j denote the ranking or order number. Therefore,$j = 1$ for the smallest and $j = n$ for the largest. The $q {-} q$ plot is based on the fact that $y_{j}$ is an estimate of the $( j - 1 / 2 ) / n$ quantile of X. In other words,

Now, suppose that we have chosen a distribution with cdf $F$ as a possible representation of the distribution of X. If $F$ is a member of an appropriate family of distributions, then a plot of $y_{j}$ versus $F^{- 1} ( ( j - 1 / 2 ) / n )$) will be approximately a straight line. If F is from an appropriate family of distributions and also has appropriate parameter values, then the line will have slope 1. On the other hand, if the assumed distribution is inappropriate, the points will deviate from a straight line,

usually in a systematic manner. The decision about whether to reject some hypothesized model is subjective.

Example 7: Normal Q-Q Plot

A robot is used to install the doors on automobiles along an assembly line. It was thought that the installation times followed a normal distribution. The robot is capable of measuring installation times accurately. A sample of 20 installation times was automatically taken by the robot, with the following results, where the values are in seconds:

The sample mean is 99.99 seconds, and the sample variance is (0.2832)^2 seconds^2. These values can serve as the parameter estimates for the mean and variance of the normal distribution. The observations are now ordered from smallest to largest as follows:

The ordered observations are then plotted versus $F^{- 1} ( ( j - 1 / 2 ) / 20 )$, for $j = 1 , 2 , \dots , 20$, where F is the cdf of the normal distribution with mean 99.99 and variance (0.2832)^2, to obtain a q-q plot. The plotted values are shown in Figure 6, along with a histogram of the data. Notice that it is difficul to tell whether the data are well represented by a normal distribution from looking at the histogram, but the general perception of a straight line is quite clear in the q-q plot and supports the hypothesis of a normal distribution.

In the evaluation of the linearity of a q-q plot, the following should be considered:

a. The observed values will never fall exactly on a straight line.

b. The ordered values are not independent; they have been ranked. Hence, if one point is above a straight line, it is likely that the next point will also lie above the line. And it is unlikely that the points will be evenly scattered about the line.

c. The variances of the extremes (largest and smallest values) are much higher than the variances in the middle of the plot. Greater discrepancies can be accepted at the extremes. The linearity of the points in the middle of the plot is more important than the linearity at the extremes.

Figure 6 A q-q plot and histogram of the installation times.

Modern data-analysis software often includes tools for generating q-q plots, especially for the normal distribution. The $q {-} q$ plot can also be used to compare two samples of data to see whether they can be represented by the same distribution (that is, that they are homogeneous).$\operatorname{I f} x_{1} , x_{2} , \ldots , x_{n}$ are a sample of the random variable $X ,$, and $z_{1} , z_{2} , \ldots , z_{n}$ are a sample of the random variable $Z ,$ then plotting the ordered values of X versus the ordered values of Z will reveal approximately a straigh line if both samples are well represented by the same distribution (Chambers, Cleveland, and Tukey [1983]).

3 Parameter Estimation

After a family of distributions has been selected, the next step is to estimate the parameters o the distribution. Estimators for many useful distributions are described in this section. In addition, many software packages—some of them integrated into simulation languages—are now available to compute these estimates.

3.1 Preliminary Statistics: Sample Mean and Sample Variance

In a number of instances, the sample mean, or the sample mean and sample variance, are used to estimate the parameters of a hypothesized distribution; see Example 7. In the following paragraphs, three sets of equations are given for computing the sample mean and sample variance. Equations (1) and (2) can be used when discrete or continuous raw data are available. Equations (3) and (4) are used when the data are discrete and have been grouped in a frequency distribution. Equations (5) and (6) are used when the data are discrete or continuous and have been placed in class intervals. Equations (5) and (6) are approximations and should be used only when the raw data are unavailable.

If the observations in a sample of size n are $X_{1} , X_{2} , \ldots , X_{n}$, the sample mean $\bar {X}$ is defined by

and the sample variance $S^{2}$ is defined by

If the data are discrete and have been grouped in a frequency distribution, Equations (1) and (2) can be modified to provide for much greater computational efficiency. The sample mean can be computed as

and the sample variance as

where k is the number of distinct values of X and $f_{j}$ is the observed frequency of the value $X_{j}$ of X.

Example 8: Grouped Data

The data in Table 1 can be analyzed to obtain $n = 100 , f_{1} = 12 , X_{1} = 0 , f_{2} = 10 , X_{2} = 1 , \ldots ,$$\textstyle \sum_{j = 1} ^ {k} f_{j} X_{j} = 364$, and $\begin{array} {r} {\sum_{j = 1} ^ {k} = f_{j} X_{j} ^ {2} = 2080} \end{array}$. From Equation (3),

and, from Equation (4),

The sample standard deviation S is just the square root of the sample variance. In this case,${\boldsymbol {S}} =$$\sqrt {7 . 63} = 2 . 76$. Equations (1) and (2) would have yielded exactly the same results for $\bar {X}$ and $S^{2}$ It is preferable to use the raw data, if possible, when the values are continuous. However, data sometimes are received after having been placed in class intervals. Then it is no longer possible to obtain the exact sample mean and variance. In such cases, the sample mean and sample variance are approximated from the following equations:

and

where $f_{j}$ is the observed frequency in the jth class interval,$m_{j}$ is the midpoint of the jth interval, and c is the number of class intervals.

Example 9: Continuous Data in Class Intervals

Assume that the raw data on component life shown in Example 6 either was discarded or was lost. However, the data shown in Table 2 are still available. To approximate values for X¯ and $S^{2}$, Equations (5) and (6) are used. The following values are created:$f_{1} = 23 , m_{1} = 1 . 5 , f_{2} = 10 \ / \mathrm{\Omega}$$m_{2} = 4 . 5 , \ldots , \sum_{j = 1} ^ {49} f_{j} m_{j} = 614$ and $\textstyle \sum_{j = 1} ^ {49} f_{j} m_{j} ^ {2} = 37 , 226 . 5$. With $n = 50$, X is approximated from Equation (5) as

Then,$S^{2}$ is approximated from Equation (6) as

and

Applying Equations (1) and (2) to the original data in Example 6 results in $\bar {X} \ = \ 11 . 894$ and $S = 24 . 953$. Thus, when the raw data are either discarded or lost, inaccuracies could result

Table 3 Suggested Estimators for Distributions Often Used in Simulation

3.2 Suggested Estimators

Numerical estimates of the distribution parameters are needed to reduce the family of distributions to a specific distribution and to test the resulting hypothesis. Table 3 contains suggested estimators for distributions often used in simulation. Except for an adjustment to remove bias in the estimate of $\sigma^{2}$ for the normal distribution, these estimators are the maximum-likelihood estimators based on the raw data. (If the data are in class intervals, these estimators must be modified.) The reader is referred to Fishman [1973] and Law [2007] for parameter estimates for the uniform, binomial, and negative binomial distributions. The triangular distribution is usually employed when no data are available, with the parameters obtained from educated guesses for the minimum, most likely, and maximum possible values; the uniform distribution may also be used in this way if only minimum and maximum values are available.

Examples of the use of the estimators are given in the following paragraphs. The reader should keep in mind that a parameter is an unknown constant, but the estimator is a statistic (or random variable), because it depends on the sample values. To distinguish the two clearly here, if, say, a parameter is denoted by α, the estimator will be denoted by α.

Example 10: Exponential Distribution

We will use the exponential distribution to illustrate the maximum likelihood method for obtaining parameter estimators. The density function of the exponential distribution with rate parameter λ is

for $x \ge 0$. If $X_{1} , X_{2} , \ldots , X_{n}$ are i.i.d., and each have this exponential distribution, then their join distribution is

The joint distribution represents, roughly, the likelihood of any collection of $x_{i}$ values.

The maximum likelihood method turns this around in the sense that, if we have a sample of data, say $X_{1} , X_{2} , \ldots , X_{n}$, but λ is not known, then the maximum-likelihood estimator (MLE) is the value of $\lambda$ that maximizes the likelihood of the sample, as measured by Equation (7). Thus, the MLE of λ for the exponential distribution maximizes the likelihood function

Finding the MLE is often a calculus exercise of maximizing the likelihood function with respect to the unknown parameters. Because the likelihood function is a product of density functions and it is positive, it is sometimes easier to maximize the log of the likelihood function, which in this case is

Taking the derivative with respect to $\lambda$ and setting it equal to $_ 0$ we get

giving the MLE $\begin{array} {r} {\widehat {\lambda} = n / \sum_{i = 1} ^ {n} X_{i} = 1 / \bar {X}} \end{array}$ where $\bar {X}$ is the sample mean of the observations.

For example, assuming that the data in Example 6 come from an exponential distribution, the MLE of λ is obtained from $\bar {X}$ as

Suppose that a random sample of size n $. X_{1} , X_{2} , \ldots , X_{n} .$, has been taken and that the observations are assumed to come from a Weibull distribution. The likelihood function can be shown to be

The maximum-likelihood estimates are those values of $\widehat {\alpha}$ and $\widehat {\beta}$ that maximize $L ( \alpha , \beta )$ or, equivalently, maximize ln $L ( \alpha , \beta )$, denoted by $l ( \alpha , \beta )$. The maximum value of $l ( \alpha , \beta )$ is obtained by taking the partial derivatives $\partial l ( \alpha , \beta ) / \partial \alpha$ and $\partial l ( \alpha , \beta ) / \partial \beta$, setting each to zero, and solving the resulting equations, which, after substitution, become

and

where

The maximum-likelihood estimates,$\widehat {\alpha}$ and ${\widehat {\beta}} ,$, are the solutions of Equations (9) and (10). First,$\widehat {\beta}$is found via the iterative procedure explained below. Then $\widehat {\alpha}$ is found from Equation (10), with $\beta = {\widehat {\beta}}$Equation (9) is nonlinear, so it is necessary to use a numerical-analysis technique to solve it. In Table $^ {3 ,}$, a suggested iterative method for computing $\widehat {\beta}$ is given as

Equation (12) employs Newton’s method in reaching ${\widehat {\beta}} ,$ where ${\widehat {\beta}} _ {j}$ is thejth iteration, beginning with an initial estimate for ${\widehat {\beta}} _ {0}$, given in Table 3, as follows:

If the initial estimate,${\widehat {\beta}} _ {0} .$, is sufficiently close to the solution ${\widehat {\beta}} ,$, then ${\widehat {\beta}} _ {j}$ approaches $\widehat \beta \operatorname{a s} j \longrightarrow \infty$ In Newton’s method,$\widehat {\beta}$ is approached through increments of size $f ( \widehat {\beta} _ {j - 1} ) / f^{\prime} ( \widehat {\beta} _ {j - 1} )$. Equation (11) is used to compute $f ( \widehat {\beta} _ {j - 1} )$ and Equation (14) is used to compute $f^{\prime} ( \widehat {\beta} _ {j - 1} )$, as follows:

Equation (14) can be derived from Equation (11) by differentiating $f ( \beta )$ with respect to $\beta .$. The iterative process continues until $f ( \widehat {\beta} _ {j} ) \doteq 0$, that is, unti $| f ( \widehat {\beta} _ {j} ) | \leq 0 . 001 $ Consider the data given in Example 6. These data concern the failure of electronic components and look to come from an exponential distribution. In Example 10, the parameter $\widehat {\lambda}$ was estimated on the hypothesis that the data were from an exponential distribution. If this hypothesis is rejected, an alternative hypothesis is that the data come from a Weibull distribution. The Weibull distribution is suspected because the data pertain to electronic component failures, which occur suddenly.

Equation (13) is used to compute $\widehat {\beta} _ {0}$. For the data in Example $6 , n = 50 , \bar {X} = 11 . 894 , \bar {X} ^ {2} =$$141 . 467$, and $\begin{array} {r} {\sum_{i = 1} ^ {50} X_{i} ^ {2} = 37 , 575 . 850 ;} \end{array}$; so $S^{2}$ is found by Equation (2) to be

and $S = 24 . 953$. Thus,

To compute ${\widehat {\beta}} _ {1}$ by using Equation (12) requires the calculation of $f ( \widehat {\beta_{0}} ) \mathrm{a n d} f^{\prime} ( \widehat {\beta_{0}} )$ from Equations (11) and (14). The following additional values are needed:$\textstyle \sum_{i = 1} ^ {50} X_{i} ^ {\widehat {\beta_{0}}} = 115 . 125 , \textstyle \sum_{i = 1} ^ {50} \ln X_{i} = 38 . 294$$\begin{array} {r} {\sum_{i = 1} ^ {50} X_{i} ^ {\widehat {\beta_{0}}} \ln X_{i} = 292 . 629} \end{array}$, and $\textstyle \sum_{i = 1} ^ {50} X_{i} ^ {\widehat {\beta_{0}}} ( \ln X_{i} ) ^ {2} = 1057 . 781$. Thus,

and

Then, by Equation (12),

After four iterations,$| f ( \widehat {\beta} _ {3} ) | \le 0 . 001 $, at which point $\widehat {\beta} \doteq \widehat {\beta} _ {4} = 0 . 525$, is the approximate   solution to Equation (9). Table 4 contains the values needed to complete each iteration.

Now,$\widehat {\alpha}$ can be computed from Equation (10) with $\beta = \widehat {\beta} = 0 . 525$, as follows:

Table 4 Iterative Estimation of Parameters of the Weibull Distribution

If ${\widehat {\beta}} _ {0}$ is sufficiently close to ${\widehat {\beta}} ,$ the procedure converges quickly, usually in four to five iterations.  However, if the procedure appears to be diverging, try other initial guesses for ${\widehat {\beta}} _ {0} ;$ for example, one-half the initial estimate or twice the initial estimate.

The difficult task of estimating parameters for the Weibull distribution by hand emphasizes the value of having software support for input modeling.

Example 12: Poisson Distribution

Assume that the arrival data in Table 1 require analysis. An examination of Figure 4 suggests a Poisson distribution assumption with unknown parameter α. From Table 3, the estimator of α is ${\bar {X}} ,$ which was found in Example 8. Thus,$\widehat {\alpha} = 3 . 64$. Recall that the true mean and variance are equal for the Poisson distribution. In Example 8, the sample variance was estimated as $S^{2} = 7 . 63$. However, it should never be expected that the sample mean and the sample variance will be precisely equal, because each is a random variable.

Example 13: Lognormal Distribution

The rates of return on 10 investments in a portfolio are 18.8, 27.9, 21.0, 6.1, 37.4, 5.0, 22.9, 1.0, 3.1, and 8.3 percent. To estimate the parameters of a lognormal model of these data, we first take the natural log of the data and obtain 2.9, 3.3, 3.0, 1.8, 3.6, 1.6, 3.1, 0, 1.1, and 2.1. Then we set $\widehat {\mu} = \bar {X} = 2 . 3$ and $\widehat {\sigma} ^ {2} = S^{2} = 1 . 3$

Example 14: Normal Distribution

The parameters of the normal distribution,$\mu$ and $\sigma^{2}$, are estimated by X¯ and $S^{2}$, as shown in Table 3. The $q {-} q$ plot in Example 7 leads to a distributional assumption that the installation times are normal. From Equations (1) and (2), the data in Example 7 yield $\widehat {\mu} = \bar {X} = 99 . 9865$ and $\widehat {\sigma} ^ {2} = S^{2} = ( 0 . 2832 ) ^ {2}$ seconds^2.

Example 15: Gamma Distribution

The estimator $\widehat {\beta}$ for the gamma distribution is chosen by referring to Table A.9. Table A.9 requires the computation of the quantity $1 / M$, where

Also, it can be seen in Table 3 that $\widehat {\theta}$ is given by

Lead time is often gamma-distributed. Suppose that the lead times (in days) associated with 20 orders have been accurately measured as follows:

To estimate $\widehat {\beta}$ and ${\widehat {\theta}} ,$ it is first necessary to compute M from Equation (15). Here,$\bar {X}$ is found, from Equation (1), to be

Then

Next

Then

and

By interpolation in Table A.9 in the Appendix,${\widehat {\beta}} = 3 . 728 .$. Finally, Equation (16) results in

Example 16: Beta Distribution

The percentage of customers each month who bring in store coupons must be between 0 and 100 percent. Observations at a store for eight months gave the values 25%, 74%, 20%, 32%, 81%, 47%, 31%, and 8%. To fit a beta distribution to these data, we first need to rescale it to the interval [0, 1] by dividing all the values by 100, which yields 0.25, 0.74, 0.20, 0.32, 0.81, 0.47, 0.31, 0.08.

The maximum-likelihood estimators of the parameters $\beta_{1} , \beta_{2}$ solve the system of equations shown in Table 9.3. Such equations can be solved by modern symbolic/numerical calculation programs, such as Maple; a Maple procedure for the beta parameters is shown in Figure 7. In this case, the solutions are $\widehat {\beta} _ {1} = 1 . 47$ and $\widehat {\beta} _ {2} = 2 . 16$```txt betaMLE := proc(X, n) local G1, G2, beta1, beta2, eqns, solns; G1 := product(X[i], i=1..n)^(1/n); G2 := product(1-X[i], i=1..n)^(1/n); eqns := {Psi(beta1) - Psi(beta1 + beta2) = ln(G1), Psi(beta2) - Psi(beta1 + beta2) = ln(G2)}; solns := fsolve(eqns, {beta1=0..infinity, beta2=0..infinity}); RETURN(solns); end;

Figure 7 Maple procedure to compute the maximum likelihood estimates for the beta distri bution parameters.

## 4 Goodness-of-Fit Tests

In this section, the Kolmogorov–Smirnov test and the chi-square test are applied to hypotheses abou distributional forms of input data.


Goodness-of-fit tests provide helpful guidance for evaluating the suitability of a potential inpu model; however, there is no single correct distribution in a real application, so you should not be a slave to the verdict of such a test. 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, and rejecting an input model as only one piece of evidence against the choice.

### 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 and for both discrete and continuous distributiona 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 

$$
\chi_{0} ^ {2} = \sum_{i = 1} ^ {k} \frac {(O_{i} - E_{i}) ^ {2}}{E_{i}}\tag{17}
$$

 where \(O_{i}\) is the observed frequency in the ith class interval and \(E_{i}\) is the expected frequency in that class interval. The expected frequency for each class interval is computed as \(E_{i} = n p_{i}\), where \(p_{i}\) is the theoretical, hypothesized probability associated with the ith class interval.

It can be shown that \(\chi_{0} ^ {2}\) approximately follows the chi-square distribution with \(k - s - 1\) degrees of freedom, where s represents the number of parameters of the hypothesized distribution estimated by the sample statistics. The hypotheses are the following:\(H_{0} {\mathrm{:}}\) The random variable, X, conforms to the distributional assumption with the parameter(s) given by the parameter estimate(s).\(H_{1} {\mathrm{:}}\): The random variable X does not conform

The critical value \(\chi_{\alpha , k - s - 1} ^ {2}\) is found in Table A.6. The null hypothesis \(H_{0}\) is rejected if \(\chi_{0} ^ {2} > \chi_{\alpha , k - s - 1} ^ {2} .\) When applying the test, if expected frequencies are too small,\(\chi_{0} ^ {2}\) will reflect not only the departure of the observed from the expected frequency, but also the smallness of the expected frequency as well. Although there is no general agreement regarding the minimum size of \(E_{i} ,\) values of \(3 , 4 .\), and 5 have been widely used. If an \(E_{i}\) value is too small, it can be combined with expected frequencies in adjacent class intervals. The corresponding \(O_{i}\) values should also be combined, and k should be reduced by 1 for each cell that is combined.

If the distribution being tested is discrete, each value of the random variable 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,

$$
p_{i} = p (x_{i}) = P (X = x_{i})
$$

 Otherwise,\(p_{i}\) is found by summing the probabilities of appropriate adjacent cells.

Table 5 Recommendations for Number of Class Intervals for Continuous Data

| Sample Size n | Number of Class Intervals k |
| --- | --- |
| 20 | Do not use the chi-square test |
| 50 | 5 to 10 |
| 100 | 10 to 20 |
| >100 |\(\sqrt{n}\) to \(n/5\)|

If the distribution being tested is continuous, the class intervals are given by \([ a_{i - 1} , a_{i} )\), where \(a_{i - 1}\) and \(a_{i}\) are the endpoints of the ith class interval. For the continuous case with assumed pdff(x), or assumed cdf \(F ( x ) , p_{i}\) can be computed as 

$$
p_{i} = \int_ {a_{i - 1}} ^ {a_{i}} f (x) d x = F (a_{i}) - F (a_{i - 1})
$$

 For the discrete case, the number of class intervals is determined by the number of cells resulting after combining adjacent cells as necessary. However, 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 5 are made to aid in determining the number of class intervals for continuous data.

### Example 17: Chi-Square Test Applied to Poisson Assumption

In Example 12, the vehicle-arrival data presented in Example 5 were analyzed. The histogram of the data, shown in Figure \(^ {4 ,}\) appeared to follow a Poisson distribution; hence the parameter \(\widehat {\alpha} = 3 .\).64 was found. Thus, the following hypotheses are formed:\(H_{0} {\mathrm{:}}\) The random variable is Poisson distributed.\(H_{1}\): The random variable is not Poisson distributed.

The pmf for the Poisson distribution:

$$
p (x) = \left\{\begin{array}{l l} \frac {e^{- \alpha} \alpha^{x}}{x !}, & x = 0, 1, \ldots \\ 0, & \text{otherwise} \end{array} \right.\tag{18}
$$

 For \(\widehat {\alpha} = 3 . 64\), the probabilities associated with various values of x are obtained from Equation (18):

$$
\begin{array}{l l} \mathrm{p(0)=0.026} & \mathrm{p(6)=0.085} \\ \mathrm{p(1)=0.096} & \mathrm{p(7)=0.044} \\ \mathrm{p(2)=0.174} & \mathrm{p(8)=0.020} \\ \mathrm{p(3)=0.211} & \mathrm{p(9)=0.008} \\ \mathrm{p(4)=0.192} & \mathrm{p(10)=0.003} \\ \mathrm{p(5)=0.140} & \mathrm{p(\geq 11) = 0.001} \end{array}
$$

 Table 6 Chi-Square Goodness-of-Fit Test for Example 17

|\(x_i\)| Observed Frequency \(O_i\)| Expected Frequency \(E_i\)|\((O_i - E_i)^2\)\(E_i\)|  |
| --- | --- | --- | --- | --- |
| 0 | 12 | 22 | 2.6 | 7.87 |
| 1 | 10 | 9.6 |  |  |
| 2 | 19 | 17.4 | 0.15 |  |
| 3 | 17 | 21.1 | 0.80 |  |
| 4 | 10 | 19.2 | 4.41 |  |
| 5 | 8 | 14.0 | 2.57 |  |
| 6 | 7 | 8.5 | 0.26 |  |
| 7 | 5 | 17 | 4.4 | 11.62 |
| 8 | 5 | 2.0 |  |  |
| 9 | 3 | 0.8 |  |  |
| 10 | 3 | 0.3 |  |  |
| ≥11 | 1 | 0.1 |  |  |
|  | 100 | 100.0 | 27.68 |  |

From this information, Table 6 is constructed. The value of \(E_{1}\) is given by \(n p_{0} = 100 ( 0 . 026 ) =\)2.6. In a similar manner, the remaining \(E_{i}\) values are computed. Since \(E_{1} = 2 . 6 < 5 , E_{1}\) and \(E_{2}\) are combined. In that case,\(O_{1}\) and \(O_{2}\) 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 \(\chi_{0} ^ {2}\) is 27.68. The degrees of freedom for the tabulated value of \(\chi^{2}\) is \(k - s - 1 =\)\(7 - 1 - 1 = 5\). Here,\(s = 1\), since one parameter,\({\widehat {\alpha}} ,\), was estimated from the data. At the 0.05 level of significance, the critical value \(\chi_{0 . 05 , 5} ^ {2}\)is 11.1. Thus,\(H_{0}\) would be rejected at level of significance \(0 . 05\). The analyst, therefore, might want to search for a better-fitting model or use the empirical distribution of the data.

### 4.2 Chi-Square Test with Equal Probabilities

If a continuous distributional assumption is being tested, class intervals that are equal in probability rather 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, 2007; 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 or lost.

Unfortunately, there is no easy method for figuring out the probability associated with each interval that would maximize 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_{i} = 1 / k\) We recommend 

$$
E_{i} = n p_{i} \geq 5
$$


so, substituting for \(p_{i}\) yields 

$$
\frac {n}{k} \geq 5
$$

 and solving for k yield 

$$
k \leq \frac {n}{5}\tag{19}
$$

 Equation (19) was used in coming up with the recommendations for maximum number of class intervals in Table 5.

If the assumed distribution is normal, exponential, or Weibull, the method described in this section is straightforward. Example 18 indicates how the procedure is accomplished for the exponential distribution. If the assumed distribution is gamma (but not Erlang) or certain other distributions, then the computation of endpoints for class intervals is complex and could require numerica integration of the density function. Statistical-analysis software is very helpful in such cases.

### Example 18: Chi-Square Test for Exponential Distribution

In Example 10, the failure data presented in Example 6 were analyzed. The histogram of the data, shown in Figure 5, appeared to follow an exponential distribution, so the parameter \(\widehat {\lambda} = 1 / \bar {X} = 0 . 084\) was computed. Thus, the following hypotheses are formed:\(H_{0} {\mathrm{:}}\) The random variable is exponentially distributed.\(H_{1} ;\): 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 found. Equation (19) indicates that the number of intervals should be less than or equal to \(n / 5 .\). Here,\(n = 50\), and so \(k \leq 10\). In Table 5, it is recommended that 5 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, as follows:

$$
F (a_{i}) = 1 - e ^ {- \lambda a_{i}}\tag{20}
$$

 where \(a_{i}\) represents the endpoint of the ith interval,\(i = 1 , 2 , \ldots , k\). Since \(F ( a_{i} )\) is the cumulative area from zero to \(a_{i} , F ( a_{i} ) = i p\), so Equation (20) can be written as 

$$
i p = 1 - e ^ {- \lambda a_{i}}
$$

 or 

$$
e ^ {- \lambda a_{i}} = 1 - i p
$$

 Taking the logarithm of both sides and solving for \(a_{i}\) gives a general result for the endpoints of k equiprobable intervals for the exponential distribution:

$$
a_{i} = - \frac {1}{\lambda} \ln (1 - i p), \quad i = 0, 1, \dots , k\tag{21}
$$

 Table 7 Chi-Square Goodness-of-Fit Test for Example 18

| Class Interval | Observed Frequency \(O_i\)| Expected Frequency \(E_i\)|\((O_i - E_i)^2\)\(E_i\)|
| --- | --- | --- | --- |
| [0, 1.590) | 19 | 6.25 | 26.01 |
| [1.590, 3.425) | 10 | 6.25 | 2.25 |
| [3.425, 5.595) | 3 | 6.25 | 0.81 |
| [5.595, 8.252) | 6 | 6.25 | 0.01 |
| [8.252, 11.677) | 1 | 6.25 | 4.41 |
| [11.677, 16.503) | 1 | 6.25 | 4.41 |
| [16.503, 24.755) | 4 | 6.25 | 0.81 |
| [24.755,∞) | 6 | 6.25 | 0.01 |
|  | 50 | 50 | 39.6 |

Regardless of the value of \(\lambda ,\) Equation (21) will always result in \(a_{0} = 0\) and \(a_{k} = \infty\). With \(\widehat {\lambda} = 0 . 084\) and \(k = 8 , a_{1}\) is computed from Equation (21) as 

$$
a_{1} = - \frac {1}{0 . 084} \ln (1 - 0. 125) = 1. 590
$$

 Continued application of Equation (21) for \(i = 2 , 3 , \dots , 7\) results in \(a_{2} , \ldots , a_{7}\) as 3.425, 5.595, 8.252, 11.677, 16.503, and 24.755. Since \(k = 8 , a_{8} = \infty\). 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 observations, the expectations, and the contributions to the calculated value of \(\chi_{0} ^ {2}\) are shown in Table 7. The calculated value of \(\chi_{0} ^ {2}\) is 39.6. The degrees of freedom are given by \(k - s - 1 = 8 - 1 - 1 = 6 . \mathrm{A t} \alpha = 0 . 05\), the tabulated value of \(\chi_{0 . 05 , 6} ^ {2}\) is 12.6. Since \(\chi_{0} ^ {2} > \chi_{0 . 05 , 6} ^ {2} ,\) the null hypothesis is rejected. (The value of \(\chi_{0 . 01 , 6} ^ {2}\) is 16.8, so the null hypothesis would also be rejected at level of significance \(\alpha = 0 . 01 . )\)
### 4.3 Kolmogorov--Smirnov Goodness-of-Fit Test

The chi-square goodness-of-fit test can accommodate the estimation of parameters from the data with a resultant decrease in the degrees of freedom (one for each parameter estimated). The chisquare test requires that the data be placed in class intervals; in the case of a continuous distributiona assumption, this grouping is arbitrary. Changing the number of classes and the interval width affects the value of the calculated and tabulated chi-square. A hypothesis could be accepted when the data are grouped one way, but rejected when they are grouped another way. Also, the distribution of the chi-square test statistic is known only approximately, and the power of the test is sometimes rather low. As a result of these considerations, goodness-of-fit tests other than the chi-square, are desired. The Kolmogorov–Smirnov test formalizes the idea behind examining a q-q plot.


The Kolmogorov–Smirnov test can be presented to test for the uniformity of numbers. Both of these uses fall into the category of testing for goodness of fit. Any continuous distributiona assumption can be tested for goodness of fit.

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 critica 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 (α) errors than those specified. The exact value of α can be worked out in some instances, as is 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, so we may use Table A.8.)

Example 19: Kolmogorov–Smirnov Test for Exponential Distribution
Suppose that 50 interarrival times (in minutes) are collected over a 100-minute interval (arranged in order of occurrence):

| 0.44 | 0.53 | 2.04 | 2.74 | 2.00 | 0.30 | 2.54 | 0.52 | 2.02 | 1.89 | 1.53 | 0.21 |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| 2.80 | 0.04 | 1.35 | 8.32 | 2.34 | 1.95 | 0.10 | 1.42 | 0.46 | 0.07 | 1.09 | 0.76 |
| 5.55 | 3.93 | 1.07 | 2.26 | 2.88 | 0.67 | 1.12 | 0.26 | 4.57 | 5.37 | 0.12 | 3.19 |
| 1.63 | 1.46 | 1.08 | 2.06 | 0.85 | 0.83 | 2.44 | 1.02 | 2.24 | 2.11 | 3.15 | 2.90 |
| 6.58 | 0.64 |  |  |  |  |  |  |  |  |  |  |

The null hypothesis and its alternate are formed as follows:\(H_{0} {\mathrm{:}}\) The interarrival times are exponentially distributed.\(H_{1}\): The interarrival times are not exponentially distributed.

The data were collected over the interval from 0 to \(T = 100\) minutes. It can be shown that, if the underlying distribution of interarrival times \(\{T_{1} , T_{2} , \dots \}\) is exponential, the arrival times are uniformly distributed on the interval [0, T]. The arrival times \(T_{1} , T_{1} + T_{2} , T_{1} + T_{2} + T_{3} , . . . ,\)\(T_{1} + \cdots + T_{50}\) are obtained by adding interarrival times. The arrival times are then normalized to a [0, 1] interval so that the Kolmogorov–Smirnov test can be applied. On a [0, 1] interval, the points will be \([ T_{1} / T , ( T_{1} + T_{2} ) / T , \dots , ( T_{1} + \cdot \cdot \cdot + T_{50} ) / 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 | 0.1502 |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| 0.1655 | 0.1676 | 0.1956 | 0.1960 | 0.2095 | 0.2927 | 0.3161 | 0.3356 | 0.3366 | 0.3508 |
| 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 | 0.8732 | 0.9022 | 0.9680 | 0.9744 |

Following the procedure produces a \(D^{+}\) of 0.1054 and a \(D^{-}\) of 0.0080. Therefore, the Kolmogorov–Smirnov statistic is \(D \ = \ \operatorname* {m a x} ( 0 . 1054 , 0 . 0080 ) \ = \ 0 . 1054\). The critical value


of D obtained from Table A.8 for a level of significance of \(\alpha \ : = \ : 0 . 05\) and \(n = 50\) is \(D_{0 . 05} \ =\)\(1 . 36 / \sqrt {n} = 0 . 1923\); however,\(D = 0 . 1054\), so the hypothesis that the interarrival times are exponentially distributed cannot be rejected.

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 computation of the test statistic is the same, bu different tables of critical values are used. Different tables of critical values are required for different distributional 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 might wish to study Lilliefors’ original work; 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. Connover [1998] gives examples o Kolmogorov–Smirnov tests for the normal and exponential distributions. He also refers to severa other Kolmogorov–Smirnov-type tests that might 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 difference) 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 parameters have been estimated. Fortunately, this test and the Kolmogorov–Smirnov test have been implemented in a number of software packages that support simulation-input modeling.

### 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 \(H_{0} {\mathrm{:}}\): The random variable conforms to the distributiona assumption. The traditional significance levels are 0.1, 0.05, and 0.01. Prior to the availability o 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, the analyst can employ a different level of significance, say, 0.07.

However, rather than require a prespecified significance level, many software packages compute a \(p \mathrm{.}\)-value for the test statistic. The \(p \mathrm{-}\)-value is the significance level at which one would just reject \(H_{0}\) 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 17, in which a chi-square test was used to check the Poisson assumption for the vehicle-arrival data. The value of the test statistic was \(\chi_{0} ^ {2} = 27 . 68\), with 5 degrees of freedom. The \(p \mathrm{-}\) 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 have to reject it at even lower levels.)

The p-value can be viewed as a measure of fit, with larger values being better. This suggests tha we could fit every distribution at our disposal, compute a test statistic for each fit, and then choose the distribution that yields the largest p-value. We know of no input modeling software that implements this specific algorithm, but many such packages do include a “best fit” option, in which the software recommends an input model to the user after evaluating all feasible models. The software might also take into account other factors—such as whether the data are discrete or continuous, bounded or unbounded—but, in the end, some summary measure of fit, like 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 might know nothing about the physical basis of the data, whereas that information can suggest distribution families that are appropriate. (See the list in Section 2.2.) Remember that the goal of input modeling is often to fill in gaps or smooth the data, rather than to 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 and that 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 and a better summary measure of fit. But again, close conformance to the data does 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.

## 5 Fitting a Nonstationary Poisson Process

Fitting a nonstationary Poisson process (NSPP) to arrival data is a difficult problem, in general, because we seldom have knowledge about the appropriate form of the arrival rate function λ(t). One approach is to choose a very flexible model with lots of parameters and fit it with a method such as maximum likelihood; see Johnson, Lee, and Wilson [1994] for an example of this approach. A second method, and the one we consider here, is to approximate the arrival rate as being constant over some basic interval of time, such as an hour, or a day, or a month, but varying from time interval to time interval. The problem then becomes choosing the basic time interval and estimating the arriva rate within each interval.

Suppose we need to model arrivals over a time period, say [0, T]. The approach that we describe is most appropriate when it is possible to observe the time period [0, T] repeatedly and count arrivals. For instance, if the problem involves modeling the arrival of e-mail throughout the business day (8 A.M. to 6 P.M.), and we believe that the arrival rate is approximately constant over half-hour intervals, then we need to be able to count arrivals during half-hour intervals for several days. If it is possible to record actual arrival times, rather than counts, then actual arrival times are clearly better since they can later be grouped into any interval lengths we desire. However, we will assume from here on that only counts are available.

Table 8 Monday E-mail Arrival Data for NSPP Example

| Time Period | Number of Arrivals | Estimated Arrival Rate (arrivals/hour) |  |  |
| --- | --- | --- | --- | --- |
| Day 1 | Day 2 | Day 3 |  |  |
| 8:00–8:30 | 12 | 14 | 10 | 24 |
| 8:30–9:00 | 23 | 26 | 32 | 54 |
| 9:00–9:30 | 27 | 19 | 32 | 52 |
| 9:30–10:00 | 20 | 13 | 12 | 30 |

Divide the time period [0, T] into k equal intervals of length \(\Delta t = T / k\). For instance, if we are considering a 10-hour business day from 8 A.M. to 6 P.M. and if we allow the rate to change every half hour, then \(T = 10 , k = 20\), and \(\Delta t = 1 / 2\). Over n periods of observation (e.g., n days), let \(C_{i j}\) be the number of arrivals that occurred during the ith time interval on the jth period of observation. In our example,\(C_{23}\) would be the number of arrivals from 8:30 A.M. to 9 A.M. (second half-hour period) on the third day of observation.

The estimated arrival rate during the ith time period,\(( i - 1 ) \Delta t < t \leq i \Delta t .\), is then just the average number of arrivals scaled by the length of the time interval:

$$
\widehat {\lambda} (t) = \frac {1}{n \Delta t} \sum_{j = 1} ^ {n} C_{i j}\tag{22}
$$

 After the arrival rates for each time interval have been estimated, adjacent intervals whose rates appear to be the same can be combined.

For instance, consider the e-mail arrival counts during the first two hours of the business day on three Mondays, shown in Table 8. The estimated arrival rate for 8:30–9:00 is 

$$
\frac {1}{3 (1 / 2)} (23 + 26 + 32) = 54 \text{arrivals / hour}
$$

 After seeing these results, we might consider combining the interval 8:30–9:00 with the interva 9:00–9:30, because the rates are so similar. Note also that the goodness-of-fit tests described in the previous section can be applied to the data from each time interval individually, to check the Poisson approximation.

## 6 Selecting Input Models without Data

Unfortunately, it is often necessary in practice to develop a simulation model—perhaps for demonstration purposes 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.


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 provided by the manufacturer (for example, the mean time to failure of a disk drive is 10,000 hours; a laser printer can produce 8 pages/minute; the cutting speed of a tool is 1 cm/second; etc.). Company rules might 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 might also be able to say whether the process is nearly constant or highly variable, and they might 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 could 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 descriptions of the distributions in Section 2.2 can be used tojustify 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 jus as likely as the central values in real processes. If, in addition to upper and lower bounds, a mostlikely 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. If a beta distribution is used, then be sure to plot the density function of the selected distribution; the beta can take unusual shapes.

A useful refinement is obtained when a minimum, a maximum, and one or more “breakpoints” can be given. A breakpoint is an intermediate value together with a probability of being less than or equal to that value. The following example illustrates how breakpoints are used.

### Example 20

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). Given 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 summarizes this information. Notice that the chances of exceeding certain sales goals have been translated into the cumulative probability of being less than or equal to those goals.

Table 9 Summary of Sales Information

| i | Interval (sales) | Cumulative Frequency \(c_i\)|
| --- | --- | --- |
| 1 |\(1000 \le x \le 2000\)| 0.10 |
| 2 |\(2000 < x \le 3500\)| 0.75 |
| 3 |\(3500 < x \le 4500\)| 0.99 |
| 4 |\(4500 < x \le 5000\)| 1.00 |

When input models have been selected without data, it is especially important to test the sensitivity of simulation results to the distribution chosen. Check sensitivity not only to the center of the distribution but also to the variability or limits. Extreme sensitivity of output results to the input model provides a convincing argument against making critical decisions based on the results and in favor of undertaking data collection

## 7 Multivariate and Time-Series Input Models

In Sections 1–4, the input random variables were considered to be independent of any other variables within the context of the problem. When this is not the case it is critical that input models be used that account for dependence; otherwise, a highly inaccurate simulation may result.

### Example 21

A supply-chain simulation includes the lead time and annual demand for industrial robots.An increase in demand results in an increase in lead time: The final assembly of the robots must be made according to the specifications ofthe purchaser. Therefore, rather than treat lead time and demand as independen random variables, a multivariate input model should be developed.

### Example 22

A simulation of the web-based trading site of a stock broker includes the time between arrivals of orders to buy and sell. Investors tend to react to what other investors are doing, so these buy and sel orders arrive in bursts. Therefore, rather than treat the time between arrivals as independent random variables, a time-series model should be developed.

We distinguish multivariate input models of a fixed, finite number of random variables (such as the two random variables lead time and annual demand in Example 21) from time-series input models of a (conceptually infinite) sequence of related random variables (such as the successive times between orders in Example 22). We will describe input models appropriate for these examples after reviewing two measures of dependence, the covariance and the correlation.

### 7.1 Covariance and Correlation

Let \(X_{1}\) and \(X_{2}\) be two random variables, and let \(\mu_{i} = \operatorname{E} ( X_{i} )\) and \(\sigma_{i} ^ {2} = \mathrm{V} ( X_{i} )\) be the mean and variance of \(X_{i} ,\) respectively. The covariance and correlation are measures of the linear dependence between \(X_{1}\) and \(X_{2} .\). In other words, the covariance and correlation indicate how well the relationship between \(X_{1}\) and \(X_{2}\) is described by the model 

$$
(X_{1} - \mu_{1}) = \beta (X_{2} - \mu_{2}) + \varepsilon
$$

 where ε is a random variable with mean 0 that is independent of \(X_{2}\). If, in fact,\(( X_{1} - \mu_{1} ) = \beta ( X_{2} - \mu_{2} )\) then this model is perfect. On the other hand,\(\operatorname{i f} X_{1}\) and \(X_{2}\) are statistically independent, then \(\beta = 0\) and the model is of no value. In general, a positive value of \(\beta\) indicates that \(X_{1}\) and \(X_{2}\) tend to be above or below their means together; a negative value o \(\boldsymbol {\beta}\) indicates that they tend to be on opposite sides of their means.

The covariance between \(X_{1}\) and \(X_{2}\) is defined to be 

$$
\operatorname{cov} \left(X_{1}, X_{2}\right) = \operatorname{E} \left[ \left(X_{1} - \mu_{1}\right) \left(X_{2} - \mu_{2}\right) \right] = \operatorname{E} \left(X_{1} X_{2}\right) - \mu_{1} \mu_{2}\tag{23}
$$

 The value cov \(( X_{1} , X_{2} ) = 0\) implies \(\beta = 0\) in our model of dependence, and \(\operatorname{c o v} ( X_{1} , X_{2} ) < 0 ( > 0 )\) implies \(\beta < 0 ( > 0 )\) The covariance can take any value between and . The correlation standardizes the covariance to be between 1 and 1:

$$
\rho = \operatorname{corr} (X_{1}, X_{2}) = \frac {\operatorname{cov} (X_{1} , X_{2})}{\sigma_{1} \sigma_{2}}\tag{24}
$$

 Again, the value corr \(( X_{1} , X_{2} ) = 0\) implies \(\beta = 0\) in our model, and \(\operatorname{c o r r} ( X_{1} , X_{2} ) < 0 ( > 0 )\) implies \(\beta < 0 ( > 0 )\). The closer \(\rho\) is to 1 or 1, the stronger the linear relationship is between \(X_{1}\) and \(X_{2}\).

Now, suppose that we have a sequence of random variables \(X_{1} , X_{2} , X_{3} , . . .\). that are identically distributed (implying that they all have the same mean and variance), but could be dependent. We refer to such a sequence as a time series and to cov \(( X_{t} , X_{t + h} )\) and \(\operatorname{c o r r} ( X_{t} , X_{t + h} )\) as the lag-h autocovariance and lag-h autocorrelation, respectively. If the value of the autocovariance depends only on h and not on \(t ,\) then we say that the time series is covariance stationary. For a covariance-stationary time series, we use the shorthand notation 

$$
\rho_{h} = \mathrm{corr} (X_{t}, X_{t + h})
$$

 for the lag-h autocorrelation. Notice that autocorrelation measures the dependence between random variables that are separated by \(h - 1\) others in the time series.

### 7.2 Multivariate Input Models

If \(X_{1}\) and \(X_{2}\) each are normally distributed, then dependence between them can be modeled by the bivariate normal distribution with parameters \(\mu_{1} , \mu_{2} , \sigma_{1} ^ {2} , \sigma_{2} ^ {2}\), and \(\rho = \operatorname{c o r r} ( X_{1} , X_{2} )\)). Estimation of \(\mu_{1} , \mu_{2} , \sigma_{1} ^ {2}\), and \(\sigma_{2} ^ {2}\) was described in Section 3.2. To estimate \(\rho ,\), suppose that we have n independent and identically distributed pairs \(( X_{11} , X_{21} ) , ( X_{12} , X_{22} ) , \dots , ( X_{1 n} , X_{2 n} )\). Then the sample covariance is 

$$
\begin{array}{l} \widehat {\operatorname{cov}} (X_{1}, X_{2}) = \frac {1}{n - 1} \sum_{j = 1} ^ {n} \left(X_{1 j} - \bar {X} _ {1}\right) (X_{2 j} - \bar {X} _ {2}) \\ = \frac {1}{n - 1} \left(\sum_{j = 1} ^ {n} X_{1 j} X_{2 j} - n \bar {X} _ {1} \bar {X} _ {2}\right) \end{array}\tag{25}
$$

 where \({\bar {X}} _ {1}\) and \({\bar {X}} _ {2}\) are the sample means. The correlation is estimated by 

$$
\widehat {\rho} = \frac {\widehat {\mathrm{cov}} (X_{1} , X_{2})}{\widehat {\sigma} _ {1} \widehat {\sigma} _ {2}}\tag{26}
$$

 where \(\widehat {\sigma} _ {1}\) and \(\widehat {\sigma} _ {2}\) are the sample variances.

### Example 23: Example 21 (cont.)

Let \(X_{1}\) represent the average lead time (in months) to deliver, and \(X_{2}\) the annual demand, for industria robots. The following data are available on demand and lead time for the last ten years:

| Lead Time | Demand |
| --- | --- |
| 6.5 | 103 |
| 4.3 | 83 |
| 6.9 | 116 |
| 6.0 | 97 |
| 6.9 | 112 |
| 6.9 | 104 |
| 5.8 | 106 |
| 7.3 | 109 |
| 4.5 | 92 |
| 6.3 | 96 |

Standard calculations give \({\bar {X}} _ {1} = 6 . 14 , {\widehat {\sigma}} _ {1} = 1 . 02 , {\bar {X}} _ {2} = 101 . 80 .\), and \(\widehat {\sigma} _ {2} = 9 . 93\) as estimates of \(\mu_{1} , \sigma_{1} , \mu_{2}\), and \(\sigma_{2} .\), respectively. To estimate the correlation, we need 

$$
\sum_{j = 1} ^ {10} X_{1 j} X_{2 j} = 6328. 5
$$

 Therefore,\(\widehat {\mathrm{c o v}} = [ 6328 . 5 - ( 10 ) ( 6 . 14 ) ( 101 . 80 ) ] / ( 10 - 1 ) = 8 . 66\), and 

$$
\widehat {\rho} = \frac {8 . 66}{(1 . 02) (9 . 93)} = 0. 86
$$


Clearly, lead time and demand are strongly dependent. Before we accept this model, however, lead time and demand should be checked individually to see whether they are represented well by normal distributions. In particular, demand is a discrete-valued quantity, so the continuous normal distribution is at best an approximation.

The following simple algorithm can be used to generate bivariate normal random variables:

Step 1. Generate \(Z_{1}\) and \(Z_{2} .\), independent standard normal random variables.

Step 2. Set \(X_{1} = \mu_{1} + \sigma_{1} Z_{1}\) Step 3. Set \(X_{2} = \mu_{2} + \sigma_{2} \left( \rho Z_{1} + \sqrt {1 - \rho^{2}} Z_{2} \right)\) Obviously, the bivariate normal distribution will not be appropriate for all multivariate-input modeling problems. It can be generalized to the k-variate normal distribution to model the dependence among more than two random variables, but, in many instances, a normal distribution is no appropriate in any form. We provide one method for handling nonnormal distributions in Section 7.4. Good references for other models are Johnson [1987] and Nelson and Yamnitsky [1998].

### 7.3 Time-Series Input Models

If \(X_{1} , X_{2} , X_{3} , \dots\). is a sequence of identically distributed, but dependent and covariance-stationary random variables, then there are a number of time-series models that can be used to represent the process. We will describe two models that have the characteristic that the autocorrelations take the form 

$$
\rho_{h} = \mathrm{corr} (X_{t}, X_{t + h}) = \rho^{h}
$$

 for \(h = 1 , 2 , \ldots\). Notice that the \(l a g .\)-h autocorrelation decreases geometrically as the lag increases, so that observations far apart in time are nearly independent. For one model to be shown below, each \(X_{t}\) is normally distributed; for the other model, each \(X_{t}\) is exponentially distributed. More general time-series input models are described in Section 7.4 and in Nelson and Yamnitsky [1998].

### AR(1) MODEL

Consider the time-series model 

$$
X_{t} = \mu + \phi (X_{t - 1} - \mu) + \varepsilon_{t}\tag{27}
$$

 for \(t = 2 , 3 , \ldots\), where \(\varepsilon_{2} , \varepsilon_{3} , \ldots .\). are independent and identically (normally) distributed with mean 0 and variance \(\sigma_{\varepsilon} ^ {2}\), and \(- 1 < \phi < 1\). If the initial value \(X_{1}\) is chosen appropriately (see below), then \(X_{1} , X_{2} , . . .\). are all normally distributed with mean \(\mu ,\), variance \(\sigma_{\varepsilon} ^ {2} / ( 1 - \phi^{2} )\), and 

$$
\rho_{h} = \phi^{h}
$$

 for \(h = 1 , 2 , \ldots\). This time-series model is called the autoregressive order-1 model, or AR(1) for short.


Estimation of the parameter \(\phi\) can be obtained from the fact that 

$$
\phi = \rho^{1} = \mathrm{corr} (X_{t}, X_{t + 1})
$$

 the lag-1 autocorrelation. Therefore, to estimate \(\phi ,\), we first estimate the lag-1 autocovariance by 

$$
\widehat {\mathrm{cov}} (X_{t}, X_{t + 1}) = \frac {1}{n - 1} \sum_{t = 1} ^ {n - 1} (X_{t} - \bar {X}) (X_{t + 1} - \bar {X})
$$


$$
\doteq \frac {1}{n - 1} \left(\sum_{t = 1} ^ {n - 1} X_{t} X_{t + 1} - (n - 1) \bar {X} ^ {2}\right)\tag{28}
$$

 and the variance \(\sigma^{2} = \operatorname{v a r} ( X )\) by the usual estimator \({\widehat {\sigma}} ^ {2}\). Then 

$$
\widehat {\phi} = \frac {\widehat {\mathrm{cov}} (X_{t} , X_{t + 1})}{\widehat {\sigma} ^ {2}}
$$

 Finally, estimate \(\mu\) and \(\sigma_{\varepsilon} ^ {2}\) by \({\widehat {\mu}} = {\bar {X}}\) and \(\widehat {\sigma} _ {\varepsilon} ^ {2} = \widehat {\sigma} ^ {2} ( 1 - \widehat {\phi} ^ {2} )\)), respectively.

  The following algorithm generates a stationary \(\mathbf{A R} ( 1 )\)) time series, given values of the parameters \(\phi , \mu\), and \(\sigma_{\varepsilon} ^ {2} \mathrm{:}\):

Step 1. Generate \(X_{1}\) from the normal distribution with mean \(\mu\) and variance \(\sigma_{\varepsilon} ^ {2} / ( 1 - \phi^{2} )\). Set \(t = 2\) Step 2. Generate \(\varepsilon_{t}\) from the normal distribution with mean \(_ 0\) and variance \(\sigma_{\varepsilon} ^ {2}\) Step 3. Set \(X_{t} = \mu + \phi ( X_{t - 1} - \mu ) + \varepsilon_{t}\) Step 4. Set \(t = t + 1\) and go to Step 2.

### EAR(1) MODEL

Consider the time-series model 

$$
X_{t} = \left\{\begin{array}{l l} \phi X_{t - 1}, & \text{with   probability} \phi \\ \phi X_{t - 1} + \varepsilon_{t}, & \text{with   probability} 1 - \phi \end{array} \right.\tag{29}
$$

 for \(t = 2 , 3 , \ldots\), where \(\varepsilon_{2} , \varepsilon_{3} , . . .\). are independent and identically (exponentially) distributed with mean \(1 / \lambda\) and \(0 \leq \phi < 1\). If the initial value \(X_{1}\) is chosen appropriately (see below), then \(X_{1} , X_{2} , \dots\) are all exponentially distributed with mean \(1 / \lambda\) and 

$$
\rho_{h} = \phi^{h}
$$

 for \(h = 1 , 2 , \ldots\). This time-series model is called the exponential autoregressive order-1 model, or \(\mathrm{E A R} ( 1 )\) for short. Only autocorrelations greater than 0 can be represented by this model. Estimation of the parameters proceeds as for the \(\mathbf{A R} ( 1 )\) by setting \({\widehat {\phi}} = {\widehat {\rho}} .\), the estimated lag-1 autocorrelation, and setting \(\widehat {\lambda} = 1 / \bar {X}\)

The following algorithm generates a stationary EAR(1) time series, given values of the parameters φ and λ:

Step 1. Generate \(X_{1}\) from the exponential distribution with mean \(1 / \lambda\). Set \(t = 2\) Step 2. Generate U from the uniform distribution on [0, 1]. If \(U \leq \phi .\), then set 

$$
X_{t} = \phi X_{t - 1}
$$

 Otherwise, generate \(\varepsilon_{t}\) from the exponential distribution with mean \(1 / \lambda\) and set 

$$
X_{t} = \phi X_{t - 1} + \varepsilon_{t}
$$

 Step 3. Set \(t = t + 1\) and go to Step 2.

### Example 24: Example 22 (cont.)

The stock broker would typically have a large sample of data, but, for the sake of illustration, suppose that the following twenty time gaps between customer buy and sell orders had been recorded (in seconds): 1.95, 1.75, 1.58, 1.42, 1.28, 1.15, 1.04, 0.93, 0.84, 0.75, 0.68, 0.61, 11.98, 10.79, 9.71, 14.02, 12.62, 11.36, 10.22, 9.20. Standard calculations give \(\bar {X} = 5 . 2\) and \({\widehat {\sigma}} ^ {2} = 26 . 7\). To estimate the lag-1 autocorrelation, we need 

$$
\sum_{j = 1} ^ {19} X_{t} X_{t + 1} = 924. 1
$$

 Thus,\({\widehat {\mathsf {c o v}}} = [ 924 . 1 - ( 20 - 1 ) ( 5 . 2 ) ^ {2} ] / ( 20 - 1 ) = 21 . 6 ,\) and 

$$
\widehat {\rho} = \frac {21 . 6}{26 . 7} = 0. 8
$$

 Therefore, we could model the interarrival times as an EAR(1) process with \(\widehat {\lambda} = 1 / 5 . 2 = 0 . 192\) and \(\widehat {\phi} = 0 . 8 .\), provided that an exponential distribution is a good model for the individual gaps.

### 7.4 The Normal-to-Anything Transformation

The bivariate normal distribution and the AR(1) and EAR(1) time-series models are useful input models that are easy to fit and simulate. However, the marginal distribution is either normal or exponential, which is certainly not the best choice for many applications. Fortunately, we can start with a bivariate normal or AR(1) model and transform it to have any marginal distributions we want (including exponential).


Suppose we want to simulate a random variable X with cdf \(F ( x )\). Let \(Z\) be a standard normal random variable (mean 0 and variance 1), and let \(\Phi ( z )\) be its cdf. Then it can be shown that 

$$
R = \Phi (Z)
$$

 is a \(U ( 0 ,\)1) random variable. If we have a \(U ( 0 , 1 )\) random variable, we can get X by using the inverse cdf transformation 

$$
X = F^{- 1} [ R ] = F^{- 1} [ \Phi (Z) ]
$$

 We refer this as the normal to anything transformation, or NORTA for short.

Of course, if all we want is \(X ,\) then there is no reason to go to this trouble; we can just generate R directly. But suppose we want a bivariate random vector \(( X_{1} , X_{2} )\) such that \(X_{1}\) and \(X_{2}\) are correlated but their distributions are not normal. Then we can start with a bivariate normal random vector \(( Z_{1} , Z_{2} )\) and apply the NORTA transformation to obtain 

$$
X_{1} = F_{1} ^ {- 1} [ \Phi (Z_{1}) ] \text{and} X_{2} = F_{2} ^ {- 1} [ \Phi (Z_{2}) ]
$$

 There is not even a requirement that \(F_{1}\) and \(F_{2}\) be from the same distribution family; for instance,\(F_{1}\) could be an exponential distribution and \(F_{2}\) a beta distribution.

The same idea applies for time series. If \(Z_{t}\) is generated by an \(\mathbf{A R} ( 1 )\) with \(N ( 0 , 1 )\) marginals, then 

$$
X_{t} = F^{- 1} [ \Phi (Z_{t}) ]
$$

 will be a time-series model with marginal distribution \(F ( x )\). To ensure that \(Z_{t}\) is \(N ( 0 , 1 )\), we set \(\mu = 0\) and \(\sigma _ {\mathrm {{s}}} ^ {2} = 1 - \phi^{2}\) in the AR(1) model.

Although the NORTA method is very general, there are two technical issues that must be addressed to implement it:

1. The NORTA approach requires being able to evaluate that standard normal cdf,\(\Phi ( z )\), and the inverse cdf of the distributions of interest,\({\cal F} ^ {- 1} ( R )\). There is no closed-form expression for \(\Phi ( z )\) and no closed-form expression for \({\cal F} ^ {- 1} ( R )\) for many distributions. Therefore, numerical approximations are required. Fortunately, these functions are built into many symbolic calculation and spreadsheet programs, and we give one example next. In addition, Bratley, Fox, and Schrage [1987] contains algorithms for many distributions.

2. The correlation between the standard normal random variables \(( Z_{1} , Z_{2} )\) is distorted when it passes through the NORTA transformation. To be more specific, if \(( Z_{1} , Z_{2} )\)) have correlation \(\rho ,\) then in general \(X_{1} = F_{1} ^ {- 1} [ \Phi ( Z_{1} ) ]\) and \(X_{2} = F_{2} ^ {- 1} [ \Phi ( Z_{2} ) ]\)] will have a correlation \(\rho_{X} \neq \rho\) The difference is often small, but not always.

The second issue is more critical, because in input-modeling problems we want to specify the bivariate or lag-1 correlation. Thus, we need to find the bivariate normal correlation \(\rho\) that gives us the input correlation \(\rho_{X}\) that we want (recall that we specify the time-series model via the lag-1 correlation,\(\rho_{X} = \operatorname{c o r r} ( X_{t} , X_{t + 1} ) )\). There has been much research on this problem, including Cario and Nelson [1996, 1998] and Biller and Nelson [2003]. Fortunately, it has been shown that \(\rho_{X}\) is


a nondecreasing function of \(\rho ,\) and \(\rho\) and \(\rho_{X}\) will always have the same sign. Thus, we can do a relatively simple search based on the following algorithm:

Step 1. Set \(\rho = \rho_{X}\) to start.

Step 2. Generate a large number of bivariate normal pairs \(( Z_{1} , Z_{2} )\)) with correlation \(\rho ,\), and transform them into \(( X_{1} , X_{2} ) \mathbf{s}\), using the NORTA transformation.

Step 3. Compute the sample correlation between \(( X_{1} , X_{2} )\), using Equation (25), and call it \(\widehat {\rho} _ {T}\). If \(\widehat {\rho} _ {T} > \rho_{X}\), then reduce \(\rho\) and go to Step 2; if \(\widehat {\rho_{T}} < \rho_{X}\), then increase \(\rho\) and \({\bf g 0}\)to Step 2. If \(\widehat {\rho} _ {T} \approx \rho_{X}\) then stop.

### Example 25

Suppose we need \(X_{1}\) to have an exponential distribution with mean \(1 , X_{2}\) to have a beta distribution with \(\beta_{1} = 1 , \beta_{2} = 1 / 2\), and the two of them to have correlation \(\rho_{X} = 0 . 45\). Figure 8 shows a procedure in Maple that will estimate the required value of \(\rho .\). In the procedure, n is the number of sample pairs used to estimate the correlation. Running this procedure with n set to 1000 gives \(\rho = 0 . 52\)
## 8 Summary

Input-data collection and analysis require major time and resource commitments in a discrete-event simulation project. However, regardless of the validity or sophistication of the simulation model, unreliable inputs can lead to outputs whose subsequent interpretation could result in faulty recommendations.

This chapter discussed four steps in the development of models of input data: collecting the raw data, identifying the underlying statistical distribution, estimating the parameters, and testing for goodness of fit.

Some suggestions were given for facilitating the data-collection step. However, experience, such as that obtained by completing any of Exercises 1 through 5, will increase awareness of the difficulty of problems that can arise in data collection and of the need for planning.

Once the data have been collected, a statistical model should be hypothesized. Constructing a histogram is very useful at this point if sufficient data are available. A distribution based on the underlying process and on the shape of the histogram can usually be selected for further investigation.

The investigation proceeds with the estimation of parameters for the hypothesized distribution. Suggested estimators were given for distributions used often in simulation. In a number of instances, these are functions of the sample mean and sample variance.

The last step in the process is the testing of the distributional hypothesis. The \(q {-} q\) plot is a useful graphical method for assessing fit. The chi-square, Kolmogorov–Smirnov, and Anderson-Darling goodness-of-fit tests can be applied to many distributional assumptions. When a distributiona assumption is rejected, another distribution is tried. When all else fails, the empirical distribution could be used in the model


```txt
NORTARho := proc(rhoX, n)
local Z1, Z2, Ztemp, X1, X2, R1, R2, rho, rhoT, lower, upper;
randomize(123456);
Z1 := [random[normald[0,1]](n)]:
ZTemp := [random[normald[0,1]](n)]:
Z2 := [0]:
# set up bisection search
rho := rhoX:
if (rhoX < 0) then
    lower := -1:
    upper := 0:
else
    lower := 0:
    upper := 1:
fi:
Z2 := rho*Z1 + sqrt(1-rho^2)*ZTemp:
R1 := statevalf[cdf,normald[0,1]](Z1):
R2 := statevalf[cdf,normald[0,1]](Z2):
X1 := statevalf/icdf,exponential[1,0]](R1):
X2 := statevalf/icdf,beta[1,2]](R2):
rhoT := describe(linearcorrelation)(X1, X2);
# do bisection search until 5% relative error
while abs(rhoT - rhoX)/abs(rhoX) > 0.05 do
    if (rhoT > rhoX) then
    upper := rho:
    else
    lower := rho:
    fi:
    rho := evalf((lower + upper)/2):
    Z2 := rho*Z1 + sqrt(1-rho^2)*ZTemp:
    R1 := statevalf[cdf,normald[0,1]](Z1):
    R2 := statevalf[cdf,normald[0,1]](Z2):
    X1 := statevalf/icdf,exponential[1,0]](R1):
    X2 := statevalf/icdf,beta[1,2]](R2):
    rhoT := describe(linearcorrelation)(X1, X2);
end do;
RETURN(rho);
end;

Unfortunately, in some situations, a simulation study must be undertaken when there is not time or resources to collect data on which to base input models. When this happens, the analyst must use any available information—such as manufacturer specifications and expert opinion—to construct the input models. When input models are derived without the benefit of data, it is particularly importan to examine the sensitivity of the results to the models chosen.

Many, but not all, input processes can be represented as sequences of independent and identically distributed random variables. When inputs should exhibit dependence, then multivariate-input models are required. The bivariate normal distribution (and more generally the multivariate normal distribution) is often used to represent a finite number of dependent random variables. Time-series models are useful for representing a (conceptually infinite) sequence of dependent inputs. The NORTA transformation facilitates developing multivariate-input models with marginal distributions that are no normal.

REFERENCES

BILLER, B., AND B. L. NELSON [2003], “Modeling and Generating Multivariate Time Series with Arbitrary Marginals Using an Autoregressive Technique,” ACM Transactions on Modeling and Computer Simulation, Vol. 13, pp. 211–237.

BRATLEY, P., B. L. FOX, AND L. E. SCHRAGE [1987], A Guide to Simulation, 2d ed., Springer-Verlag, New York.

CARIO, M. C., AND B. L. NELSON [1996], “Autoregressive to Anything: Time-Series Input Processes for Simulation,” Operations Research Letters, Vol. 19, pp. 51–58.

CARIO, M. C., AND B. L. NELSON [1998], “Numerical Methods for Fitting and Simulating Autoregressiveto-Anything Processes,” INFORMS Journal on Computing, Vol. 10, pp. 72–81.

CHOI, S. C., AND R. WETTE [1969], “Maximum Likelihood Estimation of the Parameters of the Gamma Distribution and Their Bias,” Technometrics, Vol. 11, No. 4, pp. 683–890.

CHAMBERS, J. M., W. S. CLEVELAND, AND P. A. TUKEY [1983], Graphical Methodsfor Data Analysis, CRC Press, Boca Raton, FL.

CONNOVER, W. J. [1998], Practical Nonparametric Statistics, 3d ed., Wiley, New York.

DURBIN, J. [1975], “Kolmogorov–Smirnov Tests When Parameters Are Estimated with Applications to Tests of Exponentiality and Tests on Spacings,” Biometrika, Vol. 65, pp. 5–22.

FISHMAN, G. S. [1973], Concepts and Methods in Discrete Event Digital Simulation, Wiley, New York.

GUMBEL, E. J. [1943], “On the Reliability ofthe Classical Chi-squared Test,”Annals ofMathematical Statistics, Vol. 14, pp. 253ff.

HINES, W. W., D. C. MONTGOMERY, D. M. GOLDSMAN, AND C. M. BORROR [2002], Probability and Statistics in Engineering and Management Science, 4th ed., Wiley, New York.

JOHNSON, M. A., S. LEE, AND J. R. WILSON [1994], “NPPMLE and NPPSIM: Software for Estimating and Simulating Nonhomogeneous Poisson Processes Having Cyclic Behavior,” Operations Research Letters, Vol. 15, pp. 273–282.

JOHNSON, M. E. [1987], Multivariate Statistical Simulation, Wiley, New York.

LAW, A. M. [2007], Simulation Modeling and Analysis, 4th ed., McGraw-Hill, New York.

LILLIEFORS, H. W. [1967], “On the Kolmogorov–Smirnov Test for Normality with Mean and Variance Unknown,” Journal ofthe American Statistical Association, Vol. 62, pp. 339–402.

LILLIEFORS, H. W. [1969], “On the Kolmogorov–Smirnov Test for the Exponential Distribution with Mean Unknown,” Journal ofthe American Statistical Association, Vol. 64, pp. 387–389.

MANN, H. B., AND A. WALD [1942], “On the Choice of the Number of Intervals in the Application of the Chi-squared Test,” Annals ofMathematical Statistics, Vol. 18, p. 50ff.

NELSON, B. L.,AND M.YAMNITSKY [1998], “Input Modeling Tools for Complex Problems,” in Proceedings of the 1998 Winter Simulation Conference, D. Medeiros, E. Watson, J. Carson, and M. Manivannan, eds., Washington, DC, Dec. 13–16, pp. 105–112.

STUART, A., J. K. ORD, AND E. ARNOLD [1998], Kendall’s Advanced Theory of Statistics, 6th ed., Vol. 2, Oxford University Press, Oxford

EXERCISES

1. Go to a small store and record the interarrival- and service-time distributions. If there are severa workers, how do the service-time distributions compare to one another? Do service-time distributions need to be constructed for each type of appliance? (Make sure that the management gives permission to perform this study.)

2. Go to a cafeteria and collect data on the distributions of interarrival and service times. The distribution of interarrival times is probably different for each of the three daily meals and migh also vary during the meal—that is, the interarrival time distribution for 11:00 A.M. to 12:00 noon could be different from that for 12:00 noon to 1:00 P.M. Define service time as the time from when the customer reaches the point at which the first selection could be made until the time of exiting from the cafeteria line. (Any reasonable modification of this definition is acceptable.) The service-time distribution probably changes for each meal. Can times of the day or days of the week for either distribution be grouped to suit the homogeneity of the data? (Make sure tha the management gives permission to perform this study.)

3. Go to a major traffic intersection and record the interarrival-time distributions from each direction. Some arrivals want to go straight, some turn left, some turn right. The interarrival-time distribution varies during the day and by day of the week. Every now and then an acciden occurs. Develop an input model for your data.

4. Go to a grocery store and construct the interarrival and service distributions at the checkout counters. These distributions might vary by time of day and by day of week. Also record the number of service channels available at all times. (Make sure that the management gives permission to perform this study.)

5. Go to a laundromat and collect data to develop input models for the number of washers and dryers a customer uses (probably dependent). Also collect data on the arrival rate of customers (probably not stationary). (Make sure that the management gives permission to perform this study.)

6. Prepare four theoretical normal density functions, all on the same figure, each distribution having mean zero, but let the standard deviations be 1/4, 1/2, 1, and 2.

7. On one figure, draw the pdfs of the Erlang distribution where $\theta = 1 / 2$ and $k = 1 , 2 , 4$, and 8.

8. On one figure, draw the pdfs of the Erlang distribution where θ 2 and k 1, 2, 4, and 8.

9. Draw the pmf of the Poisson distribution that results when the parameter α is equal to the following:

10. On one figure, draw the two exponential pdfs that result when the parameter λ equals 0.6 and 1.2.

11. On one figure, draw the three Weibull pdfs that result when ν 0,$\alpha = 1 / 2 ,$, and $\beta = 1 , 2 ,$, and 4.

12. The following data are generated randomly from a gamma distribution:

Compute the maximum-likelihood estimators $\widehat {\beta}$ and ${\widehat {\theta}} .$13. The following data are generated randomly from a Weibull distribution where $\nu = 0 \mathrm{:}$| 7.936 | 5.224 | 3.937 | 6.513 | | — | — | — | — | | 4.599 | 7.563 | 7.172 | 5.132 | | 5.259 | 2.759 | 4.278 | 2.696 | | 6.212 | 2.407 | 1.857 | 5.002 | | 4.612 | 2.003 | 6.908 | 3.326 |

Compute the maximum-likelihood estimators α and ${\widehat {\beta}} .$(This exercise requires a programmable calculator, a computer, or a lot of patience.)

14. The highway between Atlanta, Georgia and Athens, Georgia has a high incidence of accidents along its 100 kilometers. Public safety officers say that the occurrence of accidents along the

highway is randomly (uniformly) distributed, but the news media say otherwise. The Georgia Department of Public Safety published records for the month of September. These records indicated the point at which 30 accidents involving an injury or death occurred, as follows (the data points representing the distance from the city limits of Atlanta):

Use the Kolmogorov–Smirnov test to discover whether the distribution of location of accidents is uniformly distributed for the month of September.

15. Show that the Kolmogorov–Smirnov test statistic for Example 19 is $D = 0 . 1054$16. Records pertaining to the monthly number of job-related injuries at an underground coal mine were being studied by a federal agency. The values for the past 100 months were as follows:

(a) Apply the chi-square test to these data to test the hypothesis that the underlying distribution is Poisson. Use the level of significance $\alpha = 0 . 05$(b) Apply the chi-square test to these data to test the hypothesis that the distribution is Poisson with mean 1.0. Again let $\alpha = 0 . 05 .$(c) What are the differences between parts (a) and (b), and when might each case arise?

17. The time required for 50 different employees to compute and record the number of hours worked during the week was measured, with the following results in minutes:

Use the chi-square test (as in Example 18) to test the hypothesis that these service times are exponentially distributed. Let the number of class intervals be k 6. Use the level of significance α 0.05.

Using input-modeling software, apply as many tests for normality as are available in the software. If the chi-square test is available, apply it with at least two different choices for the number of intervals. Do all of the tests reach the same conclusion?

19. The Crosstowner was a bus that cut a diagonal path from northeast Atlanta to southwest Atlanta. The time required to complete the route was recorded by the bus operator. The bus runs from Monday through Friday. The times of the last fifty 8:00 A.M. runs, in minutes, are as follows:

How are these run times distributed? Develop and test a suitable model.

20. The time required for the transmission of a message (in milliseconds) is sampled electronically at a communications center. The last 50 values in the sample are as follows:

How are the transmission times distributed? Develop and test an appropriate model.

21. The time (in minutes) between requests for the hookup ofelectric service was accurately recorded at the Gotwatts Flash and Flicker Company, with the following results for the last 50 requests:

How are the times between requests for service distributed? Develop and test a suitable model.

22. Daily demands for transmission overhaul kits for the D-3 dragline were maintained by Earth Moving Tractor Company, with the following results:

How are the daily demands distributed? Develop and test an appropriate model.

23. A simulation is to be conducted of a job shop that performs two operations: milling and planing, in that order. It would be possible to collect data about processing times for each operation, then generate random occurrences from each distribution. However, the shop manager says that the times might be related; large milling jobs take lots of planing. Data are collected for the next 25 orders, with the following results in minutes:

(a) Plot milling time on the horizontal axis and planing time on the vertical axis. Do these data seem dependent? (b) Compute the sample correlation between milling time and planing time. (c) Fit a bivariate normal distribution to these data.

24. Write a computer program to compute the maximum-likelihood estimators $( {\widehat {\alpha}} , {\widehat {\beta}} )$ of the Weibull distribution. Inputs to the program should include the sample size $n ;$the observations $x_{1} , x_{2} , \ldots , x_{n} ;$ a stopping criterion  (stop when $| f ( \widehat {\beta} _ {j} ) | \leq \epsilon ) $; and a print option OPT (usually set = 0). Output would be the estimates α and ${\widehat {\beta}} . {\mathrm{I f}} {\mathrm{O P T}} = 1$, additional output would be printed, as in Table 4, showing convergence. Make the program as “user friendly” as possible.

25. Examine a computer-software library or simulation-support environment to which you have access. Obtain documentation on data-analysis software that would be useful in solving exercises 7 through 24. Use the software as an aid in solving selected problems.

26. The numbers of patrons staying at a small hotel on 20 successive nights was observed to be 20, 14, 21, 19, 14, 18, 21, 25, 27, 26, 22, 18, 13, 18, 18, 18, 25, 23, 20, 21. Fit both an AR(1) and an EAR(1) model to these data. Decide which model provides a better fit by looking at a histogram of the data.

27. The following data represent the time to perform transactions in a bank, measured in minutes: 0.740, 1.28, 1.46, 2.36, 0.354, 0.750, 0.912, 4.44, 0.114, 3.08, 3.24, 1.10, 1.59, 1.47, 1.17, 1.27, 9.12, 11.5, 2.42, 1.77. Develop an input model for these data.

28. Two types ofjobs (A and B) are released to the input buffer of a job shop as orders arrive, and the arrival of orders is uncertain. The following data are available from the last week of production:

Develop an input model for the number of new arrivals of each type each day.

29. The following data are available on the processing time at a machine (in minutes): 0.64, 0.59, 1.1, 3.3, 0.54, 0.04, 0.45, 0.25, 4.4, 2.7, 2.4, 1.1, 3.6, 0.61, 0.20, 1.0, 0.27, 1.7, 0.04, 0.34. Develop an input model for the processing time.

30. For the arrival-count data shown here, estimate the arrival-rate function λ(t) for an NSPP in arrivals per hour:

31. Using the web, research some of the input-modeling software packages mentioned in this chapter (try www.informs-sim.org as a starting place). What are their features? What distributions do they include?

32. As part of a surgical admitting process orderlies are available to escort patients to the presurgery preparation area. Patients are not allowed to go by themselves as a matter of policy. A Management Engineer walked the route briskly and found it took 2 minutes to walk to the prep area. Orderlies say that it sometimes takes as long as 7 minutes to escort a patient to the area, but it usually takes 5 minutes. Construct an input model for time to escort a patient based on the information given.

33. At an insertion/sealing/inspection operation electronic components are inserted into a board, the case is assembled and sealed, and the completed unit is tested. The time to accomplish all of this is expected to be 2 minutes for type A12117c assemblies and 3 minutes for the B33433x assemblies. An average deviation of plus-or-minus 20% around these values is considered the standard allowance for variability. Construct an input model for the insertion/sealing/inspection operation based on the information given.

34. The following data are counts of calls arriving to a call center by hour of the day over five Mondays. Does it appear that the process is nonstationary? If so, then develop an NSPP inpu model for the arrival of calls between 8 A.M. and 4 P.M.

This page intentionally left blank