• 沒有找到結果。

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

N/A
N/A
Protected

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

Copied!
45
0
0

(1)

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

(2)

## 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

Simulation

model/code

### →

Random output:

Performance measures Thus, the output performance measures:

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

(3)

The Statistical Nature of Simulation Output

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

Yi = production in ith hour in a factory

Yi’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 Yi’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, y11, y12, ..., y1m

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

y21, y22, ..., y2m

. . .

yn1, yn2, ..., ynm

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

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

(4)

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

(5)

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.?

(6)

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”

Si = 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

### ∑

i

= =1

) (

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

### ∑

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:

Di = delay of ith customer in queue

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

(7)

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

### ∫

T

= 0 ( ) )

(

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

### ∫

T

= 0 (2, )

2

)) ( ( )

(

(8)

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

### ∫

T

= 0 ( ) )

(

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

(9)

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)

(10)

## 9.2 Transient and Steady-State Behavior of a Stochastic Process

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

Let Fi(y | I) = P(Yi y | I) be the transient (cumulative) distribution of the process at (discrete) time i

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

If there is a CDF F(y) such that Fi(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 Fi(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 Yi’s could still have large variance and thus “bounce around” a lot — they’re just not systematically trending any more

(11)

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):

(12)

## 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

(13)

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

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

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

(14)

## 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 Xj 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 X1, X2, ..., Xn 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)

(15)

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

Example: Delays in a queue of individual jobs: D1, D2, D3, ..., Dm Want to estimate µ = E(average delay of the m jobs)

Sample mean D m D m

m

i

### ∑

i

=

=

1

)

( is an unbiased estimator for µ

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

But “sample variance”

( )

( 1)

1

2

### ∑

=

m m m

D D

m

i

i may be severely biased for

Var(D(m)) Reason:

Corr(Di, Di + 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(Di, Di + l) > 0

Causes





### ∑

− −

=

) 1 ( )

(

1

2 m m

m D D E

m

i

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

Di + 1 is pretty much the same as Di

Di’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

(16)

### 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”

(17)

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 Xj be the performance measure from the jth replication

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

Xj = utilization of a bottleneck machine

Then X1, X2, ..., Xn are IID and unbiased for µ = E(Xj)

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

n X n

X

n

j

### ∑

j

= =1

)

( is an unbiased estimator of µ

1 )) ( (

) (

2 2

=

### ∑

n

n X X n

S

j

is an unbiased estimator of Var(Xj)

n n t S

n

X n ( )

)

( ± 1,1α/2 covers µ with approximate probability 1 – α

(tn–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)

(18)

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

Xj’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, S2(10) = 3.96, t9, 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

(19)

Why “approximate” 90% confidence interval?

Assumes Xj’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: 1

2

3

Components fail independently

Times Ti to failure in components ~ Weibull(α = 0.5, β = 1.0) Time to system failure = T = min(T1, max(T2, T3))

µ = E(time to system failure) = 0.78

n Estimated coverage (want 90%) 5

10 20 40

71%

75%

80%

84%

(20)

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 tn 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”

(21)

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

(22)

### 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

(23)

Variances (or Standard Deviations)

Interested in process variability Xj = 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, Xj = 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(Xj x) = 0.98

One approach (more details in text):

Make n replications, observe X1, X2, ..., Xn Sort the Xj’s into increasing order

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

(24)

### 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

(25)

## 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 X1

X2

X3

X4

X5

X1

Many short runs One long run

Good Simple (same as terminating) Get IID data

Less point-estimator bias No restarts

(26)

### 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 t0 (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

(27)

### 9.5.2 Replication/Deletion Approaches for Means

Assume that an appropriate warmup period has been determined

Xj = 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

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)

(28)

### 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 Y1, Y2, ..., Ym 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)

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

(29)

Batch Means

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

i

Yi

k k k k k

Y 1 Y 2 Y 3 Y 4 Y 5 m =nk Pretend: Yj(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 ) ( ) ( )

( 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)

(30)

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 thatYj(k)’s are approximately uncorrelated Otherwise, SY2(n) can be biased (usually low) for Var(Yj(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

(31)

Time-Series Models

Assume a relatively simple statistical model for the output process Discrete-time output process Y1, Y2, ... 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)):

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

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

Use simulation output to identify (estimate) p and q — get and (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, σ2) (the function g is known but messy)

Estimate Var(Y(m)) by g(φˆ1,φˆ2,...,φˆpˆ,θˆ1,θˆ2,...,θˆqˆ,σˆ2) 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:

(32)

Spectral Analysis

Assume that output process is covariance-stationary:

Cov(Yi, Yi + j) = Cj, 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

### ∑

j

=

− +

=

1

1

0 (1 / )

)) ( ( Var

Estimate Cj, 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 Cj’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 m2/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ˆm1 =

Y1Y(m)

YmY(m)

### )

, which is based on only one term in the sum

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

(33)

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 ≤ B1 < B2 < ... such that:

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

The process from index Bj on is independent of the process before Bj 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

Yi = Di = delay in queue of ith arriving customer

Bj = 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

Di

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

(34)

Comparable random variables defined on successive cycles are IID Nj = Bj + 1 – Bj = 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(Nj)

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

(35)

Variance estimation, confidence interval:

Let ΣΣ =

 

22 12

12 11

σ σ

σ

σ be the covariance matrix of Uj

Define Vj = Zjν Nj, which is unobservable, but has mean zero, variance

22 2 12 11

2 σ 2νσ ν σ

σV = − +

Central-limit theorem applied to Vj’s:

→

′ 

− ′

D

V n

n N n

Z

2

) ( )

( σ

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

Need to estimate σ : V2

Estimate ΣΣ by

1

) ( )

( )

( ˆ ) ( ˆ

) ˆ ( )

ˆ ( 1

22 12

12 11

′−

− ′

− ′

 =

 

### ∑

=

n

n n

n n

n n

n

j

T j

j U U U

U σ

σ

σ σ

Let ( ) ˆ ( ) 2ˆ( )ˆ ( )

ˆ( )

### ]

ˆ22( ) 2

12 11

2 n n n n n n

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

σ , which is strongly

consistent for σV2

Can replace σ in above CLT by V2 σˆV2(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

2 /

1 n n

z V ′ ′

′ ± σ

ν α

(36)

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

(37)

Standardized Time Series

Classical univariate statistics:

Take IID sample X1, X2, ..., Xn 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 Y1, Y2, ..., Ym 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, ...)

(38)

Specifics:

Observe process Y1, Y2, ..., Ym

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 τ2/m where τ2 = lim Var(Y(m))

m

Let

### ∑ ∑∑ ( )

= = = + 

 −

= − n

j k

s s

i

k j i

j k Y

k Y A k

1

2

1 1

) 1

312 ( ) (

As k → ∞, A/τ2χ2 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(mtn,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

(39)

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

(40)

### 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 Y1, Y2, ...

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) + 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, y0.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

(41)

## 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 Yi 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

(42)

Example:

Production process, shift = 8 hours, lunch break 4 hours into each shift Ni = 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

Seems to have a steady- state; use earlier methods

(43)

## 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

(44)

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

µ2

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

100(1−α2)% 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

(45)

## 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

It is important to allow for all students to add their ideas to the story so giving each student an area of responsibility to add to the story recipe can help prompt this. For

Now, nearly all of the current flows through wire S since it has a much lower resistance than the light bulb. The light bulb does not glow because the current flowing through it

• Examples of items NOT recognised for fee calculation*: staff gathering/ welfare/ meal allowances, expenses related to event celebrations without student participation,

This kind of algorithm has also been a powerful tool for solving many other optimization problems, including symmetric cone complementarity problems [15, 16, 20–22], symmetric

For problems 1 to 9 find the general solution and/or the particular solution that satisfy the given initial conditions:. For problems 11 to 14 find the order of the ODE and

If we want to test the strong connectivity of a digraph, our randomized algorithm for testing digraphs with an H-free k-induced subgraph can help us determine which tester should

Performance metrics, such as memory access time and communication latency, provide the basis for modeling the machine and thence for quantitative analysis of application performance..

• However, inv(A) may return a weird result even if A is ill-conditioned, indicates how much the output value of the function can change for a small change in the