1 Organisatie Essay College 1: Introductie Monte Carlo Simulatie Discrete Event Simulatie Java A.A.N. Ridder Department EOR Vrije Universiteit Amsterd...
This is called a random sample with sample size N;
I
The mean of the associated H(Xi ) variables is a random variable; i.e., a function on the sample space; called sample mean:
I
I
Make use of the last property;
I
By executing the simulation program, you generate observations (or data; realizations; values; ...) x1 , x2 , . . . , xN of the random sample;
I
Compute the average of the data: `=
N 1 X H(xi ); N i=1
I
Then (SLLN): ` ≈ `: we say that ` is an estimate of `;
N 1 X `b = H(Xi ). N i=1
I
Question: how good is the estimate?
I
Actually; it is custom in the simulation literature to just consider and analyse the estimators;
E`b = `; i.e., `b is an unbiased estimator of the performance measure `;
I
Thus; how good is the sample mean estimator?
Properties: I
I
For large N, according to SLLN, `b is not so random; i.e., it is almost a constant function; i.e., `b ≈ ` for almost all samples.
c Ad Ridder (VU)
Essay – January 2016
11 / 83
c Ad Ridder (VU)
Essay – January 2016
12 / 83
Statistics
Using the CLT
Let σ 2 =
Var[H(X)] be the variance of the output variable;
CLT D
When X1 , X2 , . . . are iid replications of X ∼ f , and `b is the sample mean based on N samples; then √ D N `b − ` → N(0, σ 2 ) D D Interpretation: `b ≈ ` + N(0, σ 2 /N) ∼ N(`, σ 2 /N); f`b(y)
`q = is q-th quantile of b distribution of `:
I
√ b ≈ σ/ N; First: the standard error of the estimator is SE[`]
I
Suppose you estimated ` by ` based on a simulated sample of size N;
I
You might consider ` as a random realisation of of the estimator `b ∼ f`b(y).
I
Thus (see graph previous slide), √ |` − `| < 1.96 σ/ N ⇔ 0.025 <
I
P(`b ≤ `q ) = q
D
P `b ≤ `
< 0.975
In other words, when you would generate many independent estimates (each based on sample N), about 95% of these would be within 1.96 standard errors of the true value.
Note: `0.025
` = `0.5
c Ad Ridder (VU)
`0.975
√ `0.975 ≈ `+1.96 σ/ N
y
Essay – January 2016
13 / 83
Random Estimation Interval
I
c Ad Ridder (VU)
Essay – January 2016
14 / 83
Confidence Interval
Generally, let α ∈ (0, 1) be your uncertainty parameter such as α = 0.05 or α = 0.1. 2
Var[Y] of the output
I
Let ` be your estimate, and suppose you know σ = variable.
I
Let zq be the q-th quantile of the standard normal distribution; i.e., Φ(zq ) = (N(0, 1) ≤ zq ) = q. We saw √ |`b − `| ≤ z1−α/2 σ/ N ≈ 1 − α | {z }
P
I
The analysis above is the basis for reporting simulation output;
I
Your estimate is ` computed as the average of N iid output data obtained by simulating N times the output variable Y = H(X);
I
You get a realisation of the random estimation interval: b ` + z1−α/2 SE[`] b I(`) = ` − z1−α/2 SE[`],
P
I
You report this interval as a 100(1 − α)% confidence interval.
I
Interpretation: when you would generate many independent estimates (each based on sample N), and when you would calculate the associated confidence intervals, then about 100(1 − α)% of these confidence intervals would cover the true value.
b SE[`] I
Rewrite: the unknown parameter ` satisfies b `b + z1−α/2 SE[`]) b P ` ∈ (|`b− z1−α/2 SE[`], {z } random estimation interval
c Ad Ridder (VU)
Essay – January 2016
≈1−α
b I(`)
15 / 83
c Ad Ridder (VU)
Essay – January 2016
16 / 83
Unknown Standard Error
Properties
I
Almost always the standard error of the estimator is unknown;
I
Equivalently, the variance σ 2 = variable is unknown;
I
Thus variance can be estimated by the sample variance
`b − ` D `b − ` D √ ∼ N(0, 1) ⇒ √ ∼ tN−1 , σ/ N S/ N
N N 1 X 2 H (Xi ) − `b2 . N − 1 N i=1
the (Student) t distribution with N − 1 degrees of freedom. (iv). Asymptotic normality: √ D N `b − ` /S → N(0, 1);
b by its estimate Hence, replace in your computations SE[`] √ c = S/ N. SE
c Ad Ridder (VU)
Essay – January 2016
17 / 83
Normal vs Student t
I
Essay – January 2016
18 / 83
E
Summary. Suppose we wish to compute a performance measure ` = [Y] where the output variable Y = H(X) can be calculated as a function H of a sequence of random input variables X. Suppose σ 2 = ar[Y] also unknown.
According to property (iii), we should use the Student t distribution in stead of the normal distribution when constructing confidence intervals: √ √ `b − tN−1,1−α/2 S/ N, `b + tN−1,1−α/2 S/ N ,
V
Monte Carlo algorithm 1. Repeat for i = 1, . . . , N: (i). Generate Xi independently of previous runs; (ii). Compute output Yi = H(Xi ).
According to property (iv), there is not much difference between the critical points for large sample size N: α t9,1−α/2 t49,1−α/2 t99,1−α/2 t249,1−α/2 t999,1−α/2 z1−α/2
c Ad Ridder (VU)
c Ad Ridder (VU)
The Monte-Carlo Algorithm
where tN−1,q is the q-th quantile for the t-distribution with N − 1 degrees of freedom. I
∀ > 0
(iii).
Mathematically equivalent: S2 =
P(|S2 − σ2 | > ) = 0
0.15 1.574 1.462 1.451 1.443 1.441 1.440
0.1 1.833 1.677 1.660 1.651 1.646 1.645
0.05 2.262 2.010 1.984 1.970 1.962 1.960
Essay – January 2016
2. Compute sample mean estimator for `:
0.025 2.685 2.312 2.276 2.255 2.245 2.241
N 1 X `b = Yi ; N i=1
3. Compute sample variance estimator for σ 2 :
S2 =
N 1 X b 2; (Yi − `) N − 1 i=1
4. Report estimate, confidence interval, and/or (estimated) standard error.
19 / 83
c Ad Ridder (VU)
Essay – January 2016
20 / 83
What is DES?
Discrete-event simulation is a computer simulation model/program of a stochastic system that evolves dynamically in time via state information which changes at discrete time epochs. The components of a DES are
Discrete Event Simulation
c Ad Ridder (VU)
Essay – January 2016
21 / 83
A queueing network example
I
Entities
I
Time variable (or simulation clock)
I
System state variables
I
Events and Event list (or calendar)
I
Global variables
I
Statistics collectors (or counter variables)
c Ad Ridder (VU)
Essay – January 2016
22 / 83
Modeling aspects
As a reference model we consider a queueing network model.
c Ad Ridder (VU)
Essay – January 2016
23 / 83
I
The queueing model consists of two stations (‘queues’) Q1 and Q2 .
I
At Qi jobs arrive according to a Poisson(λi ) process.
I
There are three servers at Q1 and two servers at Q2 .
I
The service times of the servers at Qi have distribution function Gi (·).
I
After service completion at Q1 the jobs enter Q2 .
I
After service completion at Q2 a job leaves the system with probability p or re-enters (‘feeds back’ to) Q1 with probability 1 − p.
I
Both stations have infinite waiting spaces.
I
Waiting jobs are served in order of arrival (FCFS).
c Ad Ridder (VU)
Essay – January 2016
24 / 83
Performance measures
Transient vs Steady-State
Typically we are interested in I
waiting times or system times at the two stations;
I
waiting lines or system lines at two stations;
I
time spend in the system (sojourn time);
I
throughput (per station or from the whole system);
I
utilization;
I
...
We have to specify I
whether we wish to estimate these performance measures for a finite period; for instance the queueing system operates daily, opening empty at 8.00 hr in the morning, closing at 18.00 hr in the evening.
I
or whether we wish to estimate steady-state averages; then we assume that the system operates for an infinite time.
Let’s assume that
c Ad Ridder (VU)
Essay – January 2016
25 / 83
Entities
Essay – January 2016
26 / 83
System state
I
Actually, a system is defined to be a collection of entities, e.g., people or cars or machines, that act and interact together.
I
Without entities, nothing would happen.
I
Entities have attributes, usually given as data values.
I
In our example the entities are (attributes between brackets)
I
c Ad Ridder (VU)
I
I
I
The system state is a collection of variables necessary to describe a system at a particular time. The set of all states is denoted by X ; a specific state by x ∈ X ; and the (random) state at time t by X(t) ∈ X . In our example the state comprises
I
jobs or customers (arrival time);
I
the number of jobs (x1 , x2 ) present at the two stations;
I
servers (idle/busy).
I
two vectors a1 , a2 of the arrival times of these waiting jobs;
I
two vectors b1 , b2 specifying the status of the servers (idle/busy).
Most often in a simulation study, we do not bother too much with entities.
c Ad Ridder (VU)
Essay – January 2016
27 / 83
c Ad Ridder (VU)
Essay – January 2016
28 / 83
Events
I
I
I
Time or clock variable
An event is an instantaneous occurrence that may change the state or trigger a state transition. The set of all possible events is denoted by E ; a specific event by e; the events active in state x ∈ X by E(x) ⊂ E .
I
Events occur at some point in time.
I
For this we need a variable representing the current time of the simulation.
I
This is also called the simulation clock or system time that measures the elapsed simulation time.
I
Clearly we need to specify its unit size (second or minute or ...); then we set it to zero at the start of a simulation and update it every time an event occurs.
In our example events are I
the arrival of a new job (at Q1 or at Q2 );
I
a service completion at one of the five servers.
c Ad Ridder (VU)
Essay – January 2016
29 / 83
Event list or calendar
Essay – January 2016
30 / 83
A Set of Alarm Clocks
I
The calendar for the simulation is a list of the events that are currently scheduled to occur.
I
There is only one event list and it consists of the scheduled event times, sorted by the earliest scheduled time first.
I
The event list at time t is denoted by Lev (t).
I
In our example the event list comprises I
the next arrival time of a customer at Q1 ;
I
the next arrival time of a customer at Q2 ;
I
for any busy server the next time he is ready completing the job.
c Ad Ridder (VU)
c Ad Ridder (VU)
Essay – January 2016
I
Another view of the event list is that it is a collection of alarm clocks, one for each scheduled event. The clocks are preset at different (random) times in the future. For instance at some t there are 3 events scheduled:
(
)
Lev (t) = arrival Q1
31 / 83
c Ad Ridder (VU)
arrival Q2
service completion Q2
Essay – January 2016
32 / 83
The Event-Scheduling Approach
A walk through DES In our queueing example we let time be measured in seconds, starting at 08:00. This is system time 0.
I
The discrete-event simulation runs as follows.
I
Suppose current time is tsim and state x ∈ X .
I
Each active event ei ∈ E(x) has an associated scheduled alarm time ti .