Folie 6

Simulation seltener Ereignisse

• Zufällige Eingaben in stochastischen Simulationen erzeugen zufällige Ausgaben, was eine statistische Auswertung erforderlich macht.
• Ergebnisse werden in Form von Mittelwertschätzern und Konfidenzintervallen zu einem gegebenen Konfidenzniveau α dargestellt.
• Die Simulation komplexer Systeme erfordert hohen Aufwand, um genaue Ergebnisse durch viele Replikationen oder lange Beobachtungsreihen. zu erzielen.
• Zeit- und Rechnerkapazitätsbeschränkungen verhindern oft genügend genaue Analysen.
• Techniken zur Erhöhung der Genauigkeit bei vorgegebenem Aufwand sind von großer praktischer Bedeutung.
• Die Berechnung des Konfidenzintervalls erfolgt über Ỹ ± ɛα · S/√n.
• P soll mit einer relativen Genauigkeit von 10% des ermittelten Mittelwertes zum Signifikanzniveau a=0.01 (99%) geschätzt werden.
• Wenn die Wahrscheinlichkeit des Ereignisses kleiner wird, sind mehr Beobachtungen erforderlich, um den Erwartungswert mit gleicher relativer Genauigkeit zu schätzen.
• Ereignisse mit kleiner Wahrscheinlichkeit sind in vielen Fällen von Interesse.
  • Beispielsweise soll die Ausfallwahrscheinlichkeit eines Systems in einem vorgegebenen Intervall <10-6 sein.
  • Die Wahrscheinlichkeit des Paketverlusts in einem Rechnernetz soll < 10-9 sein.
  • Es sind 664 Milliarden Beobachtungen notwendig, um eine Verlustwahrscheinlichkeit von mit 10% Genauigkeit (bei 99% Signifikanzniveau) zu schätzen.
  • Bei 10 Sekunden Dauer einer Replikation wären dies mehr als 200 000 Jahre
• Mögliche Beschleunigungsstrategien:
  • Schnellere Rechner nutzen (jedoch durch Technologie begrenzt und teuer).
  • Effizienteren Code schreiben (oft aufwändig und mit geringen Effekten).
  • Parallelisierung (mit realistischen Grenzen bei 100 bis 10000 Replikationen).
  • Reduktion der Varianz.
• Varianzreduktionstechniken (VRTs) verkleinern die Varianz des Schätzers, ohne den Mittelwert zu verändern, was die Breite der Konfidenzintervalle proportional reduziert.
• Es existieren zahlreiche VRTs, die jedoch normalerweise vom konkreten Modell abhängen und nicht blind anwendbar sind.
• Eine Abschätzung der Varianzreduktion ist oft vor der Simulation nicht möglich, und Pilotläufe sind zur Abschätzung des Effekts notwendig.
• Der Einsatz von VRTs erfordert zusätzlichen Aufwand, der in Beziehung zur erreichten Varianzreduktion gesetzt werden muss.
• Es existieren unterschiedliche Techniken der Varianzreduktion, wobei zwei Ansätze Verwendung finden:
  • Veränderung des beobachteten stochastischen Prozesses, so dass der Erwartungswert unverändert bleibt oder rekonstruiert werden kann, die Varianz aber kleiner wird.
  • Beobachtung des Modells primär in dem Bereich, in dem seltene Ereignisse auftreten und anschließendes Zurückrechnen des Erwartungswertes.
• Folgende Themen sind relevant:
  • Antithetische Variablen
  • Kontrollvariationen
  • Konditionierung
  • Importance Sampling

Antithetische Variablen

• Ziel ist die Verkleinerung der Varianz eines Schätzers bei der Analyse eines Modells durch Bestimmung von Konfidenzintervallen über Replikationen.
• Die Idee von Antithetischen Variablen basiert aus gepaarten Replikationen.
  • Kleine Beobachtungswerte im ersten Lauf entsprechen großen Beobachtungswerten im zweiten Lauf und umgekehrt (Einführung negativer Korrelation). -Mittelwerte der zwei werden im Durchschnitt näher am Erwartungswert liegen, als der Durchschnitt der Mittelwerte zweier unabhängiger Replikationen.
• Die Realisierung erfolgt durch Verwendung komplementärer Zufallszahlen (ZZs), d.h. u im ersten und 1-u im zweiten Lauf.
  • Es ist sicherzustellen, dass u und 1-u für identische Ereignisse verwendet werden (Definition von ZZ-Strömen für Ereignisse).

Mathematische Grundlagen

• Für zwei beobachtete Sequenzen v1 und v2, die durch Verwendung antithetischer Variablen entstehen, gelten folgende Annahmen:
  • Die Werte von vij und vkl sind unabhängig, falls j≠l (auch wenn i=k).
  • μ = E[vij] = E[wj] mit wj = (v1j + v2j)/2
  • vij sind Realisierungen der Zufallsvariable (ZV) Vi, wj sind Realisierungen der ZV W.
• Bei unabhängigen Beobachtungen gilt COV(V1,V2)=0.
• Falls V1 und V2 negativ korreliert sind, gilt COV(V1, V2) < 0, was zu einer kleineren Varianz des Mittelwertschätzers führt.
• Es ist nicht immer klar, ob antithetische Variablen wirklich zu einer Varianzreduktion führen, da ein monotones Modellverhalten in den antithetischen Variablen vorausgesetzt wird.
• Die Anwendung von antithetischen Variablen muss durch Probeläufe oder Interpretation des Systemverhaltens validiert werden, um sicherzustellen, dass eine Varianzreduktion eintritt.

Kontrollvariationen

• Die Ausnutzung der Korrelation zwischen verschiedenen Größen zur Reduktion der Varianz des Schätzers der gesuchten Größe wird verwendet.
• Der Erwartungswert µ der ZV V soll aus der Stichprobe (v1, ..., vn) geschätzt werden, wobei angenommen wird, dass V positiv oder negativ mit einer ZV X korreliert ist und der Erwartungswert v von X bekannt ist.
• Beispiele:
  • In einem einfachen Bediensystem steigt die Verweilzeit, wenn die Bedienzeiten steigen (positive Korrelation).
  • In einem einfachen Bediensystem sinkt die Verweilzeit, wenn die Ankunftsabstände größer werden (negative Korrelation).
  • In einem technischen System steigt die Systemverfügbarkeit, wenn die Zuverlässigkeit der Komponenten steigt.
• Einsatz des Ansatzes in der Simulation:
  • Identifiziere Resultatgrößen und interne Parameter, die positiv korreliert sind (z.B. Bedienzeit (Erwartungswert bekannt) und Verweilzeit (Erwartungswert zu schätzen)).
  • Simuliere und messe bzw. speichere die Realisierungen beider Größen.
  • Werte die Ergebnisse unter Nutzung des Zusammenhangs zwischen den beiden Größen aus.
  • Notwendige Daten sind in allen Simulatoren einfach zu ermitteln.

Mathematische Grundlagen

• Definition der „kontrollierten“ ZV: Z = V - a (E[X] – v).
• Es gilt: E[Z] = E[V] = μ und σ²(Z) = σ²(V) + a² · σ²(X) – 2a · COV(V,X)
• Die Wahl von "a" erfolgt so, dass:
  • a > 0, falls V und X positiv korreliert sind
  • a < 0, falls V und X negativ korreliert sind
• Die Varianz von Z ist kleiner als die Varianz von V, falls 2a · COV(V,X) > a² · σ²(X).
• Das Optimum wird erreicht, falls die erste Ableitung 0 wird: 2a · σ²(X) – 2 · COV(V,X) = 0
• Daraus ergibt sich die Lösung: aopt = COV(V,X) / σ²(X).
• Falls a = aopt, gilt: σ²(Z) = σ²(V) - (COV(V,X))²/σ²(X) = (1 – ρ²VX)·σ²(V).
• Dabei ist 0 ≤ ρ²VX = (COV(V,X))²/σ²(V)σ²(X) ≤ 1, was bedeutet, dass die Varianz nicht größer wird.
• Probleme bei der praktischen Umsetzung:
  • σ²(X) könnte unbekannt sein.
  • COV(V,X) ist in fast allen Fällen unbekannt.
  • aopt kann daher nicht exakt ermittelt werden und muss aus Simulationsergebnissen geschätzt werden.
• Mögliche Schätzung von aopt:
  • Seien (x1, ..., xn) Realisierungen von X und (v1, ..., vn) Realisierungen von V in der Simulation, und V und X die Mittelwertschätzer von V und X.
  • Schätzer für die Kovarianz: COV(V, X) = ∑(vj - V)(xj - X) / (n - 1).
  • Schätzer für die Varianz: S²(X) = ∑(xj - X)² / (n - 1), womit gilt aopt = COV(V, X) / S²(X).
  • Schätzer der gesuchten Größe: Z = V – aopt · (X – v)
• Problem: Die einzelnen Schätzer sind abhängig, wenn sie aus den Daten eines Simulationslaufs resultieren.
  • Es ist nicht klar, ob Z ein erwartungstreuer Schätzer ist.
  • Es ist nicht klar, ob die Varianz wirklich reduziert wird.

Konditionierung

• Die Beobachtung der „interessierenden“ Zufallsvariable V unter der Bedingung, dass eine andere Zufallsvariable X im Modell einen bestimmten Wert annimmt.
  • E[V|X=x] ist der Erwartungswert von V, wenn X den Wert x hat.
• Für eine diskrete Zufallsvariable X gilt: EX [E [V|X]] = ∑ E [V|X = xi] · p(xi), wobei p(x) die Wahrscheinlichkeit ist, dass der Wert von X gleich x ist (normalerweise unbekannt).
• Für die Varianz gilt: σ²X(E[V|X]) = σ²(X) – EX[σ²(V|X)] ≤ σ²(X). - Das Vorgehen beinhaltet:
  • Simuliere Realisierungen von X und für jedes beobachtete x simuliere Realisierungen von V|X=x.
• Das Vorgehen bringt Vorteile, falls:
  • Realisierungen von X effizient realisiert werden können.
  • E[V|X] analytisch berechnet werden kann.
  • EX[σ²(V|X)] groß ist.
• Die Methode ist stark modellabhängig.
• Typisches Beispiel: Modell eines Computersystems.
  • Ziel: Bestimmung der mittleren Verweilzeit an CPU, Platte und CD.
  • Problem: Nur 8% der Jobs greifen in einem Zyklus auf die CD zu → Simulation liefert wenige Beobachtungen der Zugriffszeit → große Varianz des Schätzers.
• Mögliche Konditionierung:
  • Annahme, jeder Job, der die CPU verlässt, würde zur CD wechseln.
  • Berechnung der „potenziellen“ Verweilzeit (analytisch, falls Bedienzeit exponentiell, auf Basis des Systemzustandes sonst).
  • Falls sich mehrere Jobs in der CD-Warteschlange befinden, wird der Systemzustand festgehalten und mehrere Replikationen gestartet.
• Weiteres Beispiel: Bestimmung der Überlaufwahrscheinlichkeit eines Puffers in einem Kommunikationssystem.
  • Falls K groß ist und die Auslastung moderat ist, ist die gesuchte Wahrscheinlichkeit sehr klein.
  • Ihre Bestimmung ist trotz großer praktischer Bedeutung (z.B. Dimensionierung von Puffern zur Vermeidung von Datenverlust mit typischen Anforderungen von 10-8 – 10-9 für die Verlustwahrscheinlichkeit).
• Vorgehen bei der Simulation zur Bestimmung der Wahrscheinlichkeit:
  • Definiere Regenerationszustand (d.h. Zustand, bei dem das zukünftige Verhalten nicht von der Vergangenheit abhängt, z.B. leeres System).
  • Simuliere R Zyklen, die jeweils im Regenerationszustand starten und enden.
  • Sei br (= Anzahl verlorene Pakete/Anzahl Pakete) Schätzer für die Verlustwahrscheinlichkeit im Zyklus r (b sind unabhängig identisch verteilt).
• Problem: Für die meisten Zyklen ist b₁=0.
• Idee eines Konditionierungsansatzes
  • Wähle L < K, simuliere bis Pufferfüllung L erreicht und speichere Zustand.
  • Beginne vom gespeicherten Zustand mit r Replikationen, die jeweils enden, wenn Warteschlangenlänge L wieder erreicht wird.
  • Nach der r-ten Replikation fahre mit der Simulation fort, bis die Replikation endet.
• Aus der Simulation kann geschätzt werden:
  • P(n=L) Wahrscheinlichkeit, dass Population L in einer Replikation erreicht wird.
  • P(n=K|no=L) Wahrscheinlichkeit, dass Population K erreicht wird, falls mit Population L begonnen wird
  • P(n=K) = P(n=K|no=L)·P(n=L)
• Beispiel: M/M/1/3-System (analytisch analysierbar)
  • Untersucht werden: Ankunftsrate λ, Bedienrate μ und Übergangswahrscheinlichkeit p = λ /(λ + μ)
  • Systemverhalten kann durch folgende Markov-Kette beschrieben werden.
  • Analyseziel: Bestimmung der Wahrscheinlichkeit, dass ausgehend von Zustand 1, Zustand 3 erreicht wurde, ohne dass Zustand 0 oder 1 vorher erreicht wird.
  • X ZV mit X=1 falls 3 erreicht wird, 0 sonst, E(X) ist der gesuchte Wert
  • Vorgehen bei der Simulation: Simuliere R Replikationen, die in Zustand 1 starten, bestimme xr = Wert von X in Replikation r, berechne X und S²(X), und das Konfidenzintervall.
• Es gilt:
  • E[X]=E[X]=p² und Varianz σ²(X) = p – p4 = p· (1 – p²).
  • Für p < 1 gilt σ²(X) ≈ p ⇒ σ²(X) ≈ p / R.
  • Vereinfachung durch Konditionierung möglich.
  • Deutliche Aufwandsreduktion.
• Natürliche Verallgemeinerung: Mehrere Punkte, an denen neue Replikationen ansetzen, unterschiedliche Punkte für ab- und aufsteigende Kurven, dynamische Bestimmung der jeweiligen Replikationszahlen
• Probleme bei der konkreten Anwendung:
  • Festlegung der Punkte, wo neue Replikationen ansetzen, Festlegung der Zahl der Replikationen und Integration in Modellierungswerkzeuge.

Importance Sampling

• Sei f(x) Dichtefunktion einer ZV X, die den Ausgang eines Simulationsexperiments beschreibt und g(x) eine Funktion.
• Ermittelt werden soll γ = Ef[g(X)] = ∫ g(x) · f(x)dx.
• Falls P(y < X) = ∫ f(x)dx < 1 ist die Ermittlung von E[g(X)] schwierig, da die meisten Experimente x < y und damit g(x)=0 liefern.
• Sei X' eine andere ZV mit (beliebiger) Dichtefunktion f'(x), so dass P(y < X') = ∫ f'(x)dx > 0 und f(x) > 0 ⇒ f'(x) > 0.
• Dann gilt γ = Ef[g(X)] = ∫ g(x) · [f(x) / f'(x)] · f'(x)dx = Ef'[g(x)L(x)], wobeiL(x) = f(x) / f'(x) (Likelihood-Funktion)
• Voraussetzung für die Verbesserung der Schätzung :
  • f'(x) muss im „interessanten“ Bereich groß sein !
  • Viele Beobachtungen führen zu relevanten Resultaten (X>y).
• Sei zur Vereinfachung g(x) = 1 falls x > y und 0 sonst
• Vorgehen in der Simulation zur Bestimmung von γ:
  • Erzeuge R Stichproben x1,...,xR gemäß Dfkt. f(x).
  • Berechne γ = 1/R Σ δ(xr > y) mit offensichtlich Ef[γ] = γ.
  • Konfidenzintervall: γ ± εα·√(γ·(1-γ)/R) .
  • Für eine Konfidenzintervallbreite, die proportional zu y sein soll, sind R∈O(1/γ) Stichproben notwendig.
  • Simulationsaufwand wächst proportional zum Kehrwert von y über alle Grenzen (unbounded relative error).
• Vorgehen beim Einsatz von „importance sampling“ zur Bestimmung von γ:
  • Erzeuge R Stichproben y1,...,yR gemäß Dfkt. f'(x) und bestimme dabei jeweils die Werte der Likelihood-Funktion L(xr) = f(xr)/f'(xr).
  • Berechne γ = 1/R Σ δ(yr > y) · L(yr) (offensichtlich Ef[γ] = γ).
  • Da f' mit höherer Wahrscheinlichkeit Werte größer y liefert, werden mehr Beobachtungswerte ≠ 0 sein, dies wird durch L(xr) ausgeglichen !

Mathematik Hintergründe

• Da (δ(x>y))² = δ(x>y), gilt für das zweite Moment:
  • E[Z²] = Ef' [δ(Y > y) · L(Y)]² = ∫ δ(x > y) · (f(x)/f'(x))² · f'(x)dx = ∫ δ(x > y) · [f(x)/f'(x)] · f(x)dx = Ef' [δ(x > y) · L(x)]
  • Zur Reduktion der Varianz muss f(x)/f'(x) klein sein, wenn x > y.
  • Da f(x) für x > y klein ist (seltenes Ereignis), sollte f'(x) in diesem Fall groß sein
  • Angestrebt wird , dass die Größe der Stichprobe, die zum Erzielen einer bestimmten relativen Genauigkeit notwendig ist, unabhängig von der Größe von y ist (bounded relative error).
• Formale Formulierung
  • Sei ε ein Seltenheitsparameter (d.h. limε→∞ γ(ε) = 0)
  • γ(ε ) ist mit vorgegebener relativer Genauigkeit α mit O(R) Replikationen für alle ε > 0 bestimmbar
• Für allgemeine Systeme vom Typ GI/GI/1 können asymptotisch optimale Werte aus den momentgenerierenden Funktionen der Zwischenankunfts- und Bedienzeitverteilung gewonnen werden.
• Für M/M/1-Systeme resultiert das asymptotisch optimale Modell aus der Vertauschung der Ankunfts- und Bedienrate (System wird instabil !)
• Mehrbedienersystem GI/GI/N lassen sich ähnlich wie GI/GI/1 behandeln
• Kaum Resultate für Warteschlangennetze (außer Tandemnetzen) oder System mit korrelierten Ankunftsströmen
• Heuristiken funktionieren nur: für die angegebenen Systeme für das angegebene Leistungsmaß für große K und potenziell unendliche Zeitintervalle
• Modelle können aus mehreren Komponenten bestehen, die ausfallen können und repariert werden. In der Regel sind Ausfallraten wesentlich kleiner als Reparaturraten .
• Fehlerwahrscheinlichkeit wird gleichmäßig für alle Fehler erhöht (Fehlerwahrscheinlichkeit = Wahrscheinlichkeit, dass innerhalb eines vorgegebenen Intervalls mindestens ein Ausfall stattfindet) Verfahren funktioniert gut, solange alle Fehler ungefähr gleich wahrscheinlich sind und ähnliche Effekte bzgl. des Systemverhaltens haben.
• Falls aber unterschiedliche Fehlerraten vorliegen, kann es zu Verfälschungen im Systemverhalten kommen, so dass Schätzer verzerrt sind.
• Insgesamt bietet importance sampling theoretisch die Möglichkeit mit festem Aufwand beliebig kleine Wahrscheinlichkeiten zu schätzen.
  • Andererseits gibt es in der Praxis ist aber eine optimale (oder auch nur relativ gute) Realisierung nur für sehr einfache Beispiele und Maße
  • Es gibt Modelle, bei denen nachweislich keine statische Parameteränderung zu guten Ergebnissen führt (bounded relative error). Fehlt der Durchbruch zur praktischen Anwendbarkeit.
  • Ist eine automatische Anwendung auf komplexere Modelle und damit auch eine Integration in Simulationswerkzeuge noch nicht absehbar