As a good rule of thumb, the maximum lag for which autocorrelations are computed should be approximately 2% of the number of n realizations, although each rj,x could be tested to determine if it is significantly different from zero. Show
Sample Size Determination: We can calculate the minimum sample size required by n = [1 + 2A ] S2 t2 / (d2 mean2) Application: A pilot run was made of a model, observations numbered 150, the mean was 205.74 minutes and the variance S2 = 101, 921.54, estimate of the lag coefficients were computed as: r1,x = 0.3301 r2,x = 0.2993, and r3,x = 0.1987. Calculate the minimum sample size to assure the estimate lies within +d = 10% of the true mean with a = 0.05. n = [(1.96)2 (101,921.54) {1 + 2 [(1-1/4) 0.3301 + (1 - 2/4) 0.2993 + (1- 3/4) 0.1987]}] / (0.1)2 (205.74)2� 1757 You may like using Statistics for Time Series, and Testing Correlation JavaScript. What Is Central Limit Theorem?For practical purposes, the main idea of the central limit theorem (CLT) is that the average of a sample of observations drawn from some population with any shape-distribution is approximately distributed as a normal distribution if certain conditions are met. In theoretical statistics there are several versions of the central limit theorem depending on how these conditions are specified. These are concerned with the types of assumptions made about the distribution of the parent population (population from which the sample is drawn) and the actual sampling procedure.One of the simplest versions of the theorem says that if is a random sample of size n (say, n larger than 30) from an infinite population, finite standard deviation , then the standardized sample mean converges to a standard normal distribution or, equivalently, the sample mean approaches a normal distribution with mean equal to the population mean and standard deviation equal to standard deviation of the population divided by the square root of sample size n. In applications of the central limit theorem to practical problems in statistical inference, however, statisticians are more interested in how closely the approximate distribution of the sample mean follows a normal distribution for finite sample sizes, than the limiting distribution itself. Sufficiently close agreement with a normal distribution allows statisticians to use normal theory for making inferences about population parameters (such as the mean ) using the sample mean, irrespective of the actual form of the parent population. It is well known that whatever the parent population is, the standardized variable will have a distribution with a mean 0 and standard deviation 1 under random sampling. Moreover, if the parent population is normal, then it is distributed exactly as a standard normal variable for any positive integer n. The central limit theorem states the remarkable result that, even when the parent population is non-normal, the standardized variable is approximately normal if the sample size is large enough (say > 30). It is generally not possible to state conditions under which the approximation given by the central limit theorem works and what sample sizes are needed before the approximation becomes good enough. As a general guideline, statisticians have used the prescription that if the parent distribution is symmetric and relatively short-tailed, then the sample mean reaches approximate normality for smaller samples than if the parent population is skewed or long-tailed. In this lesson, we will study the behavior of the mean of samples of different sizes drawn from a variety of parent populations. Examining sampling distributions of sample means computed from samples of different sizes drawn from a variety of distributions, allow us to gain some insight into the behavior of the sample mean under those specific conditions as well as examine the validity of the guidelines mentioned above for using the central limit theorem in practice. Under certain conditions, in large samples, the sampling distribution of the sample mean can be approximated by a normal distribution. The sample size needed for the approximation to be adequate depends strongly on the shape of the parent distribution. Symmetry (or lack thereof) is particularly important. For a symmetric parent distribution, even if very different from the shape of a normal distribution, an adequate approximation can be obtained with small samples (e.g., 10 or 12 for the uniform distribution). For symmetric short-tailed parent distributions, the sample mean reaches approximate normality for smaller samples than if the parent population is skewed and long-tailed. In some extreme cases (e.g. binomial) samples sizes far exceeding the typical guidelines (e.g., 30) are needed for an adequate approximation. For some distributions without first and second moments (e.g., Cauchy), the central limit theorem does not hold. '' Central Limit Theorem Justification '' Preamble Define X, XBAR as a real 1-dimensional arrays Define I, J, M, N as integer variables End Main Open 3 for output, Name = "CLT.OUT" Use 3 for output LET N=50 LET M = 1000 Reserve X(*) as N Reserves XBAR(*) as M For J = 1 to M DO For I = 1 to N DO X(I) = Uniform.f(1., 0., 1) Compute XBAR (j) as the average of X(i) Loop Compute Ave as the average of XBAR(j) Compute ST as the standard deviation of XBAR(j) Loop LL = Ave -1.96*ST/SQRT.F(real.f(M)) UL = Ave + 1.96*ST/SQRT.F(real.f(M)) For J =1 to M Do IF XBAR(j) < LL AND XBAR(j) > XU Let Count = Count + 1 Always Loop Print 1 line with Count/M thus The P-value is = ***** END What Is a Least Squares Model?Many problems in analyzing data involve describing how variables are related. The simplest of all models describing the relationship between two variables is a linear, or straight-line, model. The simplest method of fitting a linear model is to "eye-ball'' a line through the data on a plot. A more elegant, and conventional method is that of "least squares", which finds the line minimizing the sum of distances between observed points and the fitted line.Realize that fitting the "best'' line by eye is difficult, especially when there is a lot of residual variability in the data. Know that there is a simple connection between the numerical coefficients in the regression equation and the slope and intercept of regression line. Know that a single summary statistic like a correlation coefficient does not tell the whole story. A scatter plot is an essential complement to examining the relationship between the two variables. ANOVA: Analysis of VarianceThe tests we have learned up to this point allow us to test hypotheses that examine the difference between only two means. Analysis of Variance or ANOVA will allow us to test the difference between 2 or more means. ANOVA does this by examining the ratio of variability between two conditions and variability within each condition. For example, say we give a drug that we believe will improve memory to a group of people and give a placebo to another group of people. We might measure memory performance by the number of words recalled from a list we ask everyone to memorize. A t-test would compare the likelihood of observing the difference in the mean number of words recalled for each group. An ANOVA test, on the other hand, would compare the variability that we observe between the two conditions to the variability observed within each condition. Recall that we measure variability as the sum of the difference of each score from the mean. When we actually calculate an ANOVA we will use a short-cut formulaThus, when the variability that we predict (between the two groups) is much greater than the variability we don't predict (within each group) then we will conclude that our treatments produce different results. Exponential Density FunctionAn important class of decision problems under uncertainty concerns the chance between events. For example, the chance of the length of time to next breakdown of a machine not exceeding a certain time, such as the copying machine in your office not to break during this week.Exponential distribution gives distribution of time between independent events occurring at a constant rate. Its density function is: f(t) = l exp(-lt), where l is the average number of events per unit of time, which is a positive number. The mean and the variance of the random variable t (time between events) are 1/ l, and 1/l2, respectively. Applications include probabilistic assessment of the time between arrival of patients to the emergency room of a hospital, and arrival of ships to a particular port. Comments: Special case of both Weibull and gamma distributions. You may like using Exponential Applet to perform your computations. You may like using the following Lilliefors Test for Exponentially to perform the goodness-of-fit test. Poisson ProcessAn important class of decision problems under uncertainty is characterized by the small chance of the occurrence of a particular event, such as an accident. Gives probability of exactly x independent occurrences during a given period of time if events take place independently and at a constant rate. May also represent number of occurrences over constant areas or volumes. The following statements describe the Poisson Process:
Poisson process are often used, for example in quality control, reliability, insurance claim, incoming number of telephone calls, and queuing theory. An Application: One of the most useful applications of the Poisson Process is in the field of queuing theory. In many situations where queues occur it has been shown that the number of people joining the queue in a given time period follows the Poisson model. For example, if the rate of arrivals to an emergency room is l per unit of time period (say 1 hr), then: P ( n arrivals) = ln e-l / n! The mean and variance of random variable n are both l . However if the mean and variance of a random variable having equal numerical values, then it is not necessary that its distribution is a Poisson. Applications: P ( 0 arrival) = e-l P ( 1 arrival) = l e-l / 1! P ( 2 arrival) = l2 e-l / 2! and so on. In general: P ( n+1 arrivals ) = l Pr ( n arrivals ) / n. You may like using Poisson Applet to perform your computations. Goodness-of-Fit for PoissonReplace the numerical example data with your up-to-14 pairs of Observed values & their frequencies, and then click the Calculate button. Blank boxes are not included in the calculations. In entering your data to move from cell to cell in the data-matrix use the Tab key not arrow or enter keys. For Technical Details, Back to: Statistical Thinking for Decision Making Uniform Density FunctionApplication: Gives probability that observation will occur within a particular interval when probability of occurrence within that interval is directly proportional to interval length.Example: Used to generate random numbers in sampling and Monte Carlo simulation. Comments: Special case of beta distribution. The mass function of geometric mean of n independent uniforms [0,1] is: P(X = x) = n x(n - 1) (Log[1/xn])(n -1) / (n - 1)!. zL = [UL-(1-U)L] / L is said to have Tukey's symmetrical l-distribution. You may like using Uniform Applet to perform your computations. Test for Binomial: NPAR TEST BINOMIAL(p)=GENDER(0, 1) Gooness-of-fit for discrete r.v.: NPAR TEST CHISQUARE=X (1,3)/EXPECTED=20 30 50 Needed information to perform the t-test: DISCRIPTIVES X /STATISTICS= 1 2 Two population t-test T-TEST GROUPS=GENDER(1,2)/VARIABLES=X Plot x vs y: PLOT FORMAT=REGRESSION/SYMBOLS='*' /TITLE='PLOT OF Y ON X' /VERTICAL='Y' /HORIZONTAL='X' /PLOT=Y WITH X RANDOM VARIATES GENERATORS: LOOP #I = 1 to 100. (normal with mean = 0 and std = 1) COMPUTE XNORM = RV.NORMAL(0,1) (chi-square with 2 d.f.) COMPUTE XCHISQ=RV.CHISQ(2) (exponential with mean 2) COMPUTE XEXPON = RV.EXP(1/2) (binomial n = 10 and p = .50 COMPUTE XBINOM = RV.BINOM(10, 0.5). END CASE END LOOP END FILE END INOUT PROGRAM EXAMINE VARS=ALL /STATISTICS /HISTOGRAM(NORMSAL)= XNORM /HISTOGRAM(NORMSAL)= XCHISQ /HISTOGRAM(NORMSAL)= XEXPON /HISTOGRAM(NORMSAL)= XBINOM NORMAL RANOM VARIATE GENERATOR, K-S, AND RUNS TESTS: SPSS/OUTPUT=HW3.OUT TITLE 'GENERATING FROM NORMAL 0,1' INPUT PROGRAM LOOP I=1 TO 50 COMPUTE X2=NORMAL(1) END CASE END LOOP END FILE END INPUT PROGRAM VAR LABLE X2 'NORMAL VARIATE' LIST CASE CASE=50/VARIABLE=ALL// CONDESCRIPTIVE X2(ZX2) STATISTICS ALL FREQUENCIES VARIABLE=ZX2/FORMAT=NOTABLE/ HISTOGRAM MIN(-3.0) MAX(+3.0) INCREMENT(0.2)/ NPAR TESTS RUNS(MEAN)=ZX2/ NPAR TESTS K-S(NORMAL,0.0,1.0)=ZX2/ SAMPLE 10 FROM 50 LIST CASE CASE=10/VARIABLES=X2,ZX2/ FINISH K-S LILLIEFORS TEST FOR NORMALITY: $SPSS/OUTPUT=L.OUT TITLE 'K-S LILLIEFORS TEST FOR NORMALITY' DATA LIST FREE FILE='L.DAT'/X VAR LABELS X 'SAMPLE VALUES' LIST CASE CASE=20/VARIABLES=ALL CONDESCRIPTIVE X(ZX) LIST CASE CASE=20/VARIABLES=X ZX/ SORT CASES BY ZX(A) RANK VARIABLES=ZX/RFRACTION INTO CRANK/TIES=HIGH COMPUTE Y=CDFNORM(ZX) COMPUTE SPROB=CRANK COMPUTE DA=Y-SPROB COMPUTE DB=Y-LAG(SPROB,1) COMPUTE DAABS=ABS(DA) COMPUTE DBABS=ABS(DB) COMPUTE LILLSTAT=MAX(DAABS,DBABS) LIST VARIABLES=X,ZX,Y,SPROB,DA,DB LIST VARIABLES=LILLSTAT SORT CASES BY LILLSTAT(D) LIST CASES CASE=1/VARIABLES=LILLSTAT FINISH K-S TEST FOR EXPONENTIAL DATA WITH MEAN = 1 DATA LIST FREE FILE='Ex.DAT'/ X DESCRIPTIVES VARIABLES=X STATISTICS MEAN COMPUTE Y=1.-EXP(-X) NPAR TESTS K-S(UNIFORM, 0.0, 1.0) Y (For large sample size) FINISH For more SPSS programs useful to simulation input/output analysis, visit Data Analysis Routines. Classical uniform random number generators have some major defects, such as, short period length and lack of higher dimension uniformity. However, nowadays there are a class of rather complex generators which is as efficient as the classical generators while enjoy the property of a much longer period and of a higher dimension uniformity. Computer programs that generate "random" numbers use an algorithm. That means if you know the algorithm and the seedvalues you can predict what numbers will result. Because you can predict the numbers they are not truly random - they are pseudorandom. For statistical purposes "good" pseudorandom numbers generators are good enough. real function random() c c Algorithm AS 183 Appl. Statist. (1982) vol.31, no.2 c c Returns a pseudo-random numbers with rectangular distribution. c between 0 and 1. The cycle length is 6.95E+12 (See page 123 c of Applied Statistics (1984) vol.33), not as claimed in the c original article. c c IX, IY and IZ should be set to integer values between 1 and c 30000 before the first entry. c c Integer arithmetic up to 30323 is required. c integer ix, iy, iz common /randc/ ix, iy, iz c ix = 171 * mod(ix, 177) - 2 * (ix / 177) iy = 172 * mod(iy, 176) - 35 * (iy / 176) iz = 170 * mod(iz, 178) - 63 * (iz / 178) c if (ix .lt. 0) ix = ix + 30269 if (iy .lt. 0) iy = iy + 30307 if (iz .lt. 0) iz = iz + 30323 c c If integer arithmetic up to 5212632 is available, the preceding c 6 statements may be replaced by: c c ix = mod(171 * ix, 30269) c iy = mod(172 * iy, 30307) c iz = mod(170 * iz, 30323) c random = mod(float(ix) / 30269. + float(iy) / 30307. + + float(iz) / 30323., 1.0) return end c c c c real function uniform() c c Generate uniformly distributed random numbers using the 32-bit c generator from figure 3 of: L'Ecuyer, P., 1988. c The cycle length is claimed to be 2.30584E+18 c Seeds can be set by calling the routine set_uniform c It is assumed that the Fortran compiler supports long variable c names, and integer*4. c integer*4 z, k, s1, s2 common /unif_seeds/ s1, s2 save /unif_seeds/ c k = s1 / 53668 s1 = 40014 * (s1 - k * 53668) - k * 12211 if (s1 .lt. 0) s1 = s1 + 2147483563 c k = s2 / 52774 s2 = 40692 * (s2 - k * 52774) - k * 3791 if (s2 .lt. 0) s2 = s2 + 2147483399 c z = s1 - s2 if (z .lt. 1) z = z + 2147483562 c uniform = z / 2147483563. return end subroutine set_uniform(seed1, seed2) c c Set seeds for the uniform random number generator. c integer*4 s1, s2, seed1, seed2 common /unif_seeds/ s1, s2 save /unif_seeds/ s1 = seed1 s2 = seed2 return end The Random Number Generator RANECU A FORTRAN code for a generator of uniform random numbers on [0,1]. RANECU is multiplicative linear congruential generator suitable for a 16-bit platform. It combines three simple generators, and has a period exceeding 81012. It is constructed for more efficient use by providing for a sequence of such numbers, LEN in total, to be returned in a single call. A set of three non-zero integer seeds can be supplied, failing which a default set is employed. If supplied, these three seeds, in order, should lie in the ranges [1,32362], [1,31726] and [1,31656] respectively. SUBROUTINE RANECU (RVEC,LEN) C Portable random number generator for 16 bit computer. C Generates a sequence of LEN pseudo-random numbers, returned in C RVEC. DIMENSION RVEC(*) SAVE ISEED1,ISEED2, ISEED3 DATA ISEED1,ISEED2,ISEED3/1234, 5678, 9876/ C Default values, used if none suppliedvia an ENTRY C call at RECUIN DO 100 I = 1,LEN K=ISEED1/206 ISEED1 = 157 * (ISEED1 - K * 206) - K * 21 IF(ISEED1.LT.0) ISEED1=ISEED1+32363 K=ISEED2/217 ISEED2 = 146 * (ISEED2 - K*217) - K* 45 IF(ISEED2.LT.O) ISEED2=ISEED2+31727 K=ISEED3/222 ISEED3 = 142 * (ISEED3 - K *222) - K * 133 IF(ISEED3.LT.0) ISEED3=ISEED3+31657 IZ=ISEED1-ISEED2 IF(IZ.GT.706)IZ = Z - 32362 IZ = 1Z+ISEED3 IF(IZ.LT.1)IZ = 1Z + 32362 RVEC(I)=REAL(IZ) * 3.0899E - 5 100 CONTINUE RETURN ENTRY RECUIN(IS1, IS2, IS3) ISEED1=IS1 ISEED2=IS2 ISEED3=IS3 RETURN ENTRY RECUUT(IS1,IS2,IS3) IS1=ISEED1 IS2=ISEED2 IS3=ISEED3 RETURN END The Shuffling Routine in Visual Basic The Square Histogram Method A simple such histogram, layed flat, might be: a:********************************32 b:**************************27 c:************************26 d:************12 e:***3 The idea is to cut the bars into pieces then reassemble them into a square histogram, all heights equal, with each final bar having a lower part, as well as an upper part indicating where it came from. A single uniform random variable U can then be used to choose one of the final bars and to indicate whether to use the lower or upper part. There are many ways to do this cutting and reassembling; the simplest seems to be the Robin Hood Algorithm: Take from richest to bring the poorest up to average. STEP 1: The original (horizontal) histogram, average "height" 20: a:******************************** 32 b:************************** 27 c:************************ 26 d:************ 12 e:*** 3 Take 17 from strip 'a' to bring strip 'e' up to average. Record donor and use old 'poor' level to mark lower part of donee: a:*************** 15 b:************************** 27 c:************************ 26 d:************ 12 e:**|****************| 20 (a) Then bring 'd' up to average with donor 'b'. Record donor and use old 'poor' level to mark lower part of donee: a:*************** 15 b:****************** 19 c:************************ 26 d:***********|*******| 20 (b) e:*****|*************| 20 (a) Then bring 'a' up to average with donor 'c'. Record donor and use old 'poor' level to mark lower part of donee: a:***************|***| 20(c) b:******************* 19 c:********************* 21 d:***********|*******| 20(b) e:*****|*************| 20(a) Finally, bring 'b' up to average with donor 'c'. Record donor and use old 'poor' level to mark lower part of donee: a:**************|****| 20(c) b:******************|| 20(c) c:*******************| 20 d:***********|*******| 20(b) e:*****|*************| 20(a) We now have a "squared histogram", i.e., a rectangle with 4 strips of equal area, each strip with two regions. A single uniform variate U can be used to generate a,b,c,d,e with the required probabilities, .32, .27, .26, .12 .06. Setup: Make tables, V[1]=a K[1]=c T[1]=0+16/20 V[2]=b K[2]=c T[2]=1+19/20 V[3]=c K[3]=c T[3]=2+20/20 V[4]=d K[4]=b T[4]=3+12/20 V[5]=e K[5]=a T[5]=4+ 6/20 Generation Process: Let j be the integer part of 1+5*U, with U uniform in (0,1). If U < T[j] return V[j], else return V[K[j]]. In many applications no V table is necessary: V[i]=i and the generating procedure becomes If U < T[j] return j, else return K[j]. For more, visit the Web site: Modeling & Simulation Resources. References & Further Readings: L'Ecuyer P., Uniform random number generation, Ann. Op. Res., 53, 77-120, 1994. L'Ecuyer P., Random number generation. In Handbook on Simulation, J. Banks (ed.), Wiley, 1998. Maurer U., A universal statistical test for random bit generators, J. Cryptology, 5, 89-105, 1992. Sobol' I., and Y. Levitan, A pseudo-random number generator for personal computers, Computers & Mathematics with Applications, 37(4), 33-40, 1999. Tsang W-W., A decision tree algorithm for squaring the histogram in random number generation, Ars Combinatoria, 23A, 291-301, 1987 We need to test for both randomness as well as uniformity. The tests can be classified in 2 categories: Empirical or statistical tests, and theoretical tests. Theoretical tests deal with the properties of the generator used to create the realization with desired distribution, and do not look at the number generated at all. For example, we would not use a generator with poor qualities to generate random numbers. Statistical tests are based solely on the random observations produced. Test for Randomness: A. Test for independence: B. Runs tests.(run-ups, run-downs): Test based on Normal approximation: C. Correlation tests: Frequency or Uniform Distribution Test: References & Further Readings: Modeling & SimulationSimulation in general is to pretend that one deals with a real thing while really working with an imitation. In operations research the imitation is a computer model of the simulated reality. A flight simulator on a PC is also a computer model of some aspects of the flight: it shows on the screen the controls and what the "pilot" (the youngster who operates it) is supposed to see from the "cockpit" (his armchair). Why to use models? To fly a simulator is safer and cheaper than the real airplane. For precisely this reason, models are used in industry commerce and military: it is very costly, dangerous and often impossible to make experiments with real systems. Provided that models are adequate descriptions of reality (they are valid), experimenting with them can save money, suffering and even time. When to use simulations? Systems that change with time, such as a gas station where cars come and go (called dynamic systems) and involve randomness. Nobody can guess at exactly which time the next car should arrive at the station, are good candidates for simulation. Modeling complex dynamic systems theoretically need too many simplifications and the emerging models may not be therefore valid. Simulation does not require that many simplifying assumptions, making it the only tool even in absence of randomness. How to simulate? Suppose we are interested in a gas station. We may describe the behavior of this system graphically by plotting the number of cars in the station; the state of the system. Every time a car arrives the graph increases by one unit while a departing car causes the graph to drop one unit. This graph (called sample path), could be obtained from observation of a real station, but could also be artificially constructed. Such artificial construction and the analysis of the resulting sample path (or more sample paths in more complex cases) consists of the simulation. Types of simulations: Discrete event. The above sample path consisted of only horizontal and vertical lines, as car arrivals and departures occurred at distinct points of time, what we refer to as events. Between two consecutive events, nothing happens - the graph is horizontal. When the number of events are finite, we call the simulation "discrete event." In some systems the state changes all the time, not just at the time of some discrete events. For example, the water level in a reservoir with given in and outflows may change all the time. In such cases "continuous simulation" is more appropriate, although discrete event simulation can serve as an approximation. Further consideration of discrete event simulations. How is simulation performed? Simulations may be performed manually. Most often, however, the system model is written either as a computer program (for an example click here) or as some kind of input into simulator software. System terminology: State: A variable characterizing an attribute in the system such as level of stock in inventory or number of jobs waiting for processing. Event: An occurrence at a point in time which may change the state of the system, such as arrival of a customer or start of work on a job. Entity: An object that passes through the system, such as cars in an intersection or orders in a factory. Often an event (e.g., arrival) is associated with an entity (e.g., customer). Queue: A queue is not only a physical queue of people, it can also be a task list, a buffer of finished goods waiting for transportation or any place where entities are waiting for something to happen for any reason. Creating: Creating is causing an arrival of a new entity to the system at some point in time. Scheduling: Scheduling is the act of assigning a new future event to an existing entity. Random variable: A random variable is a quantity that is uncertain, such as interarrival time between two incoming flights or number of defective parts in a shipment. Random variate: A random variate is an artificially generated random variable. Distribution: A distribution is the mathematical law which governs the probabilistic features of a random variable. A Simple Example: Building a simulation gas station with a single pump served by a single service man. Assume that arrival of cars as well their service times are random. At first identify the: states: number of cars waiting for service and number of cars served at any moment events: arrival of cars, start of service, end of service entities: these are the cars queue: the queue of cars in front of the pump, waiting for service random realizations: interarrival times, service times distributions: we shall assume exponential distributions for both the interarrival time and service time. Next, specify what to do at each event. The above example would look like this: At event of entity arrival: Create next arrival. If the server is free, send entity for start of service. Otherwise it joins the queue. At event of service start: Server becomes occupied. Schedule end of service for this entity. At event of service end: Server becomes free. If any entities waiting in queue: remove first entity from the queue; send it for start of service. Some initiation is still required, for example, the creation of the first arrival. Lastly, the above is translated into code. This is easy with an appropriate library which has subroutines for creation, scheduling, proper timing of events, queue manipulations, random variate generation and statistics collection. How to simulate? Besides the above, the program records the number of cars in the system before and after every change, together with the length of each event. A typical stochastic system has a large number of control parameters that can have a significant impact on the performance of the system. To establish a basic knowledge of the behavior of a system under variation of input parameter values and to estimate the relative importance of the input parameters, sensitivity analysis applies small changes to the nominal values of input parameters. For systems simulation, variations of the input parameter values cannot be made infinitely small. The sensitivity of the performance measure with respect to an input parameter is therefore defined as (partial) derivative. Sensitivity analysis is concerned with evaluating sensitivities (gradients, Hessian, etc.) of performance measures with respect to parameters of interest. It provides guidance for design and operational decisions and plays a pivotal role in identifying the most significant system parameters, as well as bottleneck subsystems. I have carried out research in the fields of sensitivity analysis and stochastic optimization of discrete event systems with an emphasis on computer simulation models. This part of lecture is dedicated to the estimation of an entire response surface of complex discrete event systems (DES) from a single sample path (simulation), such as the expected waiting time of a customer in a queuing network, with respect to the controllable parameters of the system, such as service rates, buffer sizes and routing probabilities. With the response surfaces at hand, we are able to perform sensitivity analysis and optimization of a DES from a single simulation, that is, to find the optimal parameters of the system and their sensitivities (derivatives), with respect to uncontrollable system parameters, such as arrival rates in a queuing network. We identified three distinct processes. Descriptive Analysis includes: Problem Identification & Formulation, Data Collection and Analysis, Computer Simulation Model Development, Validation, Verification and Calibration, and finally Performance Evaluation. Prescriptive Analysis: Optimization or Goal Seeking. These are necessary components for Post-prescriptive Analysis: Sensitivity, and What-If Analysis. The prescriptive simulation attempts to use simulation to prescribe decisions required to obtain specified results. It is subdivided into two topics- Goal Seeking and Optimization. Recent developments on "single-run" algorithms for the needed sensitivities (i.e. gradient, Hessian, etc.) make the prescriptive simulation feasible. Click on the image to enlarge it and THEN print it. Problem Formulation: Identify controllable and uncontrollable inputs. Identify constraints on the decision variables. Define measure of system performance and an objective function. Develop a preliminary model structure to interrelate the inputs and the measure of performance. Click on the image to enlarge it and THEN print it. Data Collection and Analysis: Regardless of the method used to collect the data, the decision of how much to collect is a trade-off between cost and accuracy. Simulation Model Development: Acquiring sufficient understanding of the system to develop an appropriate conceptual, logical and then simulation model is one of the most difficult tasks in simulation analysis. Model Validation, Verification and Calibration: In general, verification focuses on the internal consistency of a model, while validation is concerned with the correspondence between the model and the reality. The term validation is applied to those processes which seek to determine whether or not a simulation is correct with respect to the "real" system. More prosaically, validation is concerned with the question "Are we building the right system?". Verification, on the other hand, seeks to answer the question "Are we building the system right?" Verification checks that the implementation of the simulation model (program) corresponds to the model. Validation checks that the model corresponds to reality. Calibration checks that the data generated by the simulation matches real (observed) data. Validation: The process of comparing the model's output with the behavior of the phenomenon. In other words: comparing model execution to reality (physical or otherwise) Input and Output Analysis: Discrete-event simulation models typically have stochastic components that mimic the probabilistic nature of the system under consideration. Successful input modeling requires a close match between the input model and the true underlying probabilistic mechanism associated with the system. The input data analysis is to model an element (e.g., arrival process, service times) in a discrete-event simulation given a data set collected on the element of interest. This stage performs intensive error checking on the input data, including external, policy, random and deterministic variables. System simulation experiment is to learn about its behavior. Careful planning, or designing, of simulation experiments is generally a great help, saving time and effort by providing efficient ways to estimate the effects of changes in the model's inputs on its outputs. Statistical experimental-design methods are mostly used in the context of simulation experiments. Performance Evaluation and What-If Analysis: The `what-if' analysis is at the very heart of simulation models. Sensitivity Estimation: Users must be provided with affordable techniques for sensitivity analysis if they are to understand which relationships are meaningful in complicated models. Optimization: Traditional optimization techniques require gradient estimation. As with sensitivity analysis, the current approach for optimization requires intensive simulation to construct an approximate surface response function. Incorporating gradient estimation techniques into convergent algorithms such as Robbins-Monroe type algorithms for optimization purposes, will be considered. Gradient Estimation Applications: There are a number of applications which measure sensitivity information, (i.e., the gradient, Hessian, etc.), Local information, Structural properties, Response surface generation, Goal-seeking problem, Optimization, What-if Problem, and Meta-modelling Report Generating: Report generation is a critical link in the communication process between the model and the end user. A Classification of Stochastic ProcessesA stochastic process is a probabilistic model of a system that evolves randomly in time and space. Formally, a stochastic process is a collection of random variables {X(t), t Î T} all defined on a common sample (probability) space. The X(t) is the state while (time) t is the index that is a member of set T.Examples are the delay {D(i), i = 1, 2, ...} of the ith customer and number of customers {Q(t), T ³ 0} in the queue at time t in an M/M/1 queue. In the first example, we have a discrete- time, continuous state, while in the second example the state is discrete and time in continuous. The following table is a classification of various stochastic processes. The man made systems have mostly discrete state. Monte Carlo simulation deals with discrete time while in discrete even system simulation the time dimension is continuous, which is at the heart of this site.
Simulation Output Data and Stochastic ProcessesTo perform statistical analysis of the simulation output we need to establish some conditions, e.g. output data must be a covariance stationary process (e.g. the data collected over n simulation runs).Stationary Process (strictly stationary): A stationary stochastic process is a stochastic process {X(t), t Î T} with the property that the joint distribution all vectors of h dimension remain the same for any fixed h. First Order Stationary: A stochastic process is a first order stationary if expected of X(t) remains the same for all t. For example in economic time series, a process is first order stationary when we remove any kinds of trend by some mechanisms such as differencing. Second Order Stationary: A stochastic process is a second order stationary if it is first order stationary and covariance between X(t) and X(s) is function of t-s only. Again, in economic time series, a process is second order stationary when we stabilize also its variance by some kind of transformations such as taking square root. Clearly, a stationary process is a second order stationary, however the reverse may not hold. In simulation output statistical analysis we are satisfied if the output is covariance stationary. Covariance Stationary: A covariance stationary process is a stochastic process {X(t), t Î T} having finite second moments, i.e. expected of [X(t)]2 be finite. Clearly, any stationary process with finite second moment is covariance stationary. A stationary process may have no finite moment whatsoever. Since a Gaussian process needs a mean and covariance matrix only, it is stationary (strictly) if it is covariance stationary. Two Contrasting Stationary Process: Consider the following two extreme stochastic processes: - A sequence Y0, Y1,....., of independent identically distributed, random-value sequence is a stationary process, if its common distribution has a finite variance then the process is covariance stationary. - Let Z be a single random variable with known distribution function, and set Z0 = Z1 = ....Z. Note that in a realization of this process, the first element, Z0, may be random but after that there is no randomness. The process {Zi, i = 0, 1, 2, ..} is stationary if Z has a finite variance. Output data in simulation fall between these two type of process. Simulation outputs are identical, and mildly correlated (how mild? It depends on e.g. in a queueing system how large is the traffic intensity r). An example could be the delay process of the customers in a queueing system. Techniques for the Steady State SimulationUnlike in queuing theory where steady state results for some models are easily obtainable, the steady state simulation is not an easy task. The opposite is true for obtaining results for the transient period (i.e., the warm-up period).Gather steady state simulation output requires statistical assurance that the simulation model reached the steady state. The main difficulty is to obtain independent simulation runs with exclusion of the transient period . The two technique commonly used for steady state simulation are the Method of Batch means, and the Independent Replication. None of these two methods is superior to the other in all cases. Their performance depend on the magnitude of the traffic intensity. The other available technique is the Regenerative Method, which is mostly used for its theoretical nice properties, however it is rarely applied in actual simulation for obtaining the steady state output numerical results.
Suppose you have a regenerative simulation consisting of m cycles of size n1, n2,�nm, respectively. The cycle sums is: yi = S xij / ni, the sum is over j=1, 2, ..,ni The overall estimate is:Estimate = Syi / S ni, the sums are over i=1, 2, ..,m The 100(1-a/2)% confidence interval using the Z-table (or T-table, for m less than, say 30), is: Estimate ± Z. S/ (n. m�) where, n = S ni /m, the sum is over i=1, 2, ..,m and the variance is:S2 = S (yi - ni . Estimate)2/(m-1), the sum is over i=1, 2, ..,m Method of Batch Means: This method involves only one very long simulation run which is suitably subdivided into an initial transient period and n batches. Each of the batch is then treated as an independent run of the simulation experiment while no observation are made during the transient period which is treated as warm-up interval. Choosing a large batch interval size would effectively lead to independent batches and hence, independent runs of the simulation, however since number of batches are few on cannot invoke the central limit theorem to construct the needed confidence interval. On the other hand, choosing a small batch interval size would effectively lead to significant correlation between successive batches therefore cannot apply the results in constructing an accurate confidence interval.
Suppose you have n equal batches of m observations each. The means of each batch is: meani = S xij / m, the sum is over j=1, 2, ..,m The overall estimate is: Estimate = Smeani / n, the sum is over i=1, 2, ..,nThe 100(1-a/2)% confidence interval using the Z-table (or T-table, for n less than, say 30), is: Estimate ± Z. S where the variance is: S2 = S (meani - Estimate)2/(n-1), the sum is over i=1, 2, ..,nMethod of Independent Replications: This method is the most popularly used for systems with short transient period. This method requires independent runs of the simulation experiment different initial random seeds for the simulators' random number generator. For each independent replications of the simulation run it transient period is removed. For the observed intervals after the transient period data is collected and processed for the point estimates of the performance measure and for its subsequent confidence interval.
Suppose you have n replications with of m observations each. The means of each replication is: meani = S xij / m, the sum is over j=1, 2, ..,m The overall estimate is:Estimate = Smeani / n, the sum is over i=1, 2, ..,n The 100(1-a/2)% confidence interval using the Z-table (or T-table, for n less than, say 30), is: Estimate ± Z. S where the variance is: S2 = S (meani - Estimate)2/(n-1), the sum is over i=1, 2, ..,n Further Reading: Determination of the Warm-up PeriodTo estimate the long-term performance measure of the system, there are several methods such as Batch Means, Independent Replications and Regenerative Method.Batch Means is a method of estimating the steady-state characteristic from a single-run simulation. The single run is partitioned into equal size batches large enough for estimates obtained from different batches to be approximately independent. In the method of Batch Means, it is important to ensure that the bias due to initial conditions is removed to achieve at least a covariance stationary waiting time process. An obvious remedy is to run the simulation for a period large enough to remove the effect of the initial bias. During this warm-up period, no attempt is made to record the output of the simulation. The results are thrown away. At the end of this warm-up period, the waiting time of customers are collected for analysis. The practical question is "How long should the warm-up period be?". Abate and Whitt provided a relatively simple and nice expression for the time required (tp) for an M/M/1/ queue system (with traffic intensity r) starting at the origin (empty) to reach and remain within 100p% of the steady- state limit as follows: tp(r) = 2C(r) Ln {1/[(1-p)(1+2C(r))]}/(1-r)2 where C(r)=[2+ r + ( r2 + 4r )�] / 4. Some notions of tp(r) as a function of r and p, are given in following table:
Although this result is developed for M/M/1 queues, it has already been established that it can serve as an approximation for more general; i.e., GI/G/1 queues. Further Reading: Determination of the Desirable Number of Simulation RunsThe two widely used methods for experimentation on simulation models are method of bath means, and independent replications. Intuitively one may say the method of independent replication is superior in producing statistically a "good" estimate for the system's performance measure. In fact, not one method is superior in all cases and it all depends on the traffic intensity r.After deciding what method is more suitable to apply, the main question is determination of number of runs. That is, at the planning stage of a simulation investigation of the question of number of simulation runs (n) is critical. The confidence level of simulation output drawn from a set of simulation runs depends on the size of data set. The larger the number of runs, the higher is the associated confidence. However, more simulation runs also require more effort and resources for large systems. Thus, the main goal must be in finding the smallest number of simulation runs that will provide the desirable confidence. Pilot Studies: When the needed statistics for number of simulation runs calculation is not available from existing database, a pilot simulation is needed. For large pilot simulation runs (n), say over 30, the simplest number of runs determinate is: [(Za/2)2 S2] / d2 where d is the desirable margin of error (i.e., the absolute error), which is the half-length of the confidence interval with 100(1- a)% confidence interval. S2 is the variance obtained from the pilot run. One may use the following sample size determinate for a desirable relative error D in %, which requires an estimate of the coefficient of variation (C.V. in %) from a pilot run with n over 30:
These sample size determinates could also be used for simulation output estimation of unimodal output populations, with discrete or continuous random variables provided the pilot run size (n) is larger than (say) 30. The aim of applying any one of the above number of runs determinates is at improving your pilot estimates at feasible costs. You may like using the following Applet for determination of number of runs. Further Reading: At the planning stage of a simulation modeling the question of number of simulation runs (n) is critical. The following Java applets compute the needed Runs Size based on current avialable information ontained from a pilot simulation run, to achieve an acceptable accuracy and/or risk. Enter the needed information, and then click the Calculate button. The aim of applying any one of the following number of simulation runs determinates is at improving your pilot estimates at a feasible cost. Notes: The normality condition might be relaxed for number of simulation runs over, say 30. Moreover, determination of number of simulation runs for mean could also be used for other unimodal simulation output distributions including those with discrete random variables, such as proportion, provided the pilot run is sufficiently large (say, over 30). Runs' Size with Acceptable |