Simulation Output Analysis Bruno Tuffin INRIA Rennes - Bretagne - - PowerPoint PPT Presentation

Simulation Output Analysis

Bruno Tuffin

INRIA Rennes - Bretagne Atlantique

PEV: Performance EValuation M2RI - Networks and Systems Track Rennes

Bruno Tuffin (INRIA) Output Analysis PEV - 2010 1 / 47

Outline

1

Introduction

2

Very basic tool: Central Limit Theorem

3

Accuracy and Efficiency of an estimator

4

Performance measures

5

Transient Simulations

6

Steady-state Simulations

7

Acceleration techniques

8

Some open problems

Bruno Tuffin (INRIA) Output Analysis PEV - 2010 2 / 47

References

1

• S. Asmussen and P.W. Glynn. Stochastic Simulation. Algorithms and Analysis.

Stochastic Modelling and Applied Probability Series, Springer Verlag, 2007.

2

• C. Alexopoulos and A.F. Seila. Advanced Methods for Simulation Output Analysis.

In the Proceedings of the 1998 Winter Simulation Conference, D.J. Medeiros, E.F. Watson, J.S. Carson and M.S. Manivannan, eds. 1998.

3

• M. Nakayama. Simulation Output Analysis. In the Proceedings of the 2002 Winter

Simulation Conference, E. Yucesan, C.H. Chen, J.L. Snowdown and J.M. Charnes,

• eds. 2002.

4

• B. Schmeiser. Simulation Output Analysis: a Tutorial Based on one research
• Thread. In the Proceedings of the 2004 Winter Simulation Conference, R.G.

Ingalls, M.D. Rossetti, J.S. Smith and B.A. Peters, eds. 2004.

5

• G. Rubino and B. Tuffin. Simulations et m´

ethodes de Monte Carlo. In : Techniques de l’Ing´ enieur, 2007.

6

• B. Tuffin. La simulation de Monte Carlo. Editions Herm`

es, 2010.

Bruno Tuffin (INRIA) Output Analysis PEV - 2010 3 / 47

Outline

1

Introduction

2

Very basic tool: Central Limit Theorem

3

Accuracy and Efficiency of an estimator

4

Performance measures

5

Transient Simulations

6

Steady-state Simulations

7

Acceleration techniques

8

Some open problems

Bruno Tuffin (INRIA) Output Analysis PEV - 2010 4 / 47

Introduction: what should we get from simulation results

Analytical results provide an exact result for a performance measure; Numerical analysis techniques provide an appraoch result, and hopefully an idea of the error; Simulation make use of random numbers:

◮ 2 different simulations give 2 different results; ◮ We need, and there exist efficient tools providing an idea of the error.

A simulation just giving a number as result is disappointing: no

• utput analysis.

◮ Example: opinion polls in media. Bruno Tuffin (INRIA) Output Analysis PEV - 2010 5 / 47

What we should get: a confidence interval

Assume that we want to estimate a performance measure µ. While standard numerical analysis provide a deterministic approximation for µ the quantity and (potentially) a strict error bound... ... simulation does provide a Confidence Interval (A, B) with confidence level 1 − α:

◮ It means that we can only say that µ is in (A, B) with probability

1 − α;

◮ No insurance that it is true (100α% chances to be out of it). ◮ We can increase the confidence level 1 − α, but at the expense of the

interval width B − A.

◮ Usually, the more we simulate, the smaller the interval width.

The goal of this course is to explain how to build such intervals.

Bruno Tuffin (INRIA) Output Analysis PEV - 2010 6 / 47

Outline

1

Introduction

2

Very basic tool: Central Limit Theorem

3

Accuracy and Efficiency of an estimator

4

Performance measures

5

Transient Simulations

6

Steady-state Simulations

7

Acceleration techniques

8

Some open problems

Bruno Tuffin (INRIA) Output Analysis PEV - 2010 7 / 47

Basic tool in statistics: Central Limit Theorem (CLT)

Most (if not all) performance measures µ can be expressed as µ = E[X] for some random variable X (ex: mean throughput, mean delay or mean number of packets at a router...) Central Limit Theorem: Let X1, · · · , Xn n i.i.d copies of X with (finite) expectation µ and variance σ2. The CLT gives the behaviour

• f arithmetical average

¯ Xn = 1 n

n

• i=1

Xi when n is large. More precisely, normalized r.v. √n ¯ Xn − µ σ (so that mean is 0 and variance 1) has a distribution N(0, 1) (Gaussian law with mean 0 and variance 1) as n → ∞.

Bruno Tuffin (INRIA) Output Analysis PEV - 2010 8 / 47

Confidence interval for µ

From tables provided in textbooks (or on the web), we can get values z1−α/2 such that I P[N(0, 1) ≤ z1−α/2] = 1 − α/2. Then 1 − α = I P[−z1−α/2 ≤ N(0, 1) ≤ z1−α/2] ≈ I P

• −z1−α/2 ≤ √n

¯ Xn − µ σ ≤ z1−α/2

• =

I P

• ¯

Xn − z1−α/2 σ √n ≤ µ ≤ ¯ Xn + z1−α/2 σ √n

• This gives for µ a confidence interval at confidence level 1 − α
• ¯

Xn − z1−α/2 σ √n, ¯ Xn + z1−α/2 σ √n

• .

Standard values for

◮ 1 − α = 90% gives z1−α/2 = 1.64, ◮ 1 − α = 95% gives z1−α/2 = 1.96, ◮ 1 − α = 99% gives z1−α/2 = 2.58. Bruno Tuffin (INRIA) Output Analysis PEV - 2010 9 / 47

Properties and accuracy

Confidence interval width 2z1−α/2 σ √n decreasing at rate O((n−1/2). This is much better than numerical analysis in high dimension for which it is often n−k/d (for some k) for d-dimensional problems. Increasing the sample size n improves the accuracy, but at this rate O(n−1/2).

Bruno Tuffin (INRIA) Output Analysis PEV - 2010 10 / 47

But σ2 not known in general: how to proceed?

Major issue: variance σ2 generally unknown. σ2 estimated by the (unbiased) sample variance S2(n) = 1 n − 1

n

• i=1

(Xi − ¯ Xn)2 = 1 n − 1 n

• i=1

X 2

i − n(¯

Xn)2

• .

Easy to implement, for the simulation, we only need two counters, one for n

i=1 Xi, and the other for n i=1 X 2 i .

This gives for µ a confidence interval at confidence level 1 − α

• ¯

Xn − z1−α/2 S(n) √n , ¯ Xn + z1−α/2 S(n) √n

• .

If n is moderate, z1−α/2 can be found from Student distribution with n − 1 degrees fo freedom (which converges to a Gaussian law quite quickly). Textbooks provide such values also.

Bruno Tuffin (INRIA) Output Analysis PEV - 2010 11 / 47

Example: simulating the mean number of customer at fixed time T for an M/M/1 queue

Consider the simulation of an M/M/1 queue up to time T, starting from an empty queue at t = 0 Simulation is realized by using discrete-event simulation (see previous course), by schedulling arrivals and departures and keeping track of the state given by the number of customers in the queue. A run is a single simulation up to T, and gives X: the number of customers in the queue at time T. We realize n i.i.d. copies of X, X1, . . . , , Xn, and use the above framework.

Bruno Tuffin (INRIA) Output Analysis PEV - 2010 12 / 47

Outline

1

Introduction

2

Very basic tool: Central Limit Theorem

3

Accuracy and Efficiency of an estimator

4

Performance measures

5

Transient Simulations

6

Steady-state Simulations

7

Acceleration techniques

8

Some open problems

Bruno Tuffin (INRIA) Output Analysis PEV - 2010 13 / 47

Total error

Statistical quality of an estimator ¯ Xn of µ is measured by the mean squared error MSE(¯ Xn, µ) = E[(¯ Xn − µ)2] = Bias2( ¯ Xn, µ) + σ2( ¯ Xn). Indeed, it may happen that E[¯ Xn] = µ, because

1

some good estimators are inherently biased,

2

it is much ”cheaper” to sample from a close but not exact distribution

3

modelling errors...

In many cases though, Bias(¯ Xn, µ) = 0.

Bruno Tuffin (INRIA) Output Analysis PEV - 2010 14 / 47

Efficiency

It may appear that variance is THE measure of efficiency of an estimator: the best estimator is the one with smallest variance (for same sample size). However, we also need to consider CPU times. The problem is rather: what is , between ¯ X and ¯ X ′ the estimator with the smallest variance in a computational budget c? If T and T ′ are the times required to generate one independent replication of X and X ′ when computing ¯ X and ¯ X ′, the number of replications will be respectively n = c/T and n′ = c/T ′. Thus the best estimator is ¯ X if σ2(X)T < σ2(X ′)T ′. Note that σ2(X)T can be interpreted as the variance for a unit of time.

Bruno Tuffin (INRIA) Output Analysis PEV - 2010 15 / 47

Sequential procedure for absolute error criterion

Assume that we want to have an estimator ¯ Xn of µ such that | ¯ Xn − µ| < ε with probability 1 − α. From the CLT: ε ≈ z1−α/2 σ √n so that we need n >

• (z1−α/2)2σ2

ε2

• .

Since σ2 unknown and estimated, we need to proceed as follows

◮ Use a sample of size n0 (typically n0 ≥ 50), estimate σ2, then compute

n = (z1−α/2)2S(n0)2

ε2

• .

◮ Then generate a sample of size n, for which absolute error should

approximately be ε.

Can also be performed in a sequential way.

Bruno Tuffin (INRIA) Output Analysis PEV - 2010 16 / 47

Sequential procedure for relative error criterion

Assume that we want to have an estimator ¯ Xn of µ such that | ¯ Xn − µ| < ε|µ| with probability 1 − α. From the CLT: ε|µ| ≈ z1−α/2 σ √n so that we need n >

• (z1−α/2)2σ2

µ2ε2

• .

Since σ2 unknown and estimated, we need to proceed as follows

◮ Use a sample of size n0 (typically n0 ≥ 50), estimate σ2, then compute

n = (z1−α/2)2S(n0)2

¯ Xn0ε2

• .

◮ Then generate a sample of size n, for which absolute error should

approximately be ε.

Can also be performed in a sequential way.

Bruno Tuffin (INRIA) Output Analysis PEV - 2010 17 / 47

Outline

1

Introduction

2

Very basic tool: Central Limit Theorem

3

Accuracy and Efficiency of an estimator

4

Performance measures

5

Transient Simulations

6

Steady-state Simulations

7

Acceleration techniques

8

Some open problems

Bruno Tuffin (INRIA) Output Analysis PEV - 2010 18 / 47

Performance measures: transient or steady-state

The performance measures of interest can be decomposed into two sets

◮ Transient measures, which look at a system evolving over a finite time.

Ex: in a network, we may be interested in the delay at a particular time, or over a fixed interval of time.

◮ Steady-state measures which look at a system evolving over an infinite

time horizon. It typically does not depend upon the initial conditions. Ex: we are interested in the mean delay at a network in the long run.

Bruno Tuffin (INRIA) Output Analysis PEV - 2010 19 / 47

Outline

1

Introduction

2

Very basic tool: Central Limit Theorem

3

Accuracy and Efficiency of an estimator

4

Performance measures

5

Transient Simulations

6

Steady-state Simulations

7

Acceleration techniques

8

Some open problems

Bruno Tuffin (INRIA) Output Analysis PEV - 2010 20 / 47

Transient Simulations (also called finite-horizon simulations)

Transient simulations are performed using the previously described techniques based on the CLT, depending on whether we have

◮ a pre-specified sample size n ◮ a pre-spcecified absolute or relative error ε ◮ a pre-specified CPU budget c.

We thus just apply the CLT on independent runs using a sample of size n computed according to the selected criterion.

Bruno Tuffin (INRIA) Output Analysis PEV - 2010 21 / 47

Outline

1

Introduction

2

Very basic tool: Central Limit Theorem

3

Accuracy and Efficiency of an estimator

4

Performance measures

5

Transient Simulations

6

Steady-state Simulations

7

Acceleration techniques

8

Some open problems

Bruno Tuffin (INRIA) Output Analysis PEV - 2010 22 / 47

General comments on steady-state simulations

Things are a little more complicated if we want to deal with steady-state measures. Indeed, we cannot sample during an infinite time, we have to got to a time T, meaning that results are biased (so-called initial bias). Ex: simulation of a queue. Several options are possible (explained below)

◮ Replication-deletion techniques ◮ Regenerative techniques ◮ Batch Means techniques ◮ Standardized Time Series techniques Bruno Tuffin (INRIA) Output Analysis PEV - 2010 23 / 47

Stochastic processes at hand

There are two cases A discrete time process X = (Xi)i∈I

N, output process with

inter-valued time index, such that ¯ Xm = 1 m

m−1

• i=0

Xi → µ as m → ∞.

◮ Ex: Xi Number of users of a wifi spot in day i.

A continuous time process X = (Xt)t≥0, output process with inter-valued time index, such that ¯ Xt = 1 t t Xsds → µ as t → ∞.

◮ Ex: Xt number of customers in a queue at time t. ◮ How is is computed in practice? Just look at instants

t0 < t1 < · · · < tn = t of change of state and compute ¯ Xt = 1 t

m−1

• i=0

(ti+1 − ti)X +

ti .

Bruno Tuffin (INRIA) Output Analysis PEV - 2010 24 / 47

On getting a confidence interval for a long run

We know that the estimator converges if we use a long run (m or t → ∞). But how do we get a confidence interval? How to compute it? For almost all systems possessing a unique steady-state, there is a contral limit theorem: there exists ˜ σ (called time-average variance constant) such that √m ¯ Xm − µ ˜ σ has a distribution N(0, 1) as m → ∞. This is valid for the discrete and the continuous process (this last one by replacing m by t). From the CLT, we can get a confidence interval (similarly to the first slides). But how to have an estimation of ˜ σ?

Bruno Tuffin (INRIA) Output Analysis PEV - 2010 25 / 47

Replication-deletion techniques

We can use independent replications of shorter runs: instead of m steps, we use n′ replications of m′ = m/n′ steps. The idea is then to use the standard CLT over the m′ replications of average value X ′ =

1 m′

m′−1

i=0

Xi (independent replications X ′

j of that

average, leading to estimator 1

n′

n′

j=1 X ′ j ).

Easy estimation of variance, an good confidence interval, but larger initial bias, because of shorter runs... A way to attenuate the problem is to delete the first ℓ steps of each run, because too far from steady-state, and rather use X ′ = 1 m′ − ℓ

m′−1

• j=ℓ

Xi. But then n′ℓ data unused with respect to a single run... not that efficient.

Bruno Tuffin (INRIA) Output Analysis PEV - 2010 26 / 47

Regenerative techniques

Method developped in the 70s. Assume that there are (random) indices T0 < T1 < · · · such that portion (XTi+j)≥0 of the process have the same distribution ∀i, and are independent of the portion of the process prior to Ti. Suche a process is called a regenerative process, and time Ti a regeneration time. Usual case: if (Xn)n≥0 is a Markov chain, then entrance time into a given state is a regeneration time

◮ An M/M/1 queue is a regenerative process: entrance in any specific

fixed state give regeneration times.

◮ Even an M/G/1 queue is a regenerative process, entrance time into the

state ”empty queue” is a regeneration time.

A portion between two successive regeneration epochs is called a cycle.

Bruno Tuffin (INRIA) Output Analysis PEV - 2010 27 / 47

Regenerative techniques (2)

Start at T0 = 0 from regenerative state. Define Yi = Ti+1−1

Ti

Xj and Zi = Ti+1 − Ti and assume that E[Zi] < ∞, i.e. the mean regeneration time is finite. We then have µ = E[Y1] E[Z1]. The idea is to simulate the process (Xi)i≥0 over n cycles and compute ¯ µn = ¯ Yn ¯ Zn . Note that developping ¯ µ, we have ¯ Yn = Tn−1

i=0

Xi/n and ¯ Zn = Tn−1

i=0 (Ti+1 − Ti)/n = (Tn − 1)/n, leading to

¯ µn = Tn−1

i=0

Xi Tn − 1 , that is the long run estimator.

Bruno Tuffin (INRIA) Output Analysis PEV - 2010 28 / 47

Regenerative techniques (3)

Cycles are of particular interest to build a confidence interval: √n ¯ µn − µ σ′/¯ Zn is asymptotically distributed as N(0, 1), where (σ′)2 = Var(Y1) − 2µCov(Y1, Z1) + µ2Var(Z1). Proof: look at the CLT applied to independent r.v. Yi − µZi. (Exercise.) Parameter σ′ unknown, but easily estimated by

◮ µ ≈ ¯

µn,

◮ Var(Y1) ≈

1 n−1

n

i=1(Yi − ¯

Yn)2,

◮ Var(Z1) ≈

1 n−1

n

i=1(Zi − ¯

Zn)2,

◮ Cov(Y1, Z1) ≈

1 n−2

n

i=1(Yi − ¯

Yn)(Zi − ¯ Zn).

On the other hand, the method requires a sufficiently short regeneration time to have enough cycles for building the confidence

• interval. This is not always the case (multidimensional states...)

Bruno Tuffin (INRIA) Output Analysis PEV - 2010 29 / 47

Batch Means techniques: general idea

Basic principle: it is natural to think that as k grows, Xi and Xi+k will be less and less correlated, close to independent. Therefore, the sequence of observations X1, . . . Xn of a long-run simulation is divided in blocks that are used to estimate the variance (and µ). Decompose the process in groups of size b, i-th group being made of X(i−1)b, . . . , Xib−1. The average for i-th block is Yi(b) = 1 b

b−1

• j=0

X(i−1)b+j.

Bruno Tuffin (INRIA) Output Analysis PEV - 2010 30 / 47

Batch Means techniques: general idea (2)

If the process is in steady-state, Cov(Yi(b), Yi+ℓ(b)) = 1 b

b−1

• j=1−b
• 1 − |j|

b

• Cov(X0, Xℓb+j)

does ot depend on i. If Cov(Yi(b), Yi+ℓ(b)) → 0 as b → ∞, the Yi(b) are almost independent (and normally distributed) for b sufficiently large. the CLT can be applied on the Yi(b) to estimate µ. To determine and adapt the size of blocks, correlation tests can be applied.

Bruno Tuffin (INRIA) Output Analysis PEV - 2010 31 / 47

Some rules to define batch sizes

Fixed number of sizes: fix the number of batches to k and use ⌊n/k⌋. It is known to converge, but has some limitations

◮ A tendency to yield a wider confidence interval ◮ statistical fluctuations for the confidence interval d not diminish

relative to the statistical fluctuation in the sample mean.

Square root (SQRT) rule: use b = ⌊√n⌋ and sample size ⌊√n⌋. Limitation:

◮ A tendency to under-estimate the variance Bruno Tuffin (INRIA) Output Analysis PEV - 2010 32 / 47

Dynamic rules for batching

LBATCH procedure: at time nl, if a correlation test is negative between batches, we then use the fixed number of batches rule and so

• n, up to when it fails, then the SQRT rule is applied for all future

reviews. ABATCH procedure: at time nl, if a correlation test is negative between batches, we then use the fixed number of batches rule; when it fails the SQRT rule is applied... Otherwise, we can use overlapping batch means (less sensitive to batch size for instance).

Bruno Tuffin (INRIA) Output Analysis PEV - 2010 33 / 47

Standardized Time Series method

Define series Tn(t) = ⌊nt⌋( ¯ Xn − ¯ X⌊nt⌋ σ√n , 0 ≤ t ≤ 1 with ⌊x⌋ interger part of x. Under some stationnarity assumptions (Xi and Xi+j being close to independent for j large), couple (√n( ¯ Xn − µ), σTn) converges in law to (σW (1), σB) with B = (B(t))t≥0 standard Brownian bridge and W standard Brownian motion.

◮ Remind that standard Brownian motion is W = (W (t))t≥0 such that

(a) if s ≥ 0 and t > 0, X(s + t) − X(s) is N(0, t); (b) forall all intervals (t1, t2], . . . , ]t2k−1, t2k], increments X(t2) − X(t1), . . . , X(t2k) − X(t2k−1) are mutually independent.

◮ Standard Brownian bridge: a standard Brownian motion, but

condtionally to B(0) = B(1) = 0.

Bruno Tuffin (INRIA) Output Analysis PEV - 2010 34 / 47

Standardized Time Series method (2)

If A = 1

0 σB(t)dt, E[A2] = σ2/12 and σ2 can be computed by

estimating E[A2]:

◮ Separate X1, . . . , Xn in k = n/b blocks of size b. ◮ For n sufficiently large, the Ai = b−1

j=0 ((n + 1)/2 − j)X(i−1)b+j

(i = 1, . . . , k) are approximately i.i.d and normally distributed.

◮ An estimator of E[A2] is

ˆ S2

A =

1 (b3 − b)k

k

• i=1

A2

i .

We then use the CLT: ¯ Yk − µ

• 12ˆ

S2

A/n

converges in law to N(0, 1) as k → ∞, with ¯ Yk average over k blocks Yi(b) = b−1 b−1

j=0 X(i−1)b+j.

Easier to implement that batch means technique, with better asymptotic properties. But may require much longer simulation times.

Bruno Tuffin (INRIA) Output Analysis PEV - 2010 35 / 47

Outline

1

Introduction

2

Very basic tool: Central Limit Theorem

3

Accuracy and Efficiency of an estimator

4

Performance measures

5

Transient Simulations

6

Steady-state Simulations

7

Acceleration techniques

8

Some open problems

Bruno Tuffin (INRIA) Output Analysis PEV - 2010 36 / 47

Acceleration techniques: what can be done

There are different ways to improve the efficiency of an estimator:

◮ to reduce the variance of a single replication ◮ or to reduce the average time to generate such a replication.

Many variance reduction techniques exist, but average time reduction can as efficient and is less known. Ex: Reliability analysis of a network (connection between two nodes s and t) using disjoint paths from s to t. Method developped in Dionysos group.

Bruno Tuffin (INRIA) Output Analysis PEV - 2010 37 / 47

Control variates

When we want to estimation µ = E[X], this can be re-written as µ = E[Z] + E[X − Z] for some r.v. Z. if Z is simple enugh so that its expectation is known, we can only deal with the second expectation. IF Z is close enough to X, its variations follows those of X and this can lead to a drastic variance reduction.

Bruno Tuffin (INRIA) Output Analysis PEV - 2010 38 / 47

Stratified sampling

The idea is to separate the sampling domain in sub-domains, and simulate

• ver those domains in proportion to their probability.

Let λ distribution of r.v. X over the set of possibilities Ω. Ω is partitionned into Ω = ∪k

j=1Cj.

Over each Cj (1 ≤ j ≤ k), we use a sample X (j)

1 , · · · , X (j) nj

µj-distributed, where µj is λ conditionnally to being on Cj. µ = E[X] =

k

• j=1
• Cj

xdλ(x) =

k

• j=1

λ(Cj)

• Cj

xdµj(x) ≈

k

• j=1

λ(Cj) nj

nj

• i=1

X (j)

i

. Variance

k

• j=1

λ(Cj) nj

• Cj
• x −

1 λ(Cj)

• Cj

xdλ(x)

• 2

dλ(x). It is always smaller if nj = λ(Cj)n.

Bruno Tuffin (INRIA) Output Analysis PEV - 2010 39 / 47

Antithetic variates

In the simplest setting: Assume that we want to find E[f (X)] with X uniform r.v. over [0, 1] and f monotone function. Then Z = 1 2(f (X) + f (1 − X)) is such that σ2(Z) < σ2(f (X))/2. This principle can be generalized to multi-dimensional cases. A last technique, probably the main used, importance sampling, is briefly described in the case of rare events.

Bruno Tuffin (INRIA) Output Analysis PEV - 2010 40 / 47

Rare-event simulation: context

Rare events appear and are critical in many contexts:

◮ reliability analysis (ex: computer systems, a network, a nuclear plant) ◮ finance ◮ ruin of an insurance company ◮ collisions in air traffic management ◮ Telecommunications (loss of some kind of packets...) Bruno Tuffin (INRIA) Output Analysis PEV - 2010 41 / 47

Rare-event simulation: difficulties

With crude simulation technique, the resulting estimatoris highly unreliable (almost useless) when the probability γ is very small. For example, if µ = 10−10 and if we want the expected number of

• ccurrences of this event to be at least 100, we must take n = 1012

(a huge number). For n < 1010, we are likely to observe not even a single occurrence of this event. In this case, not only the estimator of γ takes the value 0 but the empirical variance as well, leading to a confidence interval (0, 0)!

Bruno Tuffin (INRIA) Output Analysis PEV - 2010 42 / 47

Robutsness of rare event estimators

The estimator ¯ Xn has the binomial distribution with parameters (n, µ). Its relative error is RE[ ¯ Xn] = (Var[¯ Xn])1/2 µ = (µ(1 − µ))1/2 nµ ≈ µ−1/2 n , which increases toward infinity when µ → 0. An alternative unbiased estimator of µ, say ˜ Xn, is said to have bounded relative error if limµ→0 RE[ ˜ Xn] < ∞. This implies that lim

µ→0

log(E[ ˜ X 2

n ])

log µ = 2. When the latter (weaker) condition holds, the estimator ˜ Xn is said to be asymptotically efficient.

Bruno Tuffin (INRIA) Output Analysis PEV - 2010 43 / 47

Rare-event simulation: a glimpse at solutions (1)

Many different ways have been proposed to deal with rare events. Two main families: Importance sampling: it consists in changing the probability laws that drive the evolution of the system, to increase the probability of the rare event, and multiplying the estimator by an appropriate likelihood ratio so that it remains unbiased.

◮ In the setting of computing an integral, it just consists in a chnage of

variable µ = Ef [h(X)] =

• h(x)f (x)dx =
• h(x) f (x)

g(x)g(x)dx = Eg[h(X)f (X)/g(X)].

◮ Ex: for an M/M/1 queue, if trying to look at the probability that the

buffer exceeds a given level, just switch arrival and service rates.

Bruno Tuffin (INRIA) Output Analysis PEV - 2010 44 / 47

Rare-event simulation: a glimpse at solutions (2)

Splitting: the probability laws of the system remain unchanged, but an artificial drift toward the rare event is created by terminating the trajectories that seem to go away from it and by splitting (cloning) those that are going in the right direction.

B_3 B_1 A=B_4 B_2

Bruno Tuffin (INRIA) Output Analysis PEV - 2010 45 / 47

Outline

1

Introduction

2

Very basic tool: Central Limit Theorem

3

Accuracy and Efficiency of an estimator

4

Performance measures

5

Transient Simulations

6

Steady-state Simulations

7

Acceleration techniques

8

Some open problems

Bruno Tuffin (INRIA) Output Analysis PEV - 2010 46 / 47

Some open problems in the area

Batching techniques require much more attention. Rare event simulation: finding efficient techniques for new kind of problems (ex: heavy tailed queues). In simulation in general: to deal with correlations, and dependent random variables. · · ·

Bruno Tuffin (INRIA) Output Analysis PEV - 2010 47 / 47