Output Data Analysis for a Single System 9.1 Introduction................................................................................................2

CHAPTER 9

Output Data Analysis for a Single System

9.1 Introduction ...2

9.2 Transient and Steady-State Behavior of a Stochastic Process ...10

9.3 Types of Simulations with Regard to Output Analysis ...12

9.4 Statistical Analysis for Terminating Simulations...14

9.4.1 Estimating Means ...16

9.4.2 Estimating Other Measures of Performance...22

9.4.3 Choosing Initial Conditions ...24

9.5 Statistical Analysis for Steady-State Parameters ...25

9.5.1 The Problem of the Initial Transient...26

9.5.2 Replication/Deletion Approaches for Means ...27

9.5.3 Other Approaches for Means ...28

9.5.4 Estimating Other Measures of Performance...40

9.6 Statistical Analysis for Steady-State Cycle Parameters...41

9.7 Multiple Measures of Performance ...43

9.8 Time Plots of Important Variables...45

9.1 Introduction

Basic, most serious disadvantage of simulation:

With a stochastic simulation, don’t get exact “answers” from run Two different runs of same model ⇒ different numerical results Random input:

Random numbers Random variates

→

model/code

→

Performance measures Thus, the output performance measures:

Really observations from their probability distribution Ask questions about this distribution:

Mean (expected value): E(average WIP) Variance: Var(average WIP) Probabilities: P(average WIP > 250)

Quantiles: What value of x is such that P(avg. WIP > x) ≤ 0.02?

Interpreting simulation output: statistical analysis of output data

Failure to recognize, deal with randomness in simulation output can lead to serious errors, misinterpretation, bad decisions

Also, must take care to use appropriate statistical methods, since simulation output data are usually nonstationary, autocorrelated, and non-normal, contrary to assumptions behind classical IID statistical methods

Enhanced computer power and speed is making it much easier to carry out

The Statistical Nature of Simulation Output

Let Y₁, Y₂, ... be an output process from a single simulation run Yi = delay in queue of ith arriving customer

Y_i = production in ith hour in a factory

Y_i’s are random variables that are generally neither independent nor identically distributed (nor normally distributed), so classical IID normal-theory statistical methods don’t apply directly to the Y_i’s

Let y11, y12, ..., y1m be a realization of the random variables Y1, Y2, ..., Ym resulting from making a single simulation run of length m observations, using a particular stream of underlying U(0, 1) random numbers.

If we use a separate stream of random numbers for another simulation run of this same length, we get a realization y21, y22, ..., y2m that is independent of, but identically distributed to, y₁₁, y₁₂, ..., y_1m

Make n such independent runs, each using “fresh” random numbers, to get y11, y12, ..., y1m

y₂₁, y₂₂, ..., y_2m

. . .

y_n1, y_n2, ..., y_nm

Within a row: not IID Across the ith column: IID

realizations of the r.v. Y_i (but still not necessarly normally distributed)

Can compute a summary measure within a run, and then the summary measures across the runs are IID (but still not necessarily normally distributed)

Bank with 5 tellers, one FIFO queue, open 9am-5pm, flush out before stopping Interarrivals ~ expo (mean = 1 min.), service times ~ expo (mean = 4 min.) Summary measures from 10 runs (replications):

Clearly, there’s variation across runs;

need appropriate statistical analysis

Types of Output Performance Measures What do you want to know about the system?

Average time in system

Worst (longest) time in system Average, worst time in queue(s)

Average, worst, best number of “good” pieces produced per day

Variability (standard deviation, range) of number of parts produced per day Average, maximum number of parts in the system (WIP)

Average, maximum length of queue(s)

Proportion of time a machine is down, up and busy, up and idle Ask the same questions of the model/code

Think ahead: Asking for additional output performance measures can change how the simulation code is written, or even how the system is modeled

Simple queueing model:

Server

Customers in queue

Customer in service Customer

arrivals Customer

departures

Want: Average number of customers in queue Proportion of time server is busy

Maybe: Average time customers spend in queue

Question: How does wanting the average time in queue affect how the model is coded, data structures, etc.?

As simulation progresses through time, there are typically two kinds of processes that are observed:

Discrete-time process: There is a natural “first” observation, “second”

observation, etc.—but can only observe them when they “happen”

S_i = time in system of ith part produced, i ∈ {1, 2, ...}

Suppose there are M parts produced during the simulation

i

1 2 3 ... M

Si

Typical discrete-time output performance measures:

Average time in system:

M S M

S

M

i

∑

= =¹

) (

Maximum time in system: _i

M

i S

M S

,..., 2 , 1

max )

(

* = =

Proportion of parts that were in system more than 60 minutes:

M S I

M P

M

i

∑

= ∞

= ^{1 (60, )}

60

) ( )

( , where





≤

= >

∞ 0 if 60

60 if

) 1

)(

, 60 (

i i

i S

S S I

Other examples of discrete-time processes:

D_i = delay of ith customer in queue

Y_i = throughput (production) during ith hour B_i = 1 if caller i gets a busy signal, 0 otherwise

Continuous-time process: Can jump into system at any point in time (real, continuous time) and take a “snapshot” of something — there is no natural

“first” or “second” observation

Q(t) = number of parts in a particular queue at time t ∈ [0, ∞) Run simulation for T units of simulated time

Q(t )

0 1 2 3

t T

Typical continuous-time output performance measures:

Time-average length of queue:

T dt t Q T

Q

∫

= ⁰ ( ) )

(

Maximum length of queue: *( ) max ( )

0 Q t

T Q

T t≤

= ≤

Proportion of time that there were more than two in the queue:

T

dt t Q I

T P

∫

= ^{0 (2, )}

2

)) ( ( )

(

Another important kind of continuous-time statistic: utilizations

Let =

t t t

B 0 if server isidle at time at time busy

is server if

) 1 (

0 1

t T

B(t )

Server utilization (proportion of time busy):

T dt t B T

U

∫

= ⁰ ( ) )

(

Other examples of continuous-time processes:

N(t) = number of parts in shop at time t (WIP)

D(t) = number of machines down at time t

Typically, we want to observe several (maybe lots of) different performance measures from the same system/model

Usually low additional cost/hassle to do so, can always ignore later But not getting a particular output measure could imply rerunning

Difficulty in statistical analysis of output with several performance measures:

May want to make several simultaneous estimates, statements Be careful how this is done, what is said

Multiple-comparisons problem in statistics literature (more in Sec. 9.7)

9.2 Transient and Steady-State Behavior of a Stochastic Process

Output process (discrete-time) Y1, Y2, ...

Let F_i(y | I) = P(Y_i≤ y | I) be the transient (cumulative) distribution of the process at (discrete) time i

In general, F_i depends on both the time i and the initial condition I Corresponding transient density functions:

If there is a CDF F(y) such that F_i(y | I) → F(y) as i → ∞ for all y and for all initial conditions I, then F(y) is the steady-state distribution of the output process F(y) may or may not exist

F(y) must be independent of the initial conditions — same for all I

Roughly speaking, if there is a time index k such that for i > k F_i(y | I) ≈ F(y) in some sense, then we say that the process is “in steady state” after time k

Even though the distribution of the Yi’s after time k is not appreciably changing, observations on the Y_i’s could still have large variance and thus “bounce around” a lot — they’re just not systematically trending any more

Steady-state distribution does not depend on initial conditions, but the nature and rate of convergence of the transient distributions can depend heavily on the initial conditions

M/M/1 queue, E(delay in queue), different number of customers s present initially:

Inventory system, E(cost in month i):

9.3 Types of Simulations with Regard to Output Analysis

Terminating: Parameters to be estimated are defined relative to specific initial and stopping conditions that are part of the model

There is a “natural” and realistic way to model both the initial and stopping conditions

Output performance measures generally depend on both the initial and stopping conditions

Nonterminating: There is no natural and realistic event that terminates the model

Interested in “long-run” behavior characteristic of “normal” operation If the performance measure of interest is a characteristic of a steady-state

distribution of the process, it is a steady-state parameter of the model Theoretically, does not depend on initial conditions

Practically, must ensure that run is long enough so that initial-condition effects have dissipated

Not all nonterminating systems are steady-state: there could be a periodic

“cycle” in the long run, giving rise to steady-state cycle parameters

Focus on terminating vs. steady-state simulations:

Which is appropriate?

Depends on goals of the study

Statistical analysis for terminating simulations is a lot easier Is steady-state relevant at all? Maybe:

24 hours/day, “lights-out” manufacturing Global telecommunications

Design conservatively for peak load of infinite duration Some examples:

Physical model

Terminating estimand Steady-state estimand Single-server

queue

Expected average delay in queue of first 25 customers, given empty-and-idle initial conditions

Long-run expected delay in queue of a customer

Manufacturing system

Expected daily production, given some number of

workpieces in process initially

Expected long-run daily production

Reliability system

Expected lifetime, or

probability that it lasts at least a year, given all components initially new, working

Probably not sensible

Battlefield model

Probability that attacking force loses half its strength before defending force loses half its strength

Probably not sensible

9.4 Statistical Analysis for Terminating Simulations

Make n IID replications of a terminating simulation Same initial conditions for each replication Same terminating event for each replication Separate random numbers for each replication

Let X_j be a summary measure of interest from the jth replication

e.g., Xj = the average delay in queue of all customers in the jth replication

Then X₁, X₂, ..., X_n are IID random variables, can apply classical statistical analysis to them

Rely on central limit theorem to justify normality assumption even though it’s generally not true

So basic strategy is replication of the whole simulation some number n of times One simulation run is a sample of size one (not worth much statistically)

What About Classical Statistics?

Classical statistical methods don’t work directly within a simulation run, due to autocorrelation usually present

Example: Delays in a queue of individual jobs: D₁, D₂, D₃, ..., D_m Want to estimate µ = E(average delay of the m jobs)

Sample mean D m D m

m

i

∑

=

=

1

)

( is an unbiased estimator for µ

Need to estimate Var(D(m)) for confidence intervals on µ, test hypotheses like H0: µ = µ₀

But “sample variance”

(

) [

]

1

2 −

∑

=

m m m

D D

m

i

i may be severely biased for

Var(D(m)) Reason:

Corr(D_i, D_{i + l}) ≠ 0, in general

Unbiasedness of variance estimators follows from independence of data, which is not true within a simulation

Usual situation:

Positive autocorrelation: Corr(D_i, D_{i + l}) > 0

Causes

( ) [ ]









∑

=

) 1 ( )

(

1

2 m m

m D D E

m

i

i < Var(D(m)) — maybe far too small Intuition:

D_{i + 1} is pretty much the same as D_i

D_i’s are more stable than if they were independent Their sample variance is understated

Thus, must take care to estimate variances properly: understating variances Have too much faith in our point estimates

Believe our simulation results too much

9.4.1 Estimating Means

Want: Estimate of some parameter µ of the process Often (not always): µ = E(something)

µ = E(average delay in queue of m customers) µ = E(time-average WIP)

µ = E(machine utilization)

µ = P(average delay > 5 minutes) = E[I_(5,∞)(average delay)]

Point estimate: µˆ = 12.3

How close is µˆ to the true unknown value of µ?

Customary, useful method for assessing precision of estimator: confidence interval for µ

Pick confidence level 1 – α (typically 0.90, 0.95, etc.)

Use simulation output to construct an interval [A, B] that covers µ with probability 1 – α

Interpretation: 100(1 – α)% of the time, the interval formed in this way will cover µ

µ (unknown)

Wrong interpretation: “I’m 95% sure that µ is between 9.4 and 11.1”

Common approach to statistical analysis of simulation output:

Can’t do “classical” (IID, unbiased) statistics within a simulation run Try to modify setup, design, to get back to classical statistics

In terminating simulations, this is conceptually easy:

Make n independent replications of the whole simulation Let X_j be the performance measure from the jth replication

X_j = average of the delays in queue X_j = time-average WIP

X_j = utilization of a bottleneck machine

Then X₁, X₂, ..., X_n are IID and unbiased for µ = E(Xj)

Apply classical statistics to X_j’s, not to observations within a run Approximate 100(1 – α)% confidence interval for µ:

n X n

X

n

j

∑

= =¹

)

( is an unbiased estimator of µ

1 )) ( (

) (

2 2

−

−

=

∑

n

n X X n

S

j

is an unbiased estimator of Var(X_j)

n n t S

n

X _n ( )

)

( ± _{−1,1−α/2} covers µ with approximate probability 1 – α

(t_{n–1,1–α/2} = point below which is area (probability) 1–α/2 in Student’s t distribution with n – 1 d.f.)

Most important point:

The “basic ingredients” to the statistical analysis are the performance measures from the different, independent replications

One whole simulation run = a “sample” of size one (not worth much)

Example: n = 10 replications of single-server queue X_j = average delay in queue from jth replication

X_j’s: 2.02, 0.73, 3.20, 6.23, 1.76, 0.47, 3.89, 5.45, 1.44, 1.23 Want 90% confidence interval, i.e., α = 0.10

) 10 (

X = 2.64, S²(10) = 3.96, t_{9, 0.95} = 1.833

Approximate 90% confidence interval is 2.64 ± 1.15, or [1.49, 3.79]

Other examples in text:

Inventory model

Estimation of proportions

Why “approximate” 90% confidence interval?

Assumes X_j’s are normally distributed — never true, but does it matter?

Central-limit theorem:

As n (number of replications) grows, coverage probability → 1 – α Robustness study with this model:

Know µ = 2.12 from queueing theory Pick an n (like n = 10)

Make 500 “90%” confidence intervals (total no. runs = 10 × 500) Count % of these 500 intervals covering µ = 2.12:

n Estimated coverage (want 90%) 5

10 20 40

88%

86%

89%

92%

Bad news: actual coverages can depend (a lot) on the model Reliability model: ¹

2

3

Components fail independently

Times T_i to failure in components ~ Weibull(α = 0.5, β = 1.0) Time to system failure = T = min(T₁, max(T₂, T₃))

µ = E(time to system failure) = 0.78

n Estimated coverage (want 90%) 5

10 20 40

71%

75%

80%

84%

Obtaining a Specified Precision

If the number n of replication is chosen too small, the confidence intervals might too wide to be useful

M/M/1 example:

90% confidence interval from n = 10 replications: 2.64 ± 1.15, or [1.49, 3.79]

Half width (1.15) is 44% of point estimate (2.64) Equivalently: 2.64 ± 44%, not very precise

Half-width expression: δ(α, n) =

n n t_n S( )

2 / 1 , 1 −α

−

To decrease: α ↑: undesirable, since α = probability of missing S(n) ↓: estimates Var(X _j), which is fixed (maybe ...) n ↑: more replications

Sequential sampling: Increase n until δ(α, n) is “small enough”

Two versions of what “small enough” might mean (more details in text):

Absolute precision:

Specify β > 0, want n big enough so that δ(α, n) < β

Requires at least some knowledge of context to set meaningful β Relative precision:

Specify γ (0 < γ < 1), want n big enough so that δ(α,n) X(n)<γ

Need not know much about context to set meaningful γ Notes:

Above description leaves out a few technical details; see text

“Fixes” robustness issue: As β or γ → 0, coverage probability → 1 – α

Can be dangerous for small β or γ : Required n increases quadratically as β or γ decrease

May be difficult to automate with a simulation language, depending on modeling constructs available, what automated statistical capabilities present, and what access the user has to internal software variables

9.4.2 Estimating Other Measures of Performance

Sometimes can miss important system/model characteristics if we look only at averages

Other measures: Proportions, variances, quantiles

Proportions

Compare two operating policies for queueing system with five servers

vs.

Estimates

Performance measure Five queues Single queue Average delay in queue

Average number in queue(s) Number of delays ≥ 20 minutes

5.57 minutes 5.52

33

5.57 minutes 5.52

6

Variances (or Standard Deviations)

Interested in process variability X_j = daily production of good items Want to estimate Var(X _j)

Make n replications, compute S(n) as before

Confidence interval on Var(X _j): use chi-square distribution

Quantiles

Inventory system, X_j = maximum inventory level during the horizon Want to determine storage capacity that is sufficient with probability 0.98 Want to find x such that P(X_j ≤ x) = 0.98

One approach (more details in text):

Make n replications, observe X₁, X₂, ..., X_n Sort the X_j’s into increasing order

Estimate x to be a value below which are 98% of the X_j’s

9.4.3 Choosing Initial Conditions

For terminating simulations, the initial conditions can affect the output performance measure, so the simulations should be initialized appropriately

Example: Want to estimate expected average delay in queue of bank customers who arrive and complete their delay between noon and 1:00pm

Bank is likely to be crowded already at noon, so starting empty and idle at noon will probably bias the results low

Two possible remedies:

If bank actually opens at 9:00am, start the simulation empty and idle, let it run for 3 simulated hours, clear the statistical accumulators, and observe statistics for the next simulated hour

Take data in the field on number of customers present at noon, fit a (discrete) distribution to it, and draw from this distribution to initialize the simulation at time 0 = noon. Draw independently from this distribution to initialize multiple replications.

Note: This could be difficult in simulation software, depending on the modeling constructs available

9.5 Statistical Analysis for Steady-State Parameters

Much more difficult problem than analysis for terminating simulations Want to estimate (e.g. for discrete- and continuous-time processes)







=

= =

∞

→

∞

→

t t

Q t Q E

i D

D E

t

i i i

at time queue

in number )

( )) ( ( lim

customer th

of queue in

delay )

( ν lim

Basic question for designing runs:

Many short

runs One long

run X₁

X₂

X3

X₄

X₅

X₁

Many short runs One long run

Good Simple (same as terminating) Get IID data

Less point-estimator bias No restarts

9.5.1 The Problem of the Initial Transient

If steady-state is the goal, initial conditions will generally bias the results of the simulation for some initial period of time

Most common technique is to warm up the model, also called initial-data deletion Identify index l (for discrete-time processes) or time t₀ (for continuous-time

processes) beyond which the output appears not to be drifting any more Clear statistical accumulators at that time

Start over with data collection, and “count” only data observed past that point After warmup period, observations will still have variance, so will bounce around

— they are just not “trending” any more

Facility for doing this in most simulation software (but the user must specify the warmup period)

Challenge — identifying a warmup period that is long enough, yet no so long as to be excessively wasteful of data — see text for details and examples

Some statistical-analysis tests have been devised

Most practical (and widely-used) method is to make plots, perhaps averaged across and within replications to dampen the noise, and “eyeball” a cutoff If there are multiple output processes, and if they disagree about what the

appropriate warmup period is, a decision must be made whether to use

different warmups for each process, or to use the same one for each process

— which would have to be the maximum of the individual warmups, to be save, and so would be conservative for most of the output processes

A different approach: Try to find “smarter” initialization states or distributions that

9.5.2 Replication/Deletion Approaches for Means

Assume that an appropriate warmup period has been determined

X_j = output measure on jth replication, collected only past warmup point Proceed with statistical analysis exactly as in terminating case

Make independent replications, each warmed up

Compute mean, variance estimates across replications, confidence intervals Advantages (compared to methods to be discussed below):

Simple — aside from warmup, the same as for terminating simulations

Get truly IID observations — important not only for variance estimation and confidence-interval construction, but also for more sophisticated statistical goals and techniques to be discussed in Chap. 10–12

Disadvantages:

No completely reliable method to identify an appropriate warmup Too long ⇒ wasteful of data

Too short ⇒ point-estimator bias, which can have serious consequences, especially if used in concert with a sequential procedure:

) (n

X is no longer an unbiased estimator of ν

Confidence interval is centered in the –wrong” place— at E[X(n)] ≠ ν As n ↑, confidence interval shrinks down around the wrong point, causing

coverage to drop

n=10 n=20 n=30 n=40

Work harder,

do worse (in coverage sense)

9.5.3 Other Approaches for Means

Make just one “replication,” presumably to ameliorate initial bias

Point estimator of ν: average Y(m) of all the data Y₁, Y₂, ..., Y_m in the run Problem: How to estimate Var(Y(m)), needed to get c.i.’s, etc.?

Know one way not to do this:

(

) [

]

1

2 −

∑

=

m m m

Y Y

m

i i

Several methods to estimate Var(Y(m)):

Batch means

Time-series models Spectral analysis

Standardized time series

A different one-long-run approach (different point estimator):

Regenerative method

Two alternative modes of operation:

Fixed-sample size procedures — select run length m in advance, precision not controlled

Sequential procedures — prespecify precision, increase run length m as needed

Batch Means

Divide the run (of length m) into n contiguous “batches” of length k each (m = nk) Let Y_j(k) be the (batch) mean of the jth batch of size k

i

Y_i

k k k k k

Y ₁ Y ₂ Y ₃ Y ₄ Y ₅ m =nk Pretend: Y_j(k)’s are IID, unbiased for ν

Note: ( )

) (

1 1

m m Y

Y n

k Y

m

i i n

j j

=

=

∑

∑

, i.e., mean of batch means is the “grand” mean

Form “sample variance” among the batch means,

( )

1 ) ( ) ( )

( ¹

2

2

−

−

=

∑

=

n

m Y k Y n

S

n

j j Y

Approximate 100(1 – α)% confidence interval for ν:

n n t S

m

Y _n ^Y ( )

)

( ± _{−1,1−α/2}

(k) (k) (k) (k) (k)

Appeals of batch means:

Simple (relatively)

Often works fairly well (in terms of coverage)

Automatically implemented in some simulation software

Issues with batch means (see text for more details and references):

Choose batches big enough so thatY_j(k)’s are approximately uncorrelated Otherwise, S_Y²(n) can be biased (usually low) for Var(Y_j(k)), causing

undercoverage

How to choose batch size k? Equivalently, how many batches n?

Evidence: It may never pay to have more than n = 20 or 30 batches

Reason: Due to autocorrelation, splitting run into a larger number of smaller batches, while increasing degrees of freedom, degrades the quality

(variability) of each individual batch Generalizations of batch means

Overlapping batch means Separated batches

Weighting points within the batches

Sequential batch-means methods to control confidence-interval width Fix the number n of batches

Increase overall run length m

Increase batch size k proportionally to increase in m

Time-Series Models

Assume a relatively simple statistical model for the output process Discrete-time output process Y₁, Y₂, ... with steady-state mean ν

As with batch means, use overall average Y(m) as point estimator of ν Examples of statistical models for output process:

Autoregressive of order 2 (AR(2)):

Y_i = ν + φ1(Y_i–1 – ν) + φ2(Y_i–2 – ν) + εi, εi’s IID N(0, σ²) Autoregressive moving average of order (2, 1) (ARMA(2, 1)):

Y_i = ν + φ1(Y_i–1 – ν) + φ2(Y_i–2 – ν) + θ1εi–1 + εi, εi’s IID N(0, σ²) In general, for ARMA(p, q) model:

Use simulation output to identify (estimate) p and q — get pˆ and qˆ (several methods exist)

Estimate parameters via regression (least-squares) — (“fit the model”):

2

ˆ 2 1

ˆ 2 1

ˆ

,..., ˆ , ˆ ˆ

,..., ˆ , ˆ ˆ

σ

θ θ

θ

φ φ

φ

q p

Under this model, Var(Y(m)) = g(φ1, φ2, ..., φp, θ1, θ2, ..., θq, σ²) (the function g is known but messy)

Estimate Var(Y(m)) by g(φˆ₁,φˆ₂,...,φˆ_pˆ,θˆ₁,θˆ₂,...,θˆ_qˆ,σˆ²) Appeals of time-series models:

Physical intuition

Some theoretical support — most processes can be approximated by an ARMA(p, q) if we’re willing to admit large p and q

Issues with time-series models:

Spectral Analysis

Assume that output process is covariance-stationary:

Cov(Y_i, Y_{i + j}) = C_j, i.e., it depends only on j, and not on i Fairly mild assumption, especially after some warmup

Then

m

C m j C

m Y

m

j

∑

=

− +

=

1

1

0 (1 / )

)) ( ( Var

Estimate C_j, j = 1, 2, ..., m – 1 from output data:

( ) ( )

j m

m Y Y

m Y Y C

j m

i

j i i

j −

−

−

=

∑

= +

1

) ( )

( ˆ

Plug Cˆ ’s into formula for Var(_j Y(m)) in place of corresponding C_j’s Appeals of spectral analysis:

Using an exact variance formula

Relationship to other methods (overlapping batch means) Issues with spectral analysis:

For large m (as we’d expect), computationally burdensome to get all the Cˆ ’s — _j number of operations is about m²/2

Remedy: Computational device — Fast Fourier Transform

For j near the end (m – 1), Cˆ will be based on only a few terms, so will itself be _j highly variable—in fact, for j = m – 1, C^ˆ_m−1 =

(

)(

)

Remedy: Use different weights — spectral window — in sum for estimate of Var(Y(m))

Regenerative Method

Different approach to gathering data — do not use Y(m) as point estimator for ν Assume output process is regenerative:

Can identify a sequence of random indices 1 ≤ B₁ < B₂ < ... such that:

Starting from each index B_j, the process from then on follows the same probability distribution as does the process from any other B_j

The process from index B_j on is independent of the process before B_j This divides the process into IID cycles (or tours):

M

1 1

1 1

3 2

2

2 1

1

,..., ,

,..., ,

− +

− +

B B

B

B B

B

Y Y

Y

Y Y

Y

Most familiar example:

Single- or multiple-server queue

Any interarrival, service-time distributions

Y_i = D_i = delay in queue of ith arriving customer

B_j = index of the jth customer who arrives to find the whole system empty of other customers and all servers idle

For single-server case, a customer begins a regeneration cycle if and only if delay in queue is zero (not true for multiple servers)

i

D_i

B1=1 B2=6 B3=8 B4=13 B5=14 6=20

Comparable random variables defined on successive cycles are IID N_j = B_{j + 1} – B_j = length of jth cycle

∑

=

= + 1 1

j

j

B

B i

i

j Y

Z = sum of the observations in the jth cycle

Can show (renewal reward theory): ν = E(Zj)/E(N_j)

Make a simulation run of n′ cycles (don’t stop in the middle of a cycle) For j = 1, 2, ..., n′, let Uj = (Zj, Nj)^T; these are IID random vectors

Point estimator for ν:

) (

) ˆ (

n N

n Z

′

= ′

ν , which is biased, but strongly consistent

Variance estimation, confidence interval:

Let ΣΣ = 



 





22 12

12 11

σ σ

σ

σ be the covariance matrix of U_j

Define V_j = Z_j – ν Nj, which is unobservable, but has mean zero, variance

22 2 12 11

2 σ 2νσ ν σ

σ_V = − +

Central-limit theorem applied to V_j’s:

→

′ 

− ′

′ _D

V n

n N n

Z

2

) ( )

( σ

ν N(0, 1) as n′→ ∞

Need to estimate σ : _V²

Estimate ΣΣ by

[ ][ ]

1

) ( )

( )

( ˆ ) ( ˆ

) ˆ ( )

ˆ ( ₁

22 12

12 11

′−

− ′

− ′

 =



 





′

′

′

′

∑

=

n

n n

n n

n n

n

j

T j

j U U U

U σ

σ

σ σ

Let ^{( ) ˆ ( ) 2ˆ( )ˆ ( )}

[

]

12 11

2 n n n n n n

V ′ =σ ′ − ν ′σ ′ + ν ′ σ ′

σ , which is strongly

consistent for σ_V²

Can replace σ in above CLT by _V² σˆ_V²(n′) (continuous mapping thm.):

→

′ 

− ′

′ _D

V n

n N n

Z ˆ2

) ( )

( σ

ν N(0, 1) as n′→ ∞

Manipulate to get asymptotically valid 100(1 – α)% confidence interval for ν:

) ( ˆ ²

2 /

1 n n

z _V ′ ′

′ ± ⁻ σ

ν ^α

Appeals of regenerative method:

Firm mathematical foundation—asymptotic validity guarantees proper coverage probability (1 – α) as n′→ ∞

Issues with regenerative method:

Underlying process must be regenerative — not universally true

Must identify regeneration points (proof), code their recognition into program, gather data differently

Asymptotic validity depends on having a lot of cycles — if cycles tend to be long (as they often do in complicated models) we may be able to observe only a few cycles, and asymptotic validity doesn’t kick in

There are other ways to get a confidence interval (notably jackknifing; see text) Sequential-sampling versions have been developed—keep simulating more

cycles until the confidence interval is small enough

Standardized Time Series

Classical univariate statistics:

Take IID sample X₁, X₂, ..., X_n Want to estimate µ = E(Xi)

“Standardize” univariate data:

) 1 , 0 ( ) N

( of deviation standard

of Estimate

unknown )

( − →D

n X n

X µ

Basis for inference (confidence intervals, hypothesis tests, ...)

Observing a process (via simulation):

Observe the process Y₁, Y₂, ..., Y_m Want to estimate lim ( _i)

i E Y

∞

= →

ν

“Standardize” process data:

→

− D

i

Y Y

of i

deviation standard

of Estimate

unknown ν

Brownian bridge process (Brownian bridge process is fully understood, like N(0, 1)

Basis for inference (confidence intervals, hypothesis tests, ...)

Specifics:

Observe process Y₁, Y₂, ..., Y_m

Point estimator: Y(m), like batch means, time-series models, spectral analysis Form n batches of size k each (m = nk), as for batch means

For large m, Y(m) is approximately normal with mean ν and variance τ²/m where τ² = lim Var(Y(m))

m→∞

Let

_{∑ ∑∑} ( )

= = = + − 

 −

= − ⁿ

j k

s s

i

k j i

j k Y

k Y A k

1

2

1 1

) 1

312 ( ) (

As k → ∞, A/τ² → χ² distribution with n d.f., and is independent of Y(m) So for large k,

( )

) (

) ) (

(

2 2

mn A

m Y

n A

m m

Y ν

τ τ

ν = −

−

has an approximate t distribution with n d.f., so an asymptotically valid confidence interval for ν is Y(m)±t_n,1−α/2 A (mn)

Appeals of standardized time series:

Firm mathematical foundation—asymptotically valid intervals Assumptions are much weaker than for regenerative method Relatively simple

Issues with standardized time series:

As with all asymptotic, methods, how long is long enough?

Above is called “area” approach — A represents area under the standardized Brownian bridge

Other approaches look instead at the maximum attained by the standardized

Summary of Steady-State Estimation Procedures

Nothing will work well if computational budget is unduly limited

Batch means, spectral analysis, standardized time series display generally good performance (with respect to coverage probability, efficiency of data usage) Simplicity might argue against spectral analysis

Sequential-sampling versions exist for most of the main methods, to control confidence-interval width

Performance measures other than means (variances, proportions, quantiles) have been investigated

It may be difficult to implement these methods in the context of some existing simulation software, though some software does allow for, and even builds in and automates, some of these methods

9.5.4 Estimating Other Measures of Performance

Means do not always provide the one and only appropriate measure of performance

Probabilities: For some set B, estimate p = P(Y ∈ B), where Y is the steady-state random variable of the process Y₁, Y₂, ...

e.g., Y = delay of a message, B = (0, 5 minutes), so p is the probability that a message is delayed by less than 5 minutes

This is a special case of estimating means, if we define the indicator random variable







 ∈

= 0 otherwise if

1 Y B

Z , then p = P(Y ∈ B) = P(Z = 1) = 1×P(Z = 1) + 0×P(Z = 0) = E(Z)

Thus, can use all the methods described above for means to estimate a proportion

Quantiles: The q-quantile yq is the value such that P(Y ≤ yq) = q

e.g., Y = delay of a message, y_0.75 is the value below which are 75% of the message delays

See text for details on methods for estimating quantiles based on order statistics, batch means, and regenerative methods

9.6 Statistical Analysis for Steady-State Cycle Parameters

To dichotomize simulations into terminating vs. steady-state is a bit of an oversimplification:

Manufacturing model, has 8-hour shifts

Simulating in detail what goes on within a shift Performance might fluctuate widely within a shift

But what matters is the production (say) over the whole shift

Aggregate the data so that a “basic” observation Y_i is the production on the ith shift

Steady-state cycle (e.g., cycle = shift) parameter: steady-state mean of process defined over a cycle as a “unit” of time

Use above steady-state methods, except applied to random variables defined over a cycle, rather than individually at their finest level of detail

Example:

Production process, shift = 8 hours, lunch break 4 hours into each shift N_i = production in hour i; would expect it to drop in 4th hour; plot averaged

over 10 replications:

C

N = production in shift i; plot averaged over 10 replications: i

No steady- state

Seems to have a steady- state; use earlier methods

9.7 Multiple Measures of Performance

Usually: Want to observe several performance measures from a large simulation Average length of queue(s)

Maximum length of queue(s) Utilization(s)

Throughput(s) Difficulty:

Estimate each (expected) performance measure with a confidence interval Any of the intervals could “miss” its expected performance measure

Must be careful about overall statements of coverage (i.e., that all intervals contain their expected performance measures simultaneously)

Sometimes called the problem of multiple comparisons

Have k output performance measures, want overall (familywise) probability of at least 1 – α that the confidence intervals for them all contain their target expected performance measures

For s = 1, 2, ..., k, suppose the confidence interval for performance measure s is at confidence level 1 – αs

Then P(all intervals contain their respective performance measures) ≥ 1 –

∑

= k

s s 1

α (Bonferroni inequality)

Thus, pick αs’s so that

∑

= k

s s 1

Could pick αs = α/k for all s, or pick αs’s differently with smaller αs’s for the more important performance measures

Alternative Approach

Use multivariate statistical methods to get a joint confidence region (not necessarily a rectangle) that covers expected-performance-measure vector with probability 1 – α

Example:

k = 2 performance measures µ1 and µ2

Bonferroni approach: Separate intervals for µ1 and µ2

Multivariate approach: Ellipse containing (µ1, µ2) with probability 1 – α

µ1

µ₂

100(1−α1)% confidence interval for µ₁

100(1−α₂)% confidence interval for µ2

Bonferroni " box "

≥100(1−α)% confidence

region

100(1−α)%

confidence ellipse

Specific multivariate methods:

Multivariate batch means Multivariate spectral analysis Multivariate time-series methods Appeal of multivariate methods:

Smaller area (volume) with multivariate methods

9.8 Time Plots of Important Variables

So far, have concentrated on performance measures over the course of the whole simulation run

Averages Variances Extrema Proportions Quantiles

But these may mask important dynamic (within-run) behavior patterns Periodicities

Explosions Learning

Ways to pick up such dynamic behavior

Plot output processes (discrete or continuous) Animation