Crafts of Simulation Programming

The following sections are included:

• Dedication
• Preface
• Contents

The following sections are included:

• Introduction
• Random Numbers Generators
  • Basic Generators
  • The Need for Multiple Substreams
  • Computing (a × s) mod m
  • Computing the Jumping Matrices
  • A Random Number Package
  • Jumping Backward
• Examples of Using Random Number Generator
• Nonuniform Random Variates
  • Random Variates of Various Distributions
  • Correlated Random Variates
• Utilities
  • Numerical Approximation of Normal Distribution
  • Quantile of Normal Distribution
  • Quantile of t Distribution
  • Quantile of Chi-square Distribution
  • Standard Deviation
• Summary

Simulation can be generally classified into two types: 1) Finite-Horizon (terminating) simulation: a simulation where there is a specific starting and stopping condition that is part of the model; and 2) Steady-State (non-terminating) simulation: a simulation where there is no specific starting and ending conditions. The purpose of a steady-state simulation is the study of the long-run behavior of a system. A performance measure is called a steady-state parameter if it is a characteristic of the equilibrium distribution of an output stochastic process. Steady-state distribution does not depend on initial conditions (when the initialization bias have been properly taken care of), but the nature and rate of convergence of the transient distributions can depend heavily on the initial conditions. The length of the simulation must be determined dynamically to ensure that a representative steady-state model behavior is reached. We seek algorithms that provide good results and do not affect the performance of a given system.

Many simulation procedures are derived based on the assumption that data are i.i.d. normal; for example, ranking and selection procedures and multiple-comparison procedures. Batch means are the method of choice to manufacture data that are approximately i.i.d. normal. In this chapter, we discuss a procedure to manufacture batch means that appear to be i.i.d. normal, as determined by the von Neumman test of independence and chi-square test of normality. The procedure performs statistical analysis on sample sequences collected on strictly stationary stochastic processes. It determines the batch sizes based entirely on data and does not require any user intervention. The only required condition is that the autocorrelations of the stochastic process output sequence die off as the lag between observations increases, in the sense of ϕ-mixing.

Order statistics are among the most fundamental tools in non-parametric statistics and inference and have been widely studied and applied to many real-world issues. It is often of interest to estimate the reliability of the component/system from the observed lifetime data. In many applications we want to estimate the probability that a future observation will exceed a given high level during some specified epoch. For example, a machine may break down if a certain temperature is reached.

In this chapter, we investigate an important special case of order statistics, namely, the underlying variates are correlated normal random variables and show how it is related to ranking and selection. [Gupta et al. (1973)] point out that the study of these order statistics is important because that the joint distributions of standardized correlated variables often tend asymptotically to the multivariate normal distribution. [Arellano and Genton (2007, 2008)] have derived the density of correlated normal random variables from a multivariate normal distribution. [Liu et al. (2012)] develop a method to construct exact simultaneous confidence interval for a finite set of contrasts of three, four or five generally correlated normal means.

Simulation is often used to investigate system characteristics, i.e., the distributional properties of an output statistic calculated from the simulation's results, such as the mean or variance of system performance measures. In many applications we want to estimate the probability that a future observation will exceed a given level during some specified epoch. For example, a machine may break down if a certain temperature is reached. Among many other aspects, reliability analysis (risk analysis) studies the expected life and the failure rate of a component or a system of components linked together in some structure. It is often of interest to estimate the reliability of the component/system from the observed lifetime data. With probability modeling, normal and lognormal densities are commonly used to model certain lifetimes in reliability and survival analysis as well as risk management.

This chapter reviews metamodels, investigates the accuracy of the fitted quantile curves (response surface) and how to construct non-functional-form metamodels with a set of histograms at certain design point.

Simulation studies have been used to investigate the characteristics of systems, for example the mean and the variance of certain system performance measures like waiting times in queue. However, without knowledge of the underlying distribution, the mean and the variance provide only limited information. On the other hand, the pdf f gives a natural description of the distribution of a stationary continuous output random variable X produced by a simulation and reveals many characteristics of the underlying distributions beyond just the mean and variance (for instance, tail probabilities and quantiles).

This chapter investigates the (Behrens-Fisher) problem of interval estimation and hypothesis testing concerning the difference between the means of two normally distributed populations when the variances of the two populations maybe unequal, based on two independent samples.

Among many other aspects, reliability analysis studies the expected life and the failure rate of a component or a system of components linked together in some structure. We propose applying ranking-and-selection procedures to this analysis. That is, we are more interested in whether a given component is better than the others rather than the accuracy of the performance measures. The underling philosophy is to rank estimators through ordinal comparison while the precision of the estimates are still poor [Ho et al. (1992)], hence, increase the efficiency of reliability analysis. A detailed analysis to obtain more accurate performance measures can then be carried out on the few chosen systems.

One crucial element in R&S procedures is the event of “correct selection” (CS) of the true best system. In a stochastic simulation, the possibility of CS, denoted by P(CS), increases as the sample sizes become larger. Most indifference-zone-selection procedures are directly or indirectly based on the procedures of [Dudewicz and Dalal (1975)] and [Rinott (1978)]. However, these procedures determine the number of additional replications based on a conservative LFC assumption and do not take into account the value of sample means. If the accuracy requirement is high and the total number of designs in a decision problem is large, then the total simulation cost can easily become prohibitively high.

The common random numbers (CRN) simulation technique is a variance-reduction method in which policy alternatives are tested against the same random input streams. The CRN literature suggests that positively correlated input streams will generate positively correlated policy responses and, therefore, that the variance of CRN estimators of response differences will be smaller than the variance of independent sample estimators. However, because the assumption of independent samples across systems is used to develop the two-stage selection procedures (see, e.g., [Dudewicz and Dalal (1975)]), CRN are generally not used with those procedures for variance reduction. There are other variance-reduction techniques, e.g., antithetic variates, control variates, see [Law (2014)] for details.

Parallel and distributed simulation (PADS) studies how a network of several interconnected models work together to support decision making by distributing the execution of a discrete event simulation (DES) program over multiple computers. Parallel DES programs are executed on multiprocessor computing platforms containing multiple central processing units that interact frequently. Distributed DES programs are executed on loosely coupled systems that may be geographically distributed and require additional interaction times. However, with new computing paradigms such as clusters of workstations and grid computing, the distinction has become un-noticeable. In both cases the execution of a single simulation model, likely composed of several simulation programs, is distributed over multiple processors (computers) and can be executed concurrently. Hence, one can reduce the execution time by up to a factor equal to the number of processors that are used. Distributing the execution across multiple computers and utilizing the resources of many computer systems is also beneficial in allowing the execution of larger simulations where the capacity of one computer may not be enough to carry out the simulation. A more detailed discussion of distributed and parallel simulation can be found in [Fujimoto (2000)].

The ranking and selection studies in previous chapters have focused on selection based on a single measure of system performance. In many practical situations, however, we need to select systems based on multiple criteria (attributes, objectives, or performance measures). For example, in product-design optimization, the cost and the quality of products are two conflicting objectives. In evaluating airline flight schedules, we may want to select flight schedules in terms of minimal (flight and ground) staffs and minimal percentage of late arrivals. In this setting, the problem of selecting the best systems from a set of alternatives through simulation becomes a multi-objective ranking and selection problem. Researchers have developed procedures to address this problem. [Butler et al. (2010)] combine the multiple attribute utility theory with the [Rinott (1978)] procedure to handle multiple performance measures. Before performing selection, users transform multiple performance measures into one utility score. The procedure then finds the system that gives the highest utility. [Dudewicz and Taneja (1978)] propose a multivariate procedure by defining a multivariate normal vector composed of ω > 1 component variates with an unknown and unequal variance-covariance matrix. They redefine the indifference-zone parameter as the Euclidean distance from a mean vector to the best mean vector.

[Andradóttir and Kim (2010)] develop procedures to find the system with the best primary performance measure in the presence of a stochastic constraint on a secondary performance measure, i.e., selection with constraints. Their fully sequential procedures increase one sample for each system that is still under consideration at each iteration and eliminate systems from further simulation when the partial sum is greater than a threshold. These procedures are efficient in terms of sample sizes but not necessarily in terms of runtime. [Morrice and Butler (2006)] extend the multiple-attribute-utility theorem [Butler et al. (2010)] to perform selection with constraints. However, in some cases, it is not straightforward to assign utility scores.

The following sections are included:

• Appendix A: Tables of Critical Constants
• Bibliography
• Index