Bastian Bluhm Informationswirt 5.Semester 1236026

11.07.07

Seminarthema an der Universität Karlsruhe(TH)

„Einsatz von Varianzreduktionstechniken II : Stratified Sampling und Common Random Numbers“ im Rahmen des Seminars

„Ausgewählte technische,rechtliche und ökonomische Aspekte des Entwurfs von Fahrerassistenzsystemen(SS 07)“

betreut von Christiane Barz

Gliederung:

Teil A: Entwurf der 4-armigen Kreuzung: 1. Die Piosson-Verteilung 2. Erzeugen eines homogen Poisson-Prozesses 3. Erstellen der Kreuzung, sowie der ankommenden Autos, rechnerunterstützt mit Hilfe der Programmiersprache Java 4. Berechnung der durchschnittlichen Warteschlangenlänge, sowie der durchschnittlichen Verweildauer im System. 5. Berechnung der Standardabweichung

Teil B: Varianzreduktionstechniken 1. Common Random Numbers 2. Stratified Sampling a) Die Idee b) Beispiel anhand der Gleichverteilung c) Die Varianzreduktion mit Stratified Sampling d) Beispiel an einem Poissonprozess 3. Fazit und Quellenangabe

Teil A: Entwurf der 4-armigen Kreuzung 1. Die Piosson-Verteilung : Ist eine diskrete Wahrscheinlichkeitsverteilung, die beim mehrmaligen Durchführen eines Bernoulli-Experiment entsteht. Ein Bernoulli-Experiment kennt nur zwei mögliche Ergebnisse (z.B. „Erfolg“ und „Misserfolg“). Die Wahrscheinlichkeit für ein solches Ereignis soll möglichst klein sein. Die Poisson-Verteilung wird deshalb auch als die Verteilung der seltenen Ereignisse bezeichnet. Die mit Pλ bezeichnete Verteilungsfunktion wird durch den Ereignisrate genannten Parameter λ bestimmt, der gleichzeitig Erwartungswert und Varianz der Verteilung ist. Die Verteilungsfunktion lautet : mit k=0,1,2.. Darstellung einer Poissonverteilung mit λ = 6.

2. Erzeugen eines homogen Poisson-Prozesses: Der Poisson-Prozess ist ein Erneuerungsprozess,dessen Zuwächse (I)unabhängig und (II)Poisson-Verteilt mit Parameter λ sind. Nur wenn I und II gilt ist der PoissonProzess homogen. Pro Zeiteinheit werden genau λ Sprünge erwartet, λ bezeichnet man auch als die Intensität des Prozesses. Die Wartezeiten zwischen den Sprüngen sindexponentialverteilt.

1. 2.

Die Ankunftszeiten sollen über den Zeitraum (0,t) verteilt sein. Um dies zu erreichen wird folgender Algorithmus verwendet: Setze n=0 , t 0 =0 Erzeuge Zufallszahl z Gleichverteilt [0,1] (abgekürzt ZZ)

3. 4.

Setze t  n1=t n −1/ ln z  Sobald tn+1>t stoppe. Sonst weiter mit 1. (I) n= n + 1; 2. (II) Gehe zu Punkt 2. Dies wurde von mir in der Methode Poisson umgesetzt. Über den Zeitraum (25,26) mit λ = 5 erhält man je nach ZZ poissonverteilte Ankunftszeiten der Form: R2 Nr 107: 25.118917733198128 R2 Nr 108: 25.183858105496636 R2 Nr 109: 25.5052747595677 R2 Nr 110: 25.71186247115331 R2 Nr 111: 25.973535714116682 Anhand der Zählvariable erkennt man das das System die Poisson-Verteilung von [0,T] erstellt und diese dann ab einem eingeschwungenen Zustand abgreift. Zur Erzeugung der ZZ sowie zur Berechnung des Logarithmus wurde die Java Klasse Math() verwendet.

3. Erstellen der Kreuzung, sowie der ankommenden Autos, rechnerunterstützt mit Hilfe der Programmiersprache Java. Die Kreuzung ist 4-armig und die sich gegenüberliegenden Richtungen (1 und 3 sowie 2 und 4 ) besitzen die selben λ zum erstellen des Poissonprozesses. Der Poissonprozess stellt die Ankunftszeiten dar, die Bedienzeit(Dauer zum Überqueren der Kreuzung) beträgt konstant 5 Sekunden. Der somit erzeugte Poissonprozess wird in 4 sogenannten LinkedList abgespeichert, welche folgende Funktionen haben LinkedList.size() gibt die Länge der Liste aus LinkedList.getFirst() und LinkedList.getLast(), sowie LinkedList.get(k) sind selbsterklärend.Statt „get“ kann auch „remove“ zum Löschen genutzt werden. Mit LinkedList.add(x) werden Elemente hinzugefügt. Jedes erzeugte Auto besitzt eine Fahrtrichtung, diese ist an unterschiedliche Wahrscheinlichkeiten für geradeaus fahren bzw. rechts oder links abbiegen gekoppelt. Richtung1 und 3 fährt mit 50% Whkt. geradeaus und mit jeWhkt. 25% rechts oder links. Richtung2 und 4 fährt mit 30% Whkt. geradeaus und mit jeWhkt. 35% rechts oder links. In welche Richtung ein Auto fährt wird zufällig mit Hilfe der Random Klasse ermittelt ( liefert GV ZZ zwischen [0,1] ). Abspeicherung analog zum Poissonprozess in LinkedListen. Die Vorfahrtsrechte der einzelnen Autos ergeben sich nach einer FIFO (first-in,first-out) Strategie. Das heisst, es darf immer das Auto als nächstes fahren das sich am längsten an der Kreuzung befindet.

Über ein verschachteltes System von If-Anweisungen wird ermittelt, wer sich als erster an der Kreuzung befindet, welche Autos sich auch noch an der Kreuzung befinden und welchem sich von diesen zusätzlichen Autos noch ein Fahrtrecht ergibt. Dies stellt die Fahrerassistenzkomponente dar,quasi eine grosse Steuereinheit die die Vorfahrt regelt, um die Wartezeit und Warteschlange möglichst zu minimieren.

4.Berechnung der durchschnittlichen Warteschlangenlänge, sowie der durchschnittlichen Verweildauer im System. Um die durchschnittliche Warteschlangenlänge zu berechnen, wird die Dauer die jedes Auto in der Warteschlange verbracht hat aufsummiert und durch die beobachtete Zeit geteilt. So kann die Warteschlangenlänge sekundengenau ausgegeben werden. Die durchschnittliche Verweildauer im System berechnet sich aus der Wartezeit in der Warteschlange und der Dauer zum Überqueren der Kreuzung und ist in Sekunden angegeben. Mit eingerechnet sind auch Autos die noch in der Warteschlange stehen solange sie vor der Endzeit angekommen sind. Die angesetzten 1=9 für R1R 3 sowie 2=5 für R2 R4 hoch erwiesen und die Warteschlange ist sehr schnell gewachsen.

4

haben sich als zu

Beispielausgabe: Die Dauer der Beobachtung = 20 Minuten 100 mal Wiederholt Mit lambda1/3 = 9 und lambda 2/4 = 5 Durchschnittliche Verweildauer im System in Sekunden 161.61 Durchschnittliche Warteschlangenlänge von 0.5367 Autos pro Sekunde Durchschnittliche Warteschlangenlänge von 32.204 Autos pro Minute

Ein befriedigenderes Ergebniss konnte mit λ1= 3 für R1+R3 sowie λ2 = 2 für R2+R4 erziehlt werden. Beispielausgabe: Die Dauer der Beobachtung = 20 Minuten 100 mal Wiederholt Mit lambda1/3 = 3 und lambda 2/4 = 2 Durchschnittliche Verweildauer im System in Sekunden 24.68 Durchschnittliche Warteschlangenlänge von 0.0285 Autos pro Sekunde Durchschnittliche Warteschlangenlänge von 1.7092 Autos pro Minute

5.Berechnung der Standardabweichung: In der Simulation wird die Standardabweichung von der Warteschlangenlänge berechnet. Von jeder Ausführung der Simulation wird die Warteschlangenlänge in einer Liste abgespeichert und die Standardabweichung über die bekannte Formel

berechnet. N ist hierbei wie oft wiederholt wurde und xi die Warteschlangenlänge bzw deren Mittelwert.

Teil B: Varianzreduktionstechniken Man versucht die Effizienz zu erhöhen , indem mit weniger Simulationsaufwand die gleiche Präzision erreicht wird. Es werden im folgenden zwei Techniken vorgestellt.

1. Die Common Random Numbers Die Common Random Numbers werden zum Vergleich von zwei Modellen eingesetzt. Als Vergleichsmodell zu der im ersten Teil beschriebenen Kreuzung, wird eine „normale Vorfahrtskreuzung“ verwendet, bei der die Autos aus Richtung1 und Richtung3 die Vorfahrtsstrasse verwenden. (Da laut Aufgabenstellung die Wahrscheinlichkeit für Geradeausfahren bei 50% liegt ist die Warteschlangenlänge sowie die Wartezeit insgesamt kürzer wie in der ürsprünglichen Kreuzung.) Die Idee der Common Random Numbers ist die Verwendung gemeinsamer Zufallszahlen für verschiedene Prozesse. Mit hoher Wahrscheinlichkeit trifft dies auf folgende Zufallszahlenströme zu. Die der Ankunftszeiten und der Bedienzeiten. Damit dies korrekt umgesetzt werden kann muss eine eindeutige Zuordnung der jeweiligen Prozesse, wie Ankunftszeit zu Ankunftszeit und Bedienzeit zu Bedienzeit usw. möglich sein. Probleme können entstehen, wenn unterschiedliche Systemkonfigurationen verglichen werden, oder Zufallszahlen aus unterschiedlichen Verteilungen kommen (und damit mit unterschiedlicher Zahl von [0,1)-verteilten Werten generiert werden). Beide Probleme werden versucht zu minimieren, indem die gleichen Systemkonfigurationen genutzt werden und indem alle synchonisierten Zufallszahlen ursprünglich zwischen [0,1)- verteilt sind. Erreicht werden identische Zufallszahlen über Seedwerte, welche wiederum für jeden Simulationslauf zufällig erzeugt und dann gespeichert werden. In meiner Simulation sind die Ankunftszeiten identisch, sowie in welche Richtung ein Auto fährt. Die Bedienzeit, hier die Dauer zum überqueren der Kreuzung ist schon konstant und braucht nicht weiter synchronisiert zu werden.

Was genau bewirkt nun die Synchronisation der Zufallsströme? Hier ein Beispiel: X1j : ist die mittlere Durchlaufzeit des j-ten Simulationslaufs im System Fifo. X2j : ist die mittlere Durchlaufzeit des j-ten Simulationslaufs im System Normale Vorfahrtsstrasse. Zj = X1j – X2j : Zielgröße zum Vergleich Duch die gemeinsamen Zufallszahlen wird nun versucht eine positive Korrelation zu erzeugen ( diese wird von meinem Programm zur Veranschaulichung berechnet ) wodurch die Varianz reduziert wird. Dies wird an folgende Formel sehr deutlich : Var [Z(n)] = Var(Zj) / n = [Var(X1j) + Var(X2j) – 2*Cov(X1j,X2j) ] / n Bei unabhängigen Zufallszahlen gilt Cov(X1,X2) =0. Sind X1 und X2 hingegen positiv korreliert, gilt Cov(X1,X2) >0 , sodass die Varianz des Schätzers Z(n) reduziert werden kann. Um die Standardabweichung der Zielgröße Zj zu schätzen, muss von dieser noch der Mittelwert errechnet werden. Die Varianz lässt sich dann direkt über folgende Formel schätzen.

Folgende Ausgabe ist aus dem Buch von Kelton :

Leider kann man nicht vorhersagen wie gross die Reduktion sein wird, da dies zu stark von dem jeweiligen gewählten Modell abhängt. Aussagekräftiger ist hierbei die Monotonie. Verhalten sich die Zufallsströme mit einer gemeinsamen Zufallszahl zueinander monoton-fallend oder steigend erwartet man eine Varianzreduktion. Verhalten sich die Ströme gegensätzlich kann hingegen die Varianz vergrössert werden.

Beispiel: positive Korrelation

negative Korrelation

Die negative Korrelation wird auch als Backfire bezeichnet. Eine alternative Methode zur oberen Tabelle, um zu schauen ob CRN auch wirklich in unserem Beispiel funktioniert, ist das S 2z  nS 21 nS 22 n auch wirklich gilt. Hier ein Beispiel aus meiner Simulation: Die Dauer der Beobachtung = 475 Minuten 200 mal Wiederholt Mit lambda1/3 = 9 und lambda 2/4 = 5 Standardabweichung pro Sekunde S1^2 = 0.2393 Standardabweichung pro Sekunde S2^2 = 0.2749 Gemeinsame Varianz mit CRN S12^2 = 0.1985 ohne CRN S1^2+ S2^2 = 0.5142 Korrelation = 0.62 Ersparnis 61.0%

Aus meinen Beobachtungen folgt allerdings, das es etliche Beispiele gibt wo CRN gar nicht oder nur schlecht funktioniert. Dies wird durch unterschiedliche Lambda bzw. unterschiedlicher Beobachtungsdauer hervorgerufen. Die Dauer der Beobachtung = 10 Minuten 2000 mal Wiederholt Mit lambda1/3 = 5 und lambda 2/4 = 3 Ersparnis = 22.0% und

Mit lambda1/3 = 9 und lambda 2/4 = 5

Ersparnis = 45.0%

2. Stratified Sampling a) Die Idee: Im stratified sampling wird die Grundgesamtheit in Schichten (Strata) eingeteilt, so dass die Anzahl in den einzelnen Schichten getrennt erhöht werden kann. Es wird eine konditionierte Stichprobenziehung zb. Gleichverteilter Zufallsvariablen erstellt, um nicht das Pech zu haben, dass die realisierten Werte U alle in einem Abschnitt liegen. Deswegen unterteilen wir die Region in i disjunkte Subregionen Α1, Α2 ... Ak , mit den Wahrscheinlichkeiten p k = P(Y ∈ Αk ) mit k = 1,...,j. Die Zufallsvariablen werden nun quasi “gezwungen” sich gleichmäßig auf die Subregionen Α1, Α2 ... Ak zu verteilen.

Mit Y als Funktion erhalten wir nun : E[Y] = p 1 * E[Y| Y ∈ Α1] +

p 2 * E[Y| Y ∈ Α2] + ... +

p k * E[Y| Y ∈ Αk]

( C.1 )

J

∑ pk

=

* E[Y| Y ∈ Αk] ( C.1 )

k =1

Der einfachste Fall ist proportional Sampling, wobei hier die Anzahl gezogener Beobachtungen der jeweiligen Schicht Ak von deren Wahrscheinlichkeit pk abhängt. Beträgt die frei wählbare Samplegrösse n, werden n k =np k unabhängige Schichten erzeugt ( np k wird auf eine Ganzzahl abgerundet). Für jede Schicht und den daraus gezogenen Stichproben lässt sich ein erwartungstreuer Schätzer von E[Y| Y ∈ Αk] mit (Yk1 + ... + Yk nk) / n k angeben. J

Dieser lautet : Y =

∑ P Y ∈ Ak 

* E[Y| Y ∈ A k ]

∑ P Y ∈ Ak 

*

k =1 J

=

nk

k =1 J

=

∑ pk

nk

∑ Y ks

( C.2 ) n k =np k

s=1

J

1/n

mit n k =np k und p k = P(Y ∈ Αk) folgt

s=1

1/ n∗ p k 

*

k =1

Y=

∑ Y ks

1/ n k

nk

∑ ∑ Y ks k =1 s=1

Wir definieren nun m Zufallsvariablen mit Xj ∼X ~ Un(Aj) (eine jede ZV ist also gleichverteilt in “seiner” Subregion) mit den Wahrscheinlichkeiten pj = P(X ∈ Αj ) = 1/m (j = 1 .. m). X wird als stratifizierte Zufallszahl bezeichnet.

Damit erhalten wir:

E[Y] = p1 * E[Y| X ∈ Α1] + p2 * E[Y| X ∈ Α2] + ... + pk * E[Y| X ∈ Αk] j

=

∑ pk

* E[Y| X ∈ Αk]

( C.3 )

k =1

Die Hauptprobleme bei dieser Methode bestehen in der richtigen Auswahl der Schichten, der Bestimmung der Versuchszahl pro Schicht und – falls nicht vorher bekannt – auch in der Bestimmung der Wahrscheinlichkeiten für die Schichten.

b) Ein Beispiel anhand der Gleichverteilung : Das Intervall von (0,1) wird in Partitionen zerlegt, so genannte Stratas, mit Wahrscheinlichkeit 1/k wobei k=0,...,n die Intervallanzahl darstellt

Aus jedem Stratum wird die gleiche Anzahl an Zufallszahlen gezogen. Dies geschieht mit Hilfe einer Transformation, der im Intervall (0,1) verteilten Zufallszahlen Ui.

Wir schreiben für V i

V i =a k U bk −a k  schreiben. ( C.4 )

auch

Einteilung der Gleichverteilung in 5 Intervalle. Nun werden fünf gleichverteilte Zufallszahlen gezogen Diese in V i eingesetzt ergibt nun :

U i ( 0.47, 0.12, 0.09, 0.71, 0.34).

V i = (0.094, 0.224, 0.418, 0742, 0868) Dies kann nun beliebig wiederholt werden, um für jedes Stratum beliebig viele Zufallszahlen zu generieren. Man erkennt das die resultierende Verteilung ’eher’ der gewünschten Gleichverteilung entspricht. Simulationen mit den modifizierten Zufallszahlen haben nun eine reduzierte Varianz.

c) Die Varianzreduktion mit Stratified Sampling : wir führen folgende Abkürzungen ein ○ µ = E[Y], ○ µ = E[Y| Y ∈ A k ] ○  2 = Var [ Y | Y ∈ A k ] 1/n k Damit erhalten wir für den Erwartungswert mit (C.1) J

∑ pk

E[Y] =

k =1 J

nk

∑ p k k

=

1/n k ∑ E [Y ki ]

*

i=1

=µ

k =1

Die Varianz berechnet sich nun als : J

Var[Y] =

∑p k =1 J

=

nk

Var [1 /nk ∑ Y ks ] s=1

∑  p 2k 2k nk /n 2k k =1 J

=

2 k

∑p k =1

2 k

J

 / nk = 1/n ∑ p k  2k 2 k

( mit n k =np k ( C.5 )

k =1

Vergleichen wir dies mit der Varianz von einer normalen Stichprobe: J

E[ Y 2 ] =

J

∑

p k * E[ Y 2 | Y∈ Ak ] =

k =1

J

erhalten wir mit µ =

∑ p k k k =1

∑ p k 2k 2k  k =1

( C.6 )

J

Var[Y] = E[ Y

2

]– 

2

=

∑ pk  k =1

2 k

J

∑ pk 

+

k =1

( C.7 )

( C.8 )

2 k

J

2

-  ∑ pk  k  k=1

( C.9 )

Da C.7 und C.8 die normale Varianz darstellen gilt J

J

k =1

k=1

∑ p k 2k≥∑ pk k 

2

Man sieht das C.7 + C.8 aus C.6 kommen und die Varianz innerhalb der Schicht darstellen. Der Teil C.9 verschwindet und stellt die Varianz über die Schichten dar. Daraus folgt , dass stratified sampling mit einer Gleichverteilung die Varianz nur verringern kann.

d) Beispiel an einem Poissonprozess. Zuerst wird die Verteilungsfunktion F(x) der Poisson-Verteilung gebildet, woraus sich die Schichten ergeben. z.b. F(x) = (0,3 ; 0,66 ; 0,88 ; 0,97 ; 0,99) Jede dieser Schichten hat die Wahrscheinlichkeit

p k =a k −a k −1  .

Nun können pro Schicht wieder beliebig viele Gleichverteilte stratifizierte Zufallszahlen U~[0.1] erzeugt werden, wie im obigen Beispiel. Mit der Formel aus (C.4) erhalten wir : V =a k −1 U∗a k −ak −1  mit unterer Grenze a  k−1 und oberer Grenze a k . Aus den erstellten Schichten werden die zur Erstellung des Poissonprozess benötigten Zufallszahlen mit Wahrscheinlichkeit p k gezogen. Beispielausgabe: Die Dauer der Beobachtung = 20 Minuten 500 mal Wiederholt Mit lambda1/3 = 7 und lambda 2/4 = 4 Die modifizierte Kreuzung ohne Stratifizierung Durchschnittliche Verweildauer im System in Sekunden 91.51 Durchschnittliche Warteschlangenlänge von 0.2697 Autos pro Sekunde Durchschnittliche Warteschlangenlänge von 16.1824 Autos pro Minute Standardabweichung pro Sekunde S1^2 = 0.0057 Die modifizierte

Kreuzung mit Stratifizierung

Durchschnittliche Verweildauer im System in Sekunden 22.74 Durchschnittliche Warteschlangenlänge von 0.0799 Autos pro Sekunde Durchschnittliche Warteschlangenlänge von 4.7963 Autos pro Minute Standardabweichung pro Sekunde S2^2 = 0.0020 Ersparnis 65.0%

Man sieht das die Warteschlangenlänge sehr stark reduziert wurde. Dies muss aber nicht

zwangsläufig auch zu einer geringeren Varianz der Warteschlangenlänge führen. Vielmehr ist die Beobachtungsdauer und die Ankunftsrate interessant, wie eine andere Eingabe belegt. Nun Mit lambda1/3 = 9 und lambda 2/4 = 5 ohne Stratifizierung Durchschnittliche Warteschlangenlänge von 32.6067 Autos pro Minute Standardabweichung pro Sekunde S1^2 = 0.0078 mit Stratifizierung Durchschnittliche Warteschlangenlänge von 19.0629 Autos pro Minute Standardabweichung pro Sekunde S2^2 = 0.019

Die Varianz verschlechtert sich hier um -40%, obwohl die Warteschlangenlänge reduziert wurde. Nun Mit lambda1/3 = 4 und lambda 2/4 = 3

Ersparnis von 80 %.

Meine Erklärung wäre das sich die Poisson-Verteilung bei großen Lambda gleichmäßiger verteilt und die Einteilung sehr fein ist. Der Gewinn fällt dann nicht mehr so hoch aus, was allerdings keine Erklärung für -40% Varianz sein kann.

( bei langen Beobachtungszeiten reicht der vom System zugewiesene Speicher von 64MB für Stratified Sampling nicht mehr aus, dieser muss manuell auf 256MB erhöht werden)

3.)Fazit: Common Random Numbers sind leicht zu implementieren, allerdings wird eine zweite Vergleichskreuzung benötigt, die auch jedes mal neu berechnet werden muss. Wir haben also einen hohen Rechenaufwand. Ein großer Vorteil ist allerdings, dass man sichergehen kann, dass sich das Ergebnis verbessert wenn keine gegensätzliche Monotonie herrscht. Stratified Sampling : Nicht so einfach zu implementieren, dafür braucht man hier aber eigentlich keine zwei Kreuzungen ( nur zum Vergleich ) und damit erheblich weniger Rechenleistung. Außerdem kann man eher eine Verbesserung erwarten als bei den CRN. Diese fällt dann auch meist höher aus. Quellen : Glasserman, P.(2004). Monte Carlo Methods in financial engineering. Law, A.M., und Kelton, W:D (2000). Simulation modelling and analysis.

Ein Nachtrag zur Bedienung : ist CRN auf false gesetzt , wird automatisch das 2te Verfahren mit Stratified Sampling genutzt.