AN INTRODUCTION TO STATISTICAL ANALYSIS OF SIMULATION OUTPUTS

AN INTRODUCTION TO STATISTICAL ANALYSIS OF SIMULATION OUTPUTS

Sample Statistics n If x 1, x 2, …, xn are n observations of the value of an unknown quantity X, they constitute a sample of size n for the population on which X is defined. n Sample mean n Sample variance

Central Limit Theorem n If the n mutually independent random variables x 1, x 2, …, xn have the same distribution, and if their mean m and their variance s 2 exist then …

Central Limit Theorem n The random variable is distributed according to the standard normal distribution (zero mean and unit variance).

Estimating a mean n n Assume that we have a sample x 1, x 2, …, xn consisting of n independent * observations of a given population The sample mean xbar is an unbiased estimator of the mean m of the population * This is the critical assumption Without it, we cannot apply the formula

Large populations n n For large values of n, the (1 - )% confidence interval for m is given by with for the standard normal distribution

Explanations n n 1 – expressed in percent is the level of the confidence interval n 90% means = 0. 10 n 95% means = 0. 05 n 99% means = 0. 01 is the error probability n Probability that the true mean falls outside the confidence interval

95% Confidence Intervals n n For =. 05, za/2 =1. 96 Example: = 35, s = 4 and n = 100 n If n The 95% confidence interval for m is 35 ± 1. 96 x 4/ 10 = 35 ± 7. 84

When s is not known n We can replace s in the preceding formula by the standard-deviation s of the sample n When n < 30, we must read the value of za/2 from a table of Student's t-distribution with n - 1 degrees of freedom n When n 30, we can use the standard values

Confidence Intervals in CSIM can automatically compute confidence intervals for the mean value of any table, qtable and so on. n For everything but boxes n xyx->confidence(); n For the elapsed times in a box n bd->time_confidence(); n For the population of a box n bd->_number_confidence();

Confidence Intervals in CSIM n We get 90, 95 and 98 percent confidence intervals n Computed using batch means method See next section But only for the mean n n

Estimating a proportion n A proportion p represents the probability P(X ) for some fixed threshold n 97% of our customers have to wait less than one minute n Confidence intervals for proportions are much easier to compute than confidence intervals for quantiles

Assumptions n n n Assume that we have n independent observations x 1, x 2, …, xn of a given population variable X and that this variable has a continuous distribution Let p represent the proportion we want to estimate, say P(X ) Let k represent the number of observations that are .

Basic property n If p represent the proportion we want to estimate, say P(X ), and k represent the number of observations that are n The rv k is distributed according to a binomial distribution with mean np and variance p(1 – p) n The rv k/n is distributed according to a binomial distribution with mean p and variance p(1 – p)/n

The formula n When n > 29, we can use the Wilson’s interval where za/2 = 1. 96 for a 95% C. I.

Advantages n n Produces tighter confidence intervals then the Central Limit Theorem Works when is equal to zero

Example n Assume we have 400, 000 independent observations without noting any failure n The 95% confidence interval for the probability that the system could fail is

AUTOCORRELATION

The problem n n n All statistical analysis techniques we have discussed assume that sample values are mutually independent Generally false for quantities such as n Waiting times, response times, … Tend instead to be autocorrelated n When the waiting lines are long, everybody wait a long time

Traditional solution n n Keep the measurements sufficiently apart n Sample them every T minutes apart n Standard solution for collecting observations on a running system Not practical n Would require very long simulations

Three good solutions n n n Batch means Regenerative method Time series analysis

Batch means n n n We group consecutive observations into “batches” We compute the means of these batches We observe that autocorrelation among batch means decreases with size of batches n When size increases, each batch includes more observations that are far apart

Example n n n We collected the following values: n 4, 3, 3, 4, 5, 5, 3, 2, 2, 3 We group them into two batches of five observations: n 4, 3, 3, 4, 5 and 5, 3, 2, 2, 3 The batch means are: n 3. 8 and 3

Batch means in CSIM uses fixed-size batches n To compute confidence intervals n To control the duration of a simulation (run-length control)

Regenerative method n n Most systems with queues go through states that return it to an state identical to its original state n The system regenerates itself Examples: n Whenever a disk array is brought back to its original state n Whenever a camper rental agency has all its campers available

Key idea n Define a regeneration interval as an interval between two consecutive regeneration points: n Observations collected during the same regeneration interval can be correlated n Observations collected during different regeneration intervals are independent

Application n We group together observations that occur within the same regeneration interval We compute the means of these groups of observations These group means are independent from each other

Limitations of the approach n Not general: n System must go through regeneration points n n System must be idle Leads to complex computations n We rarely have exactly the same number of observations in two different regenerations intervals

Time series analysis n n n Treats consecutive observations as elements of a time series Estimates autocorrelation among the elements of a time series Includes this autocorrelation in the computation of all confidence intervals

RUN LENGTH CONTROL

Objective n n Accuracy of confidence intervals increases with duration of simulation n The 1/ n factor We would like to be able to stop the simulation once a given accuracy level has been reached for the confidence interval of a specific measurement

The CSIM solution (I) n Specify n Quantity of interest n Mean value from a table, a qtable, … n A relative accuracy n Maximum relative error n A confidence level (say, 0. 90 to 0. 99) n A maximum simulation duration n In seconds of CPU time

The CSIM solution (II) n n n void table: : run_length(double accuracy, double conf_level, double max_time) void qtable: : run_length(double accuracy, double conf_level, double max_time) void meter: : run_length(double accuracy, double conf_level, double max_time) void box: : time_run_length(double accuracy, double conf_level, double max_time) void box: : number_run_length(double accuracy, double conf_level, double max_time)

The CSIM solution (III) n Example: n thistable->run_length(. 01, . 95, 500) n Specifies than we want an error less than 1 percent for 95% confidence interval of mean of thistable n Stops simulation after 500 seconds n thisbox->: time_run_length(. 01, . 99, 500) n Same for time average of a box

The CSIM solution (IV) n Replace termination test by n converged. wait();

Example (I) n n n The campers We want n A maximum error of 1 percent (0. 01) n For the 95 percent confidence interval n Of average number of rented campers n A maximum simulation time of 100 s Will use number aspect of agency box

Example (II) n Add n agency->number_confidence() after activation of box agency in csim process

Example (III) n n Put main loop n while(simtime() < DURATION) { hold(exponential(MIART); customer(); } in a separate arrivals process Make it an infinite loop

Example (IV) n n Add n converged. wait(); after call to arrivals process n arrivals(): converged. wait(); Best way to let sim process generate customers and wait for termination in parallel

Example (V) n The new arrivals process void arrivals() { process(“arrivals”); // REQUIRED for(; ; ) { // forever loop hold(exponential(MIART)); customer(); } // forever } // arrivals

Warnings n Confidence intervals do not take into account model inaccuracies n n While the batch means method eliminates most effects of measurement autocorrelation, it is not always 100% effective The max_time parameter of the run_length() will not necessarily stop the simulation just after the specified CPU time n Like the emergency brake of a train

Objective n n Partition processes into different classes n Low priority n High priority Obtain separate statistics for each process priority class

Declaring and initializing process classes n n To declare a dynamic process class: n process_class *c; To initialize a process class before it can be used in any other statement. n c = new process_class("low priority")

To assign a priority class n n Add inside the process n c->set_process_class(); Processes that have not been assigned a process class belong to the “default” process class

Reporting n Can use n report() n report_classes()

Other options n n Can change the name of a process class: n c->set_name("high priority"); Can reset statistics associated with a process n c ->reset(); Can do the same for all process classes: n reset_process_classes(); Can delete a dynamic process class: n delete c;

RANDOM NUMBERS

Background n We need to distinguish between n Truly random numbers n Obtained through observations of a physical random process n Rolling dices, roulette n Atomic decay n Pseudo-random numbers n Obtained through arithmetic operations

Example n Linear Congruential Generators (LCG) n Easy to implement and fast n Defined by the recurrence relation: rn+1 = (a rn + c ) mod m where r 1, r 2, … are the random values n m is the "modulus“ n r 0 is the seed n

Two realizations n n GCC family of compilers 32 n m = 2 a = 69069 c=5 Microsoft Visual/Quick C/C++ 32 n m = 2 a = 214013 c = 2531011

Problems with pseudorandom number generators n n Much shorter periods for some seed states Lack of uniformity of distribution Correlation of successive values …

Better RNGs n n Use the Mersenne twister 19937 - 1 n Period is 2 Blum-Schub

A quote n "Any one who considers arithmetical methods of producing random digits is, of course, in a state of sin. For, as has been pointed out several times, there is no such thing as a random number– there are only methods to produce random numbers, and a strict arithmetic procedure of course is not such a method. “ John von Neumann

CSIM RNGs n n BY default, CSIM uses a single stream of random numbers Can reset the seed using n void reseed(stream *s, long n) as in n reseed(NIL, 13579)

Continuous distributions supported by CSIM (I) n n n double uniform(double min, double max) double triangular(double min, double max, double mode) double beta(double min, double max, double shape 1, double shape 2) double exponential(double mean) double gamma(double mean, double stddev) double erlang(double mean, double var)

Continuous distributions supported by CSIM (II) n n n double hyperx(double mean, double var) double weibull(double shape, double scale) double normal(double mean, double stddev) double lognormal(double mean, double stddev) double cauchy(double alpha, double beta) double hypoexponential(double mn, double var)

Continuous distributions supported by CSIM (III) n n n double pareto(double a) double zipf(long n) double zipf_sum(long n, double *sum)

Discrete distributions supported by CSIM n n n long uniform_int(long min, long max) long bernoulli(double prob_success) long binomial(double prob_success, long num_trials) long geometric(double prob_success) long negative_binomial(long success_num, double prob_success) long poisson(double mean)

Empirical distributions supported by CSIM n ? ? ?

Using multiple streams n n In campers example, sequence of RNs used to generate arrivals is affected by the numbers of campers n If agency has less campers n More customers will be lost n Lost customers do no generate any RN Better to have separate random number streams for arrivals and service times

Declaring and initialing random number streams n n Can use: n stream *s; s = new stream(); By default, streams are created with seeds that are spaced 100, 000 values apart

Reseeding a stream n Use: n s->reseed(24680); where the new seed is a long integer

Other functions n n Inspect the current state of a stream n i = s->state(); Delete a stream n delete s;

Using a specific seed n Prefix RNG function with name of seed: n s->uniform (3. 0, 7. 0)

A CASE STUDY

RAID array revisited n n Reliability of a RAID array Reliability R(t) of a system is the probability that will remain operational over a time interval [0. t ] given that it was operational at time t = 0 n Not the same as availability n Our focus is evaluating the risk of a data loss during array lifetime

Executing multiple runs n n Multiple runs provide statistically independent repetitions of original simulation n Useful for n Collecting more accurate results n Constructing confidence intervals Use rerun() function within a loop n create("sim") call must be inside that loop

Overview of sim function extern "C" void sim() { // sim process runcount = 0; while(runcount < NRUNS) { create("sim"); // make it a process // usual contents of sim process rerun(); runcount++; } // while report_hdr(); // produce statistics report } // sim

Global includes and defines #include <cpp. h> // CSIM C++ header file #define NDISKS 5 // number of disks in array #define NYEARS 5 // useful lifetime of array #define MTTF 300000. 0 // disk MTTF #define MTTR 24. 0 // disk MTTR #define NRUNS 100000 // no of runs void disk(int i);

Global declarations int nfailed; // number of failed disks int runcount = 0; // number of runs int ndatalosses = 0; //counter double lifetime = NYEARS * 365 * 24; // simulation duration

Sim function (I) extern "C" void sim() { // sim process int i; runcount = 0; while(runcount < NRUNS) { create("sim"); // make this a process dataloss = new event("dataloss"); dataloss->clear(); nfailed = 0;

Sim function (II) // create NDISKS disk processes for (i=0; i < NDISKS; i++){ disk(i); } // for dataloss->timed_wait(lifetime); if (simtime() < lifetime) { ndatalosses++; } // if rerun(); runcount++; } // while

Sim function (III) report_hdr(); // produce statistics report printf("Array lifetime %fd years", NYEARS); printf(“or %f hoursn", lifetime); printf("Sim time %fn", simtime()); printf("Disk MTTF %f hoursn", MTTF); printf("Disk MTTR %ff hoursn", MTTR); printf(“Completed runs %dn", runcount); printf(“Data lossese%dn", ndatalosses); } // sim

Disk process (I) void disk(int i) { create("disk"); while(simtime() < lifetime || dataloss) { hold(exponential(MTTF)); // disk failed nfailed++;

Disk process (II) if (nfailed == 2) { dataloss = 1; failtime = simtime(); finish->set(); terminate(); } // if hold(MTTR); // repair process nfailed--; // disk is replaced } // while } // disk

Simulation Outcome C++/CSIM Simulation Report (Version 19. 0 for Linux x 86) Tue Apr 22 10: 13: 14 2008 Ending simulation time: Elapsed simulation time: CPU time used (seconds): Array lifetime 5 years which corresponds to 43800 hours Simulated time 0 Disk MTTF 300000 hours Disk MTTR 24 hours Number of runs completed 100000 Number in runs ending in data loss 17 0. 000

Discussion n Whole program took less than 98 seconds on linux 03 server Simulated time should be equal to 43, 800 hours, not zero. n rerun() artifact? We observe 17 failures out of 100, 000 runs n Data survival rate is 99. 983 percent n Three nines

Confidence intervals n n n Data loss rate and survival rate are distributed according to a binomial law Since n = 100, 000, the distributions of both proportions are approximately normal Will use

CI for data survival rate s n n ŝ =99. 983% or 0. 99983 95% CI is 0. 999849 ± 0. 000083 [0. 999766, 0. 999932] between three and four nines

What are nines? n n n A measure of system reliability/availability n 99% is two nines n 99. 9% is three nines More formally, we compute -log 10 (1. 0 – x) Our confidence interval could then be expressed as [3. 63 nines, 4. 17 nines]

CI for data loss rate 0. 0000189 ± 0. 000088 [0. 0000106, 0. 000273]

Analytical approach n Will use a Markov model: n Works since failure rates and repair rates are distributed according to exponential laws n Obviously not true for repair rates

The model n n State label indicates number of failed drives Failure state is an absorbing state

Using ergodic hypothesis n n Assume that array is returned to original state after each data loss Failure rate will be 4 lp 1

Model equations n n 5λp 0=(μ+ 4λ)p 1 po + p 1 = 1 n Solution is n Failure rate is

Computing five year survival n n Decay is For MTTF = 300, 000 hours, MTTR = 24 hours and t = 5 years n Data survival rate is 0. 999767 or 3. 632 nines n Barely inside the C. I.

Why? n n n Markov model assumed repair times were exponentially distributed n Required by the technique Simulation model assumed deterministic repair times n Somewhat more realistic Difference illustrates a major limitation of stochastic approach

Checking the explanation (I) n We repeat the simulation using repair times that are exponentially distributed n We observe 19 data losses out of 100, 000 runs n qhat is 0. 99981 n Confidence interval is (3. 588 nines, 4. 080 nines)

Checking the explanation (II) n n We repeat a second time the simulation using repair times that are exponentially distributed and collecting more data n We observe 95 data losses out of 400, 000 runs n qhat is 0. 999763 n Confidence interval is (3. 552 nines, 3. 734 nines) Markov model predicted 3. 632 nines

Comparing both approaches n Which is the best approach in terms of n Generality and flexibility? n n Provided results? Discrete simulation and stochastic modeling complement each other

Extensions (I) n Other repair time distributions n Exponential ( to compare with results of stochastic analysis) n Ad hoc (80% of repairs within one day, 20% within two days) n Taking accounts of differences between day and night, weekdays and weekends n Will map simtime() into a calendar

Extensions (II) n Other disk arrays n Mirrored disks: NDISKS = 2 n RAID level 6: tolerate two disk failures n Triplicate disks: tolerate failure of two out of three disks n SSPi. RAL arrays n See next slide n Will require more work

SSPi. RAL arrays ABC BCD B A DAB D C CDA Resists all triple failures and most quadruple failures

CONCLUSION

CONCLUSION n n n Statistical analysis of outputs is an important aspect of simulation Standard statistical techniques require independent observations w/o autocorrelation n Simplest solution is to use batch means CSIM tools simplify constructions of confidence intervals for all sorts of means