Tilburg University
Validation in simulation: bootstrap tests Jack P.C. Kleijnen Tilburg University SACCO Gregynoc (Wales), April 10-14, 2000
Overview l How
to validate in simulation?
Real system versus simulation model l Simulation model versus ‘metamodel’ (response surface; polynomial, NN, Kriging) l
l No
Gaussian simulation output: bootstrap But: no unique bootstrap solution! l Monte Carlo estimates of Type I & II errors 25-mei-2000
Kleijnen: WSC 2000
2
Trace-driven simulation Real arrival times a(i, t) = A(i); t = 1, …, k; i = 1, …, n
Simulated service times: random numbers of rep. s Ri(r) (r = 1, …, s)
Real system
Simulation model Sim. output vi; t(r)
Real output series wi; t Real performance Xi
Sim. performance Yi(r) Compare Xi and Yi(r)
Efron & Tibshirani ’s general bootstrap method l
l
l l l l
Observations Zi are i.i.d. with i = 1, ..., n (any distribution) Select any statistic T(Z1, ..., Zi , …, Zn) that is a continuous function of Zi Resample zi with replacement: Zi* (weights 0, 1, …, n) Compute bootstrap statistic T* =T(Z1* , ..., Zi*, …, Zn* ) Repeat b times: Tj* (j = 1, …, b) Sort Tj* to estimate 5% & 95% quantiles: 90% confidence interval for T 25-mei-2000
Kleijnen: WSC 2000
4
Six validation statistics T1 =
F (2, n - 2) (means & variances X, Y)
T2 = T3 =
∑ ∑
T4 =
∑
T5 =
∑
T6 = 25-mei-2000
∫
∞
−∞
n 1
n 1 n 1
i
/ n
with
D
i
= X
i
− Yi
2
n 1
D
Di / n
Di / n (Y i / X i ) / n
Fˆ x ( z ) − Gˆ y ( z ) dz
Kleijnen: WSC 2000
5
Bootstrap if s = 2 Sim. replicate #1 is valid model of #2! Replicate 1: Y(1, 1) … Y(1, i) … Y(1, n) Replicate 2: Y(2, 1) … Y(2, i) … Y(2, n)
Bootstrap technique: l l
25-mei-2000
Resample n pairs {Y(1, i), Y(2. i)} Compute
T = * 4
∑ {Y n
*
1
* 4;
T
*
(1, i ) - Y ( 2 , i )} / n
l
Repeat b times: sort
l
Reject H0 if original T4(X, Y1) or T4(X, Y2) falls outside c.i. - use Bonferroni (α /2)
find 10% and 90% quantiles: c.i.
Kleijnen: WSC 2000
6
Bootstrap if s > 2 l
Condition on trace Ai with i = 1, …, n Resample 1 pair {Y(r, i), Y(r’, i)} from column i with r, r’ = 1, …, s; r ≠ r’ Repeat for each of n columns
l
Compute
l
Reject H0 if any original T4(X, Yi) falls outside c.i. with
T4* ; etc.
Bonferroni (α /s)
25-mei-2000
Kleijnen: WSC 2000
7
Monte Carlo estimate of type I error n = 10; s = 10; " = 0.10 k = 1,000 (Gaussian?) 0.100
k = 10
0.095 0.090
0.05 0.04
0.085
0.03
0.080
0.02
0.075
0.01
0.070
0
T1
25-mei-2000
T4
T5
T6
Kleijnen: WSC 2000
T1
T4
T5
T6
8
Monte Carlo estimate of power
0.45
k = 10
0.4 0.35 0.3
T1 T4 T5 T6
0.25 0.2 0.15 0.1 0.05 0 0.8
0.9 25-mei-2000
1
1.1
1.2
1.3
1.4
ρ~
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
T4
k =10 k = 1000
0.8
Kleijnen: WSC 2000
0.9
1.0
1.1
1.2
1.3
ρ~
1.4 9
Conclusions of Part 1 l
l
Bootstrap in general • Any distribution • Any continuous statistic • Non-unique technique: art Bootstrap & validation of trace-driven simulation • 6 statistics: T1 through T6 • Number of replicates: s = 1, 2, >2 s > 2: T4 (average X versus average Y) • Advice: replicate more than 2 times
25-mei-2000
Kleijnen: WSC 2000
10
Reminder: Reality, simulation, metamodel Real system Validation Simulation model Validation Metamodel 25-mei-2000
Kleijnen: WSC 2000
11
Part 2: Simulation model versus ‘metamodel’ Metamodel: y = X$ + e with x(i) = (1, x(i, 2), …, x(i, q - 1)) (i = 1, …, n) −1 ˆ β = ( X ' X ) X 'w l OLS estimator: with w(i, r) simulation output of replicate r (r = 1, …, s) ~ l EGLS: β non-linear g(X, w) l
25-mei-2000
Kleijnen: WSC 2000
12
Validation statistics l
Rao (‘59): Estimate residuals (fitting errors)
~ ~ e i = wi − yi
l l
with
~ ~ yi = xi′ β
Expected value 0 iff valid metamodel Rao’s statistic R( e~ ) = g[w(i, r), x(i)] R = F(n - q, s - n + q) iff valid metamodel & Gaussian simulation output w
25-mei-2000
Kleijnen: WSC 2000
13
Bootstrap of Rao’s statistic Z(i) becomes simulation I/O: x(i), w(i, r) with i = 1, …, n; r = 1, …, s l H0: Valid metamodel: E(w) = X$ l
l
Wrong bootstrap methods (H1, not H0): • Resample w(i, r) (estimate $, e, R) ~ e • Resample , etc.
25-mei-2000
Kleijnen: WSC 2000
14
Correct bootstrap method w E(w) = x
x
Simulate (n = 3; s = 2) w E(w) = x Seed 2 Seed 1
x1
x2
x3
x
Fit regression metamodel w E(w) = x
~ ~ β 0 + β1 x
x1
x2
x3
x
Compute deviation d(i, r) = w(i, r) - w (i ) w E(w) = x
~ ~ β 0 + β1 x
x1
x2
x3
x
Summary: Correct bootstrap of Rao’s R under H0 l l l l l
Compute deviation d(i, r) = w(i, r) - w (i ) Resample d(i, r): d
* i ;r
~ * Compute w = x i β + d i ;r ~* ~* * From (X, w*) compute β , e , R * i ;r
* ( 0.8 b)
Repeat b times; sort R , find c.i.; reject H0 if R > R Note: simulation outputs are correlated with non-constant variances: resample m times the n-dimensional vectors d*
25-mei-2000
*
Kleijnen: WSC 2000
19
Four examples; s = 53; " = .20; b = 1,000; m =100 l
l
l
l
Case 1: H0: E(w) = x; Gaussian w Use F-table for Rao’s R: αˆ = 0.19 Bootstrap: = 0.19 Case 2: H1: E(w) = x + x2/4; Gaussian w Use F-table for Rao’s R: = 0.57 Bootstrap: = 0.58 Case 3: H0; lognormal w Use F-table for Rao’s R: = 0.27 Bootstrap: = 0.24 Case 4: H1; lognormal w F-table for Rao’s R: = 0.68 Bootstrap: = 0.63
25-mei-2000
Kleijnen: WSC 2000
20
Conclusions of Part 2 l l
l l
l
Bootstrap is an art: resample deviations! Bootstrap versus F-table for Rao’s statistic: • Gaussian output: same type I & II errors • Lognormal output: bootstrap gives better type I error Explore other cases: s = 5, etc. Explore other statistics • Cross-validation statistics (PRESS, DEFITS, DFBETAS, CoOk’s D, RTSTUDENT) • Relative absolute fitting errors References: http://cwis.kub.nl/~few5/center/staff/kleijnen/
25-mei-2000
Kleijnen: WSC 2000
21
