Statistische Analyse von Simulationsergebnissen

• Problemstellung:
  • Modell mit stochastischen Prozessen
    • → Simulationsergebnisse sind zufällig
    • → Ergebnisse müssen statistisch analysiert werden
  • Standardtests oft nicht anwendbar
    • Ergebnisse zeitlich korreliert
    • Verteilungen nicht stationär
    • Nullhypothese "Modelldaten ≙ Systemdaten" ist sicher falsch (bei hinreichender Genauigkeit)
• Standardbeispiel M|G|1-Queuing-System:
  • Ankunftszeiten ~ Ex(λ)
  • Servicezeiten ~ N(1/μ, σ²)
  • Ausgaben
    • aktuelle Queuelänge nQ(t_i) mit Zeiten t_i
    • mittlere Queuelänge L_Q
    • mittlere Wartezeit in der Queue W_Q
    • mittlere Serverauslastung A
  • Beispiellauf mit Parametern λ = 1, μ = 1.2, σ = 0.9, t_End = 10000
    • 
  • Ergebnisse von 6 Beispielläufen mit aufeinanderfolgenden Seeds

    • Lauf	LQ	WQ	A
      1	8.0348	8.0814	0.9201
      2	9.1228	9.0291	0.9287
      3	7.3117	7.3779	0.9129
      4	5.6341	5.8311	0.8785
      5	10.2697	10.2059	0.9338
      6	15.8034	15.5841	0.9428

• Abschätzen der mittleren Queuelänge - naiv (und falsch):
  • betrachten Queuelänge X_i zu Zeiten t_i = i (i = 1 .. n)
    • aus Modellergebnissen interpoliert (mit 'previous')
  • mit üblichen Standardformeln aus den X_i berechnet
    • Mittelwert L_Q = 8.0380
    • Varianz S² = 61.5559
    • zum Nachrechnen: t_End = 10000, seed = 4
  • Problem: L_Q erscheint zu klein (angesichts der Tabelle)
  • Ursache: Korrelation der Werte innerhalb eines Laufs
  • genauer: Korrelationskoeffizient ρ jeweils aufeinanderfolgender Werte X_i und X_i+1 abschätzen
    • 
  • Ergebnis
    • ρ = 0.9859
    • Werte sehr stark korreliert
    • klar, z.B. Queue voll → aufeinander folgende X_i alle groß
• Ensemblemittelung:
  •  k komplette Läufe mit verschiedenen Seeds
    • → Daten X_ij, i = 1..k, j = 1..n
  • daraus Mittelwerte Y_i pro Lauf bestimmen
    • 
  • Gesamtmittel = Mittel über Läufe (Ensemblemittel)
    • 
  • Y_i sind i.i.d → Standardverfahren funktionieren
    • Schätzer für Varianz
    • 
    • Schätzer für (1-α)-Konfidenzintervall
    • 
• Ensemblemittelung im Beispiel:
  • k = 100 Läufe bis t_End = 100 (statt t_End = 10000)
  • Verteilung der Y_i und der Standardabweichungen pro Lauf
    • 
  • Ergebnisse
    • = 3.8414
    • S² = 9.5769
    • Konfidenzintervall: [3.2274, 4.4554] für 1-α = 95% 
  • Ergebnis seltsam (vgl. Tabelle)
    • Mittelwert anscheinend viel zu klein
    • Konfidenzintervall enthält keins von 6 Ergebnissen
  • Ursache
    • Läufe zu kurz
    • größere Queuelängen erst nach längerer Zeit
    • Gleichgewicht (falls vorhanden) noch nicht erreicht
  • Abhilfe
    • "Aufwärmzeit" abwarten, erst danach Werte nehmen
    • Problem: Wie lange dauert die?
• Ermitteln der Aufwärmzeit t_A:
  • grundsätzlich
    • häufig: Verteilungen für X(t_i) haben Limes für t→∞ (Steady State)
    • gesucht: Laufzeit t, ab der Steady State ungefähr erreicht ist
    • schwierig, da X_i auch im Steady-State stark schwanken
  • graphisches Verfahren nach Welch [12]
  • 1. mache k Läufe der Länge n → Daten X_ij, i = 1..k, j = 1..n
    • z.B. k = 10, n groß
  • 2. mittle Werte gleicher Zeit über die Läufe
    • 
    • verringert Varianz um k bei gleichem Erwartungswert
  • 3. mittle Werte über Zeitfenster der Breite w
    • 
    • glättet hochfrequente Anteile weg
    • plotte _j(w) für verschiedene w
    • wähle w so klein, dass Kurven "hinreichend glatt"
    • im Zweifelsfall k vergrößern
  • 4. plotte _j(w) und wähle t_A, ab dem Grenzwert etwa erreicht
• Vorgehensweise im Beispiel (Replication/Deletion-Verfahren):
  • Aufwärmzeit t_A nach Welsh ermitteln
    • k = 10 Läufe bis t_End = 1000
    • 
    • geglättete Kurven mit w = 250 und w = 500 sehr ähnlich
    • Schwankungen bleiben groß
    • Aufwärmbereich grob gut erkennbar
    • wähle t_A = 500 (eher etwas zu groß)
  • Achtung Matlab
    • gleitendes Mittel mit movmean benutzt am Rand einseitige Intervalle
    • damit wird der (entscheidende) Anstieg am Anfang weggemittelt
    • also: eigene Routine modmean2, die korrekte Formel verwendet
  • anschließend mit neuen Läufen X_ij bestimmen (Replication)
    • jeweils Werte vor t_A streichen (Deletion)
    • mit den restlichen Werten Mittel Y_i über Läufe bestimmen
    • Y_i mit Standardverfahren auswerten
  • Ergebnisse bei k = 10 und t_End = 1500
    • = 10.7314
    • S² = 56.1750
    • Konfidenzintervall: [5.3698, 16.0930]
  • gut verträglich mit Tabelle
  • Ergebnisse bei k = 100 und tEnd = 1500
    • =  9.4257
    • S² = 29.7782
    • Konfidenzintervall: [8.3430, 10.5085]
    • mehr Daten → deutlich höhere Genauigkeit
• Weitergehende Methoden:
  • bisher nur an der Oberfläche gekratzt
  • verschiedene Verfahren für endlichen (festen) Zeitrahmen, Steady-State-Systeme oder steady-state-zyklische Systeme
  • komplexe Methoden zum Vergleich zweier Systemkonfigurationen
  • ausgefeilte Verfahren zur Veringerung der Varianz
  • guter Einstieg: [Law, Kap 9-11]
• Aufgaben:
  • Aufgabe 20
  • Aufgabe 21