MCMC - Simulationsgerhold/pub_files/sem15/s_mueller.pdf · MCMC - Simulation Seminar f¨ur Finanz-...

23
MCMC - Simulation Seminar f¨ ur Finanz- und Versicherungsmathematik Wintersemester 2015/16 Betreuer Autor Dr. Stefan Gerhold Christoph M¨ uller e1226165

Transcript of MCMC - Simulationsgerhold/pub_files/sem15/s_mueller.pdf · MCMC - Simulation Seminar f¨ur Finanz-...

Page 1: MCMC - Simulationsgerhold/pub_files/sem15/s_mueller.pdf · MCMC - Simulation Seminar f¨ur Finanz- und Versicherungsmathematik Wintersemester 2015/16 Betreuer Autor Dr. Stefan Gerhold

MCMC - Simulation

Seminar fur Finanz- und VersicherungsmathematikWintersemester 2015/16

Betreuer AutorDr. Stefan Gerhold Christoph Muller

e1226165

Page 2: MCMC - Simulationsgerhold/pub_files/sem15/s_mueller.pdf · MCMC - Simulation Seminar f¨ur Finanz- und Versicherungsmathematik Wintersemester 2015/16 Betreuer Autor Dr. Stefan Gerhold

Inhaltsverzeichnis

1 MCMC und die Kryptographie 3

2 Einfuhrung in die Theorie der Markovketten 5

3 MCMC-Simulationen 10

4 Hard disk in a box 13

5 Mathematische Ausfuhrung 15

6 MCMC in der Statistik 17

7 Zusammenfassung und Schlussfolgerung 20

8 Literaturverzeichnis 21

9 Abbildungsverzeichnis 22

1

Page 3: MCMC - Simulationsgerhold/pub_files/sem15/s_mueller.pdf · MCMC - Simulation Seminar f¨ur Finanz- und Versicherungsmathematik Wintersemester 2015/16 Betreuer Autor Dr. Stefan Gerhold

Einleitung

Simulationen sind in der Mathematik nicht mehr wegzudenken. Viele komplexe Fragestellun-gen konnen nur numerisch berechnet werden. Die folgende Arbeit wird sich mit der MCMC-Simulation beschaftigen. Durch die Einteilung der mathematischen Fragestellung in Zustande,welche mit zufalligen Wahrscheinlichkeiten wechseln, wird versucht das Problem auf Markovket-ten zu transferieren. Mit Hilfe des Gesetzes der großen Zahlen kann mit zunehmender Stichpro-beanzahl ein genaueres Ergebnis erzielt werden. Zunachst werde ich ein Beispiel zur Anwendungdieser Methode geben und anschließend auf die Theorie und das Themengebiet genauere einge-hen.

2

Page 4: MCMC - Simulationsgerhold/pub_files/sem15/s_mueller.pdf · MCMC - Simulation Seminar f¨ur Finanz- und Versicherungsmathematik Wintersemester 2015/16 Betreuer Autor Dr. Stefan Gerhold

1 MCMC und die Kryptographie

Zuerst gibt es einen kurzen Abstecher in die Kryptographie, um mit der MCMC-Simulationvertraut zu werden. In der Kryptographie ist das Ziel verschlusselte Texte zu entschlusseln. Essoll versucht werden dieses Problem mathematisch zu losen. Hier ein Beispiel.

Abbildung 1.1: Verschhlusselter Text

Um folgenden Text zu entschlusseln benotigt man eine Funktion f mit:

f : coded → decoded

welche fur alle Zeichen im verschlusselten Text Buchstaben oder Satz- und Sonderzeichenzuruckgibt, wobei diese zwei Gruppen fortan als Symbole deklariert werden. Um dieses Pro-blem zu losen ist zunachst eine Statistik von Noten. Man nimmt sich ein Standardliteraturwerkund zahlt alle Eins-zu-Eins Ubergange von Symbolen. Man startet mit dem ersten Buchstabendes Alphabets und zahlt wie oft im Text ein A nachgestellt ist. Dies wird fur alle Symbolkombi-nationen durchgefuhrt und anschließend eine Matrix erstellt, welche so angeordnet ist, dass dieAnzahl von dem Startsymbol in der Zeile aufgeschrieben sind. Es entsteht eine Ubergangsmatrixvon Eins-zu-Eins Ubergangen. Man erhalt eine quadratische Matrix A ∈ Nn×n. In der deut-schen Sprache ist die Dimension dieser Matrix sehr groß, da man Groß- und Kleinbuchstabenunterscheiden muss sowie alle Satz- und Sonderzeichen betrachten muss. Dies ist notwendig, daman nicht weiß ob bei der Verschlusselung Zeichen Groß- und Kleinschreibung vernacchlassigtwurden oder ob die Satz- und Sonderzeichen verschlusselt wurden. Falls ersichtlich wird, dassSatzzeichen nicht verschlusselt wurden, konnte das n kleiner gemacht werden. Nun normierenwir die Zeilen der Matrix wie folgt

∀i, j ∈ 1, . . . , n : bij = aij∑iaij

damit erhalten wir eine stochastische Matrix B ∈ Rn×n, welche die Ubergangswahrscheinlichkeit

3

Page 5: MCMC - Simulationsgerhold/pub_files/sem15/s_mueller.pdf · MCMC - Simulation Seminar f¨ur Finanz- und Versicherungsmathematik Wintersemester 2015/16 Betreuer Autor Dr. Stefan Gerhold

der Eins-zu-Eins Ubergange angibt.Die Entschlusselung erfolgt nach dem Trail-and-Error Prinzip. Deshalb wird eine Moglichkeitbenotigt, um festzustellen wie gut eine Entschlusselung tatsachhlich ist. Da schon eine Ubergangsmatrixvorhanden ist, kann diese zur Klassifizierung der Entschlusselung herangzogen werden. Eine ge-eignete Funktion ware

W (f) =∑i

Bf(si),f(si+1)

wobei si alle Zeichen des verschlusselten Textes durchlauft.Nun kann folgender Algorithmus gestartet werden.

I Wahle ein fI Berechne W (f)I durch Mischen von Zuordnungen erstelle ein f∗

I Berechne W (f∗)I wenn W (f∗) >W (f) wahle f = f∗ und starte erneutI falls W (f) > W (f∗) wirf eine Munze und bei Kopf gehe zu f∗ und bei Zahl bleibe bei f

Diese Ruckschrittsmoglichkeit ist wichtig, da man damit verhindern kann, dass der Algorithmusnicht in einem lokalen Maximum stecken bleibt sondern das globale Maximum erreichen kann.Der Algorithmus versucht f mittels zufalligen Vertauschungen zu verbessern. Die Laufzeit die-ses Algorithmus kann sehr stark schwanken. In den meisten Fallen werden ein paar tausendDurchlaufe benotigt, um das globale Maximum zu erreichen. Der in Abbildung 1.1 gezeigteverschlusselte Text, wird nach Anwendung des Algorithmus entschlusselt:

Abbildung 1.2: Entschlusselter Text

In diesem Fall handelt es sich um ein interessantes Beispiel. Der Text wurde von einem Gefangnisinsassenverschlusselt und sollte verschickt werden. Die Gefangniswarter konnten den Text abfangenund er wurde nach obriger Methode entschlusselt. Die Sinnhaftigkeit der Ubersetzung kannman leicht erkennen. Der Text ist eine Mischung aus Englisch, Spanisch und einem eigenenGefangnisjargon.Im folgenden Kapitel wird es eine allgemeine Einfuhrung in die Theorie der Markovketten ge-ben und wichtige Satze und Definitionen werden formuliert, um das MCMC-Verfahren bessernachvollziehen zu konnen.

4

Page 6: MCMC - Simulationsgerhold/pub_files/sem15/s_mueller.pdf · MCMC - Simulation Seminar f¨ur Finanz- und Versicherungsmathematik Wintersemester 2015/16 Betreuer Autor Dr. Stefan Gerhold

2 Einfuhrung in die Theorie der Markov-ketten

Definiton 2.1. Sei χ ein endlicher Zustandsraum. Sei (Xt)t∈N ein zeitdiskreter stochastischerProzess. Dann ist Xt eine Markovkette genau dann, wenn

P (Xt+1 = it+1|X0 = i0, . . . , Xt = it) = P (Xt+1 = it+1|Xt = it)

gilt.

Diese Eigenschaft wird als die Markoveigenschaft bezeichnet. Eine wichtige Darstellungsformfur Markovketten ist die stochastische Matrix.

Definiton 2.2. Sei A ∈ Rn×n. A ist genau dann eine stochastische Matrix, wenn

∀i, j : Aij > 0∀j :

∑j

Aij = 1.

Nun lasst sich ein Markovprozess, wie folgt mit einer stochastischen Matrix charakterisieren:

Definiton 2.3. Sei χ endlicher Zustandsraum mit |χ| = n. So beschreibt eine stochastischeMatrix A ∈ Rn×n eine Markovkette wobei gilt, dass Aij = P (Xt+1 = j|Xt = i) ist.

Einzig eine Startverteilung muss noch angegeben werden, sodass die Markovkette wohldefiniertist. Um einen wichtigen Satz fur die weitere Arbeit einzufuhren, braucht man zunachst einewichtige Eigenschaft.

Definiton 2.4. Eine Markovkette Xt auf dem Zustandsraum χ heißt irreduzibel genau, dannwenn ∀i, j ∈ χ ein s, t ∈ N mit s < t existiert, sodass

P (Xt = j|Xs = i) > 0.

5

Page 7: MCMC - Simulationsgerhold/pub_files/sem15/s_mueller.pdf · MCMC - Simulation Seminar f¨ur Finanz- und Versicherungsmathematik Wintersemester 2015/16 Betreuer Autor Dr. Stefan Gerhold

Das bedeutet, dass jeder Zustand in endlichen Schritten erreicht werden kann. Ein wichtigerSatz fur die Betrachtung der Konvergenz von Markovketten ist der folgende Satz:

Satz 2.1. Sei Xt eine irreduzible Markovkette auf dem endlichen Zustandsraum χ und A seidie stochastische Matrix welche die Markovkette beschreibt. Dann existiert ein einfacher Eigen-wert λ = 1 zur Matrix A.

Beweis. Zunachst wird gezeigt, dass Potenzen von stochastischen Matrizen wieder eine sto-chastische Matrix ist. Wir betrachten nun

cij = (A2)ij =∑r

airarj .

Da A eine stochastische Matrix ist sind ihre Eintrage nichtnegativ. Daraus folgt auch, dass dieEintrage von A2 nichtnegativ sind. Zu zeigen bleibt noch, dass die Zeilensumme gleich 1 ist.Wir betrachten∑

j

cij =∑j

∑r

airarj =∑r

∑j

airarj =∑r

air∑j

arj = 1.

Damit ist gezeigt, dass A2 stochastisch ist. Analog kann man die Behauptung auch fur allePotenzen von A zeigen. Nach dem Satz von Perron Frobenius weiß man, dass der Spektralra-dius ein positiver und einfacher Eigenwert ist und er erfullt die Gleichung %(A) = lim

n→∞‖An‖

1n .

Wie vorhin gezeigt sind die Potenzen von A stochastisch, damit auch limn→∞

An. Falls fur dieMatrixnorm nun die Zeilensummennorm gewahlt wird ist ‖An‖ = 1 damit folgt mit dem Satzvon Perron Frobenius die Behauptung.

Damit folgt, dass ein Eigenvektor zum Eigenwert 1 existiert. Dieser wird im folgenden immer mitπ bezeichnet. Fur diesen gilt Aπ = π. Zu beachten ist, dass falls man einmal diesen Vektor alsVerteilung angibt, verlasst man die Verteilung nicht mehr, da A ∗Aπ = Aπ = π. Dieser Vektorwird als stationare Verteilung bezeichnet. Eine weitere Eigenschaft, welche die Existenz einerstationaren Verteilung impliziert, ist die Detailed-Balance-Eigenschaft. Diese lautet Aijπi =Ajiπj . Diese Eigenschaft ist eine Voraussetzung fur die Definiton von MCMC-Simuationen.Das wichtigste Resultat das man wir fur spatere Konvergenzbetrachtungen benotigen, ist, dasslimn→∞

(An)ij = πj . Dies bedeutet, dass alle Zustande gegen die stationare Verteilung konvergie-ren. Das motiviert den folgenden Satz:

Satz 2.2. Sei Xt eine Markovkette und A die zugehorige Ubergangsmatrix. Weiters sei derZustandsraum χ endlich. Dann gilt:

limn→∞

(An)ij = πj ∀i, j ∈ χ.

Beweis. Wie zuvor bewiesen, existiert ein Eigenwert λ = 1 und der dazugehorige Eigenvektorπ. Der Beweis fur die Konvergenz ist das Grundprinzip der Vektoriteration. Bei dieser wird einebeliebige Startverteilung x hergenommen und es wird x1 = Ax berechnet und iterativ fortge-setzt. Dann konvergiert xt gegen π. Dies ist die Behauptung die zu zeigen ist. Da der Beweis sehrumfangreich ist, wird er in dieser Arbeit nicht behandelt. Er kann im Skriptum ”NumerischeMathematik A”von Lothar Nannen auf Seite 106 nachgeschlagen werden.

6

Page 8: MCMC - Simulationsgerhold/pub_files/sem15/s_mueller.pdf · MCMC - Simulation Seminar f¨ur Finanz- und Versicherungsmathematik Wintersemester 2015/16 Betreuer Autor Dr. Stefan Gerhold

Als nachstes wird die Konvergenz der Markovkette betrachtet. Die Frage ist, wie schnell konver-giert die Markovkette. Dies wird im folgenden Satz behandelt. Fur die Norm um die Konvergenzzu betrachten, wird die Norm der Totalen Variation herangezogen. Diese ist wie folgt definiert:

‖Ani − π‖TV = 12∑j

|Anij − πj |.

Satz 2.3. Sei A Ubergangsmatrix und es existiert eine stationare Verteilung. Dann existierenEigenvektoren und Eigenwerte, sodass:

4‖Ani − π‖2TV 6|χ|−1∑i=0

β2ni ψ2

i .

Beweis. Wir fuhren folgendes Skalarprodukt auf dem Raum L2(π) ein

〈g, h〉 =∑x

g(x)h(x)πx,

wobei L2(π) die Menge aller g : χ −→ R ist. Dann operiert K auf dem L2 mit:

Ag(x) =∑y

g(y)Axy,

nun folgt aus:

〈Ag, h〉 =⟨∑

j

g(j)Aij , h⟩

=∑i

∑j

g(j)Aijh(i)πi

=∑i

∑j

g(j)πjAijh(i)

=∑j

g(j)πj∑i

Ajih(i)

= 〈g,Ah〉 .

Nun folgt, dass A ein selbstadjungierender Operator ist. Nach dem Spektralsatz folgt nun, dasses eine ONB von Eigenvektoren gibt und es gilt:

∀x ∈ χ : Ax =|χ|−1∑i=0

βi 〈x, ψi〉ψi,

mit Kψi = βiψi. Nun folgt mit x als das neutrale Element des Vektorraums:

A =|χ|−1∑i=0

βiψi(x)ψi(y)π(y) = π(y)|χ|−1∑i=0

βiψi(x)ψi(y),

weiters sei:

A2 = πy

|χ|−1∑i=0

βiψi(x)ψi(y)Axy = πy

|χ|−1∑i=0

β2i ψi(x)ψi(y),

7

Page 9: MCMC - Simulationsgerhold/pub_files/sem15/s_mueller.pdf · MCMC - Simulation Seminar f¨ur Finanz- und Versicherungsmathematik Wintersemester 2015/16 Betreuer Autor Dr. Stefan Gerhold

induktiv folgt nun:

An(x, y) = π(y)|χ|−1∑i=0

βni ψi(x)ψi(y),

schließlich folgt mit der Cauchy-Schwarz Ungleichung:

4‖Anx − π‖2TV 6∑y

(Anxy − πy)2

πy=|χ|−1∑i=0

β2ni ψ2

i .

Man kann nun sehen, dass die Konvergenz von den Eigenwerten und Eigenvektoren abhangt.Falls eine besser Konvergenzrate erreicht werden soll, dann mussen geeignete Manipulationenfur die EW und EV gemacht werden.Es ist einfach die Ubergangswahrscheinlichkeiten, jedoch die stationare Verteilung oftmals schwie-rig zu berechnen. Genau aus diesem Grund wird auf die MCMC-Simulation zuruckgegriffen. DieMCMC-Simulation wird im nachsten Kapitel vorgestellt.Zunachst folgt noch eine kurze Einfuhrung in die Theorie der Markovketten auf allgemeinenZustandsraumen. Bis hierher wurden nur endliche Zustandsraume betrachtet. Um Markovket-ten zu verallgemeinern, braucht man den Markovkern.

Definiton 2.5. Sei Ψ ein beliebiger Raum mit einer σ-Algebra F . Ein Kern ist eine Abbil-dung P : (Ψ,F)→ [0, 1], welche folgende Eigenschaften erfullt:

· P (x, ·) ist eine Wahrscheinlichkeit auf (Ψ,F) fur alle x ∈ Ψ,· P (·, A) ist eine messbare Funktion fur alle A ∈ F .

Der Markovkern P (x,A) gibt also eine Wahrscheinlichkeit an, mit der man von dem Zustand xzum Ereignis A gelangt. Der Markovkern kann nun wie folgt iteriert werden:

P 2(x,A) =∫P (z,A)P (x, dz).

Diese Gleichung sagt aus, dass man von x nach A in zwei Schritten gelangt, wobei uber dieWahrscheinlichkeiten aller moglichen Zwischennzustande z integriert wird. Eine stationare Ver-teilung erfullt dann folgende Gleichung:

π(A) =∫P (x,A)π(dx).

Im spateren Verlauf der Arbeit wird diese Verallgemeinerung verwendet, um komplexere Bei-spiele zu losen.

8

Page 10: MCMC - Simulationsgerhold/pub_files/sem15/s_mueller.pdf · MCMC - Simulation Seminar f¨ur Finanz- und Versicherungsmathematik Wintersemester 2015/16 Betreuer Autor Dr. Stefan Gerhold

Zur Veranschaulichung sieht man in der nachsten Grafik den Graphen einer Markovkette.

Abbildung 2.1: Markovkettengraph

Wie man in diesem Graphen sieht, sind die einzelnen Zustande mit einer Ubergangswahrscheinlichkeitverbunden. Diese Markovkette ist auch irreduzibel, da es keine getrennte Kommunikationsklassegibt. Die Ubergangsmatrix wurde wie folgt lauten:

A =

0 15

45

12 0 1

2110

910 0

.Man sieht, dass die Zeilensumme immer 1 betragt. Da sie irreduzibel ist, existiert die stationareVerteilung. Die stationare Verteilung ist π = (0.22, 0.4, 0.38). Im nachsten Kapitel wird derMCMC-Algorithmus eingefuhrt.

9

Page 11: MCMC - Simulationsgerhold/pub_files/sem15/s_mueller.pdf · MCMC - Simulation Seminar f¨ur Finanz- und Versicherungsmathematik Wintersemester 2015/16 Betreuer Autor Dr. Stefan Gerhold

3 MCMC-Simulationen

In diesem Kapitel werden die MCMC-Simulationen eingefuhrt und Motivation und Beispielebehandelt. MCMC steht fur Markov-Chains-Monte-Carlo. Die MCMC-Simulation wird verwen-det um eine Markovkette zu erzeugen, welche eine gewunschte Wahrscheinlichkeitsverteilungals stationare Verteilung hat. Diese erzeugten Markovketten erfullen alle die Detailed-Balance-Eigenschaft, welche im vorigen Kapitel vorgestellt wurde. Die Eigenschaft impliziert die Existenzeiner stationaren Verteilung.Wir betrachten nun einen Metropolis-Hastings-Algorithmus:

Satz 3.1. Sei χ ein endlicher Zustandsraum und π(x) eine Wahrscheinlichkeit. Sei Axy dieUbergangsmatrix einer Markovkette auf χ mit Axy > 0⇔ Ayx > 0. Sei Bxy := min1, πyAyx

πxAxy.

Dann sei

Kxy :=AxyBxy fur x 6= yJxy +

∑zAxz(1−Bxz) fur x = y.

Dann erfulle Kxy folgende Gleichung: πxKxy = πyKyx

Beweis. fur x = y trivial erfullt. Fur x 6= y sei oBdA Byx = 1 (ansonsten ist Bxy = 1)πxKxy = πxAxy

πyAyx

πxAxy= πyAyx = πyKyx.

Nun sieht man, dass wir Markovkette, die anfangs keine Relation zu der stationaren Verteilunghatte, zu einer neuen Markovkette umgewandelt wurde, welche die gewunschte stationare Ver-teilung aufweist.Zunachst wird auf die Entstehung und die Motivation des Metropolisalogrithmus und desMetropolis-Hastings-Algorithmus eingegangen. Der Metropolisalogrithmus geht auf NicholasMetropolis zuruck. 1953 veroffentlichte dieser eine Puplikation, in welcher er Markovketten ver-wendete um eine Boltzmann-Verteilung zu generieren. Er stellte folgenden Algorithmuus auf:

I Man startet im Ort xi, wobei dies den Ortsvektor nach dem i-ten Iterationsschrittbezeichnet.

I Es wird ein neuer Ort y vorgeschlagen welcher folgende Gleichung erfullt y = xi + rqiwobei r fester Suchradius ist und q ein Zufallsvektor mit q(i) ∈ [−1, 1] ∀i.

I Nun wird die Energiedifferenz ∆E = E(y)− E(xi) berechnet.I Falls die Energiedifferenz ∆E ≤ 0 wird y als neuer Ort akzeptiert und xi+1 = y.I ansonsten wird y mit Wahrscheinlichkeit p := min(1, exp

(−∆EκT

)), wobei T die

Temperatur des Systems und κ die Boltzmannkonstante ist, akzeptiert.

10

Page 12: MCMC - Simulationsgerhold/pub_files/sem15/s_mueller.pdf · MCMC - Simulation Seminar f¨ur Finanz- und Versicherungsmathematik Wintersemester 2015/16 Betreuer Autor Dr. Stefan Gerhold

Dieser Algorihtmus generiert eine stationare Verteilung welche der Boltzman-Verteilung ent-spricht. In der nachsten Grafik sieht man eine Boltzman-Verteilung:

Abbildung 3.1: Boltzmann-Verteilung

Wie man sieht ist die Boltzmann-Verteilung fur niedrigere Temperaturen nach links gestaucht.W. Keith Hastings verallgemeinerte diesen Algorithmus, in dem er die Vorschlagsdichte P (x|y)generiert, welche Vorschlage mit Wahrscheinlichkeit p = min(1, W (y)P (x|y)

W (x)P (y|x)) akzeptiert und auchvom nachsten Zustand abhangig ist. Wobei W (·) eine beliebige Wahrscheinlichkeitsverteilungist.Nun wird wieder der in Satz 3.1 definierte Algorithmus herangezogen. Es wurde gezeigt, dassdie Detailed-Balance-Eigenschaft gilt. Desweitern folgt nun, dass

(πK)y =∑x

πxKxy =∑x

πyKyx = πy∑x

Kyx = πy,

damit folgt, dass π der Linkseigenvektor zu K ist. Damit folgt mit Satz 2.1 und 2.2, dass diedurch den Algorithmus erzeugte Markovkette K gegen die stationare Verteilung konvergiert.Dies angewendet auf das Kryptographiebeispiel im 1. Kapitel ergibt folgenden Algorithmus.Zuvor benotigt man noch ein paar Eigenschaften. Sei χ die Menge aller Eins-zu-Eins Ubergangevom Raum der verschlusselten Zeichen zum Raum der Symbole. Die stationare Verteilung hatfolgende Form:

π(f) = z−1∏i

B(f(si), (f(si+1)),

wobei z der Normierungsfaktor ist. M ist die Matrix der Einschrittubergange der Sprache diefur das Enschlusseln verwendet wird. Der Normierungsfaktor lautet wie folgt:

z =∑f

∏i

B(f(si), (f(si+1)).

11

Page 13: MCMC - Simulationsgerhold/pub_files/sem15/s_mueller.pdf · MCMC - Simulation Seminar f¨ur Finanz- und Versicherungsmathematik Wintersemester 2015/16 Betreuer Autor Dr. Stefan Gerhold

Es wird uber alle moglichen Enschlusselungsfunktionen f summiert. Dies macht die Bestimmungvon z fast unmoglich. Ein weiteres Problem stellt die Große des Zustandsraums χ dar. Diese zweiProblematiken versucht man mit der MCMC-Simulation zu umgehen. Sei nun die Machtigkeitder Menge aller Eins-zu-Eins Ubergange im verschlusselten Text gleich m und Machtigkeit derMenge aller Eins-zu-Eins Ubergange der Symbole gleich n mit n ≥ m.Dann ist |χ| = n(n − 1) . . . (n −m + 1). Diese Zahl kann sehr groß werden, falls beispielsweisem = n = 50 ist. Fur den Algorithmus generiert man J(f, f∗) was einen zufalligen Tausch vonzwei Symbolen darstellt. Damit folgt:

J(f, f∗) := 1

n(n−1)(m−n+2)(m−n+1) wenn f, f∗sich in maximal zwei Stellen unterscheiden0 sonst.

Aufgrund von J(f, f∗) = J(f∗, f) folgt fur die Akzeptierwahrscheinlichkeit, dass B(f, f∗) =π(f∗)π(f) ist.

12

Page 14: MCMC - Simulationsgerhold/pub_files/sem15/s_mueller.pdf · MCMC - Simulation Seminar f¨ur Finanz- und Versicherungsmathematik Wintersemester 2015/16 Betreuer Autor Dr. Stefan Gerhold

4 Hard disk in a box

In diesem Kapitel wird eine weitere Anwendung fur das MCMC-Verfahren gezeigt. Bei ”Harddisk in a box”gibt es folgende Ausgangssituation. Sei n die Anzahl der Scheiben mit Radiusr in dem Einheitsquadrat. Es wird nun versucht die Teilchenbewegung mithilfe des MCMC-Algorithmus zu beschreiben. Dies wird wie folgt gemacht:

I Man startet mit x ∈ τ(n, r).I Man nehme nun den Mittelpunkt einer Scheibe mit Wahrscheinlichkeit 1

n .I Nun wird ein Punkt im Radius h ausgewahlt.I Verschiebe den Mittelpunkt zu dem Punkt.I Falls diese neue Position in τ(n, r) liegt so akzeptiere diese Bewegung.

Wobei τ(n, r) die Topologie der moglichen Positionen der Scheiben ist. Das Problem ist dasWahlen einer uniform-verteilten Stichporbe von x ∈ τ(n, r). Dieser Algorithmus verschiebtzufallig Koordinaten. Falls nun X1, X2, . . . die erfolgreichen Positionen sind so wird, Xk furkleine k und r uniform-verteilt sein. Fur große k kann man eine Funktion f mit f : τ(n, r)→ Rfinden, wobei man

∫τ(n,r)

f(x)dx mit 1k

k∑i=1

f(Xk) approximieren kann.

Um dieses f sinnvoll zu bestimmen, muss zunachst die Frage gestellt werden, wie diese aussehenkonnen. Die Motivation fur dieses Beispiel kommt von der Betrachtung von Phasendiagrammen.Diese geben den Verlauf der Aggregatzustande im Verhaltnis zu Druck und Temperatur an. Einsolches sieht man in folgender Abbildung:

Abbildung 4.1: Phasendiagramm

13

Page 15: MCMC - Simulationsgerhold/pub_files/sem15/s_mueller.pdf · MCMC - Simulation Seminar f¨ur Finanz- und Versicherungsmathematik Wintersemester 2015/16 Betreuer Autor Dr. Stefan Gerhold

Diese Phasendiagramme sind bereits in der Forschung intensiv analysiert worden. Es gibt eineKurve endlicher Lange zwischen dem flussigen und gasformigen Zustand, welche im Tripletpunktstartet. An diesem coexistieren alle drei Zustand. Der Endpunkt dieser Kurve ist der kritischePunkt. Man kann eine gegen Unendlich laufende fest-flussig Kurve erkennen. Eine interessanteFolgerung dieses Phasendiagramms ist, dass egal welche Temperatur vorgegeben wird, man mitdem richtigen Druck eine fest-flussige Transformation, oder umgekehrt erzeugen kann. Ab einemgewissen Druck sind die starken innermolekularen abstoßenden Krafte nicht mehr relevant. Inder nachsten Grafik sieht man wie die Scheibenverteilung bei gewissen Druck aussieht.

Abbildung 4.2: Scheibenverteilung

η bezeichnet die dynamische Viskositat. Dies entspricht Druck in einer Zeiteinheit. Je hoherdiese ist, desto hoher ist der Druck. Man sieht nun, dass fur hoheren Druck die Scheiben eineGitterform annehmen. Eine sinnvolle Funktion f welche oben gesucht wurde, lasst sich nunaufgrund der Gitterstruktur darstellen. Diese lautet:

f(x) =

∣∣∣∣∣∣ 1N

N∑j=1

1Nj

∑k

e6iθjk

∣∣∣∣∣∣ ,wobei uber die Anzahl N summiert wird. Nj ist die Anzahl der Nachbarn der j-ten Scheibe.θjk ist der Winkel zwischem der j-ten und k-ten Scheibe. Werte fur diese Funktion konnenmit Simulationen erzeugt werden. Im nachsten Abschnitt wird dieses Problem mathematischgenauer ausgefuhrt.

14

Page 16: MCMC - Simulationsgerhold/pub_files/sem15/s_mueller.pdf · MCMC - Simulation Seminar f¨ur Finanz- und Versicherungsmathematik Wintersemester 2015/16 Betreuer Autor Dr. Stefan Gerhold

5 Mathematische Ausfuhrung

Zunachst sei Ω ⊆ Rn ein verbundenes offenes Gebiet. p∗(x) > 0 ∀x ∈ Ω und z=∫

Ω p∗(x)dx <∞.

Dann ist z−1 die Normierungskonstante und p∗(x) = z−1p(x). Dieser Normmierungsfaktor istschwierig zu bestimmen. Deshalb sollte man einen MCMC-Algorithmus finden, bei dem mannicht z berechnen muss. Eine Moglichkeit ware folgender Algoritmus:

I Starte in x und wahle y ∈ A(h) wobei A eine messbare Menge ist und A ⊆ Ω wobeiΩ ein Lipschitzgebiet ist.

I Berechne p(y).I falls p(y) ≥ p(x) dann gehe zu y.I falls p(y) < p(x) gehe nach y mit Wahrscheinlichkeit p(y)

p(x) .I andernfalls bleibe in x.

Dieser Algoritmus setzt nicht voraus, dass man z direkt bestimmen muss. Da auf allgemei-nen Raumen gearbeitet wird, muss ein Markovkern aufgestellt, welcher den Ubergang von xzu y darstellt. Dafur mussen einige Dinge beachtet werden. Sei φ(z) = 1

V (A(h))δA(h), sodass∫φ(z)dz = 1. Sei p eine positiv beschrankte Funktion mit

∫Ωp(z)dz = 1. Fur h ∈ [0, 1] gilt nun,

dass:

Kh,p(x, y) = h−dφ

(x− yh

)min

(p(y)p(x) , 1

).

Nun lasst sich der Markovkern durch folgende Gleichung darstellen:

P (x, dy) = m(x)δx +∫Ω

Kh,p(x, y)δydy,

wobei

m(x) = 1−∫Ω

Kh,p(x, y)dy.

Falls man nun alles einsetzt, so kommt man schließlich auf folgenden Kern:

P (x, dy) = m(x)δx + h−d

V (A(h))δA(h)

(x− yh

)min

(p(y)p(x) , 1

)dy

mit

m(x) = 1−∫Rd

h−d

V (A(h))δA(h)

(x− yh

)min

(p(y)p(x) , 1

)dy,

wobei h der Radius ist, in welchem y zu x liegt. δx ist die Dichte am Ort x und V (A(h)) ist dasVolumen von A(h). Nun kann der Algorithmus wie folgt dargestellt werden.

15

Page 17: MCMC - Simulationsgerhold/pub_files/sem15/s_mueller.pdf · MCMC - Simulation Seminar f¨ur Finanz- und Versicherungsmathematik Wintersemester 2015/16 Betreuer Autor Dr. Stefan Gerhold

Wahle X0 = x ∈ A(h) −→ wahle X1 ∈ P (X0, dy) −→ wahle X2 ∈ P (X1, dy).Damit folgt aus Def. 2.5, dass

P (X2 ∈ B) =∫RdP (z,B)P (x, dz)

iterativ folgt dann weiters

P (Xk ∈ B) =∫RdP (z,B)P k−1(x, dz)

fur geeignete h konvergiert nun P (Xk ∈ B) gegen∫B p(y)dy. Der Algorithmus hat eine Konver-

genzrate die man wie folgt darstellen kann:

|P (Xk ∈ B)−∫Bp(y)dy| ≤ c1e

−c2kh2,

wobei p(x) beschrankt sein muss und B Teilmenge eines Lipschitzgebietes ist. c1, c2 sind positiveZahlen. Den Beweis fur diese Abschatzung kann man in der Ausarbeitung ”Geometric Analysisfor the Metropolis Algorithm on Lipschitz Domains” von Persi Diaconis, Gilles Lebeau undLaurent Michel nachlesen.

16

Page 18: MCMC - Simulationsgerhold/pub_files/sem15/s_mueller.pdf · MCMC - Simulation Seminar f¨ur Finanz- und Versicherungsmathematik Wintersemester 2015/16 Betreuer Autor Dr. Stefan Gerhold

6 MCMC in der Statistik

Simulationen sind eine haufig verwendete Methode zur Erzeugung von Zufallszahlen. In derStatistik ist die Manipulation von Stichproben und Betrachtung ihrer Verteilungsfunktionenein zentrales Thema. Es gibt einige Moglichkeiten solche Zufallszahlen zu erzeugen. Methodendafur sind die Inversionsmethode, oder das Box-Muller-Verfahren. In diesem Kapitel soll dassogenannte Gibbs-Sampling betrachten werden. Das Gibbs-Sampling geht auf die MathematikerS. und D. German zuruck. Benannt wurde es nach dem amerikanischen Physiker J. W. Gibbs.Nachstehend folgt eine kurze Einfuhrung in die Bayes-Statistik.

Definiton 6.1. A-priori Verteilung: Die Verteilung einer sG X ∼ f(x; θ) wobei θ ∈ Θ. Inder Bayes-Statistik wird der Faktor θ nicht als fester unbekannter Parameter betrachtet, son-dern als Realisation der sG θ mit der Dichte h(θ). Diese wird als A-priori Dichte bezeichnetund dementsprechend erhalt man durch Integration die A-priori Verteilung.

Definiton 6.2. Sei X = (X1, X2, . . . , Xn) eine Stichprobe von X ∼ f(x; θ) so ist die bedingteDichte von X|θ = θ gleich f(x|θ) =

∏if(xi|θ). Dann ist die gemeinsame Dichte von X und θ

g(x, θ) = f(x|θ)h(θ). Die Randdichte ist dann wie folgt definiert

g1(x) =∞∫−∞

g(x, θ)dθ =∞∫−∞

f(x|θ)h(θ)dθ.

Definiton 6.3. A-posteriori Verteilung: Sei die Randdichte wie oben, dann ist die A-posterioriDichte definiert als

k(θ|x) = f(x|θ)h(θ)g1(x) = f(x|θ)h(θ)

∞∫−∞

f(x|θ)h(θ)dθ.

Die A-posteriori Verteilung ist also asymptotischverteilt nach der Likelihoodfunktion mal derA-priori Verteilung. Zunachst wird ein kleines Motivationsbeispiel fur Gibbs-Sampling gezeigt.Sei X ∼ fX(x) und Y ∼ fY (y). Es wird angenommen, dass die Simulation der Beobachtungenvon X schwierig ist, jedoch die von Y und X|Y=y aber nicht.

17

Page 19: MCMC - Simulationsgerhold/pub_files/sem15/s_mueller.pdf · MCMC - Simulation Seminar f¨ur Finanz- und Versicherungsmathematik Wintersemester 2015/16 Betreuer Autor Dr. Stefan Gerhold

Dann kann man dieses Problem umgehen indem man:I Erzeuge Y ∼ fY (y).I Erzeuge X ∼ fX|Y (x|Y = y).

Das erzeugte X hat die Dichte fX(x). Die Behauptung folgt aus:

P (X ≤ t) = E(FX|Y (t)) =∞∫−∞

t∫−∞

fX|Y (x|y)dx

fY (y)dy

=t∫

−∞

∞∫−∞

fX|Y (x|y)fY (y)dy

dx=

t∫−∞

∞∫−∞

fX,Y (x, y)dy

dx=

t∫−∞

fX(x)dx.

Nun kann man den Erwartungswert einer Funktion W = g(X) von X bestimmen. Man erzeugtzunachst Paare (X1, Y1), . . . , (Xm, Ym) und es gilt nach dem Gesetz der großen Zahlen:

W = 1m

m∑i=1

g(Xi)P−−→ E(g(X)).

Das Gibbs-Sampling hat ein ahnliches Prinzip. Fur Erzeugung von sG nach f(x,y) werden wech-selweise Beobachtungen von f(x|y) und f(y|x) herangezogen. Nun wird der bivariante Fall desGibbs-Samplings betrachtet:I Erzeugung des Startwerts (x0, y0) fur i=0.I Erzeuge Xi+1 ∼ f(x|yi).I Erzeuge Yi+1 ∼ f(y|xi+1).I i=i+1 und wiederhole.

Die Paare (X1, Y1), . . . , (Xm, Ym) bilden eine Markovkette, da

P [(Xi, Yi) = (ii, ji)|(Xi−1, Yi−1) = (ii−1, ji−1), . . . , (X1, Y1) = (i1, j1)] =P [(Xi, Yi) = (ii, ji)|(Xi−1, Yi−1) = (ii−1, ji−1)].

Also ist das Paar (Xi, Yi) nur von dem vorherigen Paar (Xi−1, Yi−1) abhangig. Fur i→∞ folgtnun fur die Gibbs-Sampler:

XiD−−→ fX(x),

YiD−−→ fY (y).

Die Konvergenzgeschwindigkeit fur die Kette kann stark varieren. Ublicherweise wird zur Schatzungzunachst die ersten m Elemente nicht verwendet und dann erst die nachsten n-m Paare. Dieswird als Einbrennen der Kette bezeichnet. Der Grund dafur ist, dass diese Paare sich anfangsnoch stark von der stationaren Verteilung unterscheiden konnen. Je mehr Schritte der Algo-rithmus lauft, desto besser wird die stationare Verteilung erreicht. Deshalb lasst man die erstenm Paare weg, da diese die großte Varianz zur stationaren Verteilung besitzen. Demnach ist eingeeigneter Schatzer fur den Erwartungswert von X:

E(X) = 1n−m

n∑i=m+1

XiP−−→ E(X).

18

Page 20: MCMC - Simulationsgerhold/pub_files/sem15/s_mueller.pdf · MCMC - Simulation Seminar f¨ur Finanz- und Versicherungsmathematik Wintersemester 2015/16 Betreuer Autor Dr. Stefan Gerhold

Fur die Randdicht wurde sich dann folgender Schatzer anbieten:

fX(x) = 1n−m

n∑i=m+1

f(x|yi).

Dieser bivariante Sampler kann einfach auf einen mehrdimensionalen Fall erweitert werden. DieHauptanwendung fur den Gibbs-Sampler in der Bayes-Statistik, ist das sogenannte hierachischeModell. Bei diesem gibt es folgende Problemstellung. Da die A-priori Verteilung einen großenEinfluss auf die Ereignisse hat, versucht man eine sG zu finden, welche die A-priori Verteilungmodelliert. Das Modell sieht wie folgt aus:

X|θ ∼ f(x|θ),θ|γ ∼ h(θ|γ),

γ ∼ ψ(γ).

Der Hilfsparameter γ wird als Hyperparameter bezeichnet. In diesem Fall sehen die Gibbs-Sampler wie folgt aus:

θi|x, γi−1 ∼ g(θ|x, γi−1),γ|x, θi ∼ g(γ|x, θi).

Fur die Grenzwerte dieser Gibbs-Sampler gilt nun:

θiD−−→ k(θ|x),

γiD−−→ g(γ|x),

1n

n∑i=1

u(θi)P−−→ E(u(θ)|x).

Nach dem Einbrennen der Markovkette schatzt man dann den Erwartungswert von u(θ) wiefolgt ab:

1n−m

n∑i=1

u(θi), n > m.

19

Page 21: MCMC - Simulationsgerhold/pub_files/sem15/s_mueller.pdf · MCMC - Simulation Seminar f¨ur Finanz- und Versicherungsmathematik Wintersemester 2015/16 Betreuer Autor Dr. Stefan Gerhold

7 Zusammenfassung und Schlussfolgerung

Markovketten sind mathematische Modelle, welche fur zufallsbedingte Ereigniswechsel eine guteApproximation liefern. Ihre Theorie ist gut erforscht und dadurch ist es leicht diese anzuwenden.Es ist oft eine Erleichterung eine komplexe Fragestellung auf Markovketten transformieren zukonnen.

MCMC-Simulationen nutzen die Theorie der Markovketten aus und versuchen diese zum Losenvon komplexen Aufgabenstellungen zu verwenden. Die Anwendung der MCMC-Simulation istweit gefachert. In dieser Arbeit wurde bewiesen, dass sie zur Entschlusselung von Texten ver-wendet werden konnen. MCMC kann jedoch auch verwendet werden, um Teichlenbewegungund deren resultierenden Zustande zu beschreiben. Deshalb findet dieses Verfahren auch inder Thermodynamik Verwendung. Man konnte auch andere physikalische Phanomene damitzu beschreiben versuchen,z.B. Ladungsverteilung in elektrischen Bauteilen. Naheliegend ist dieAnwendung in der Statistik. Das Gibbs-Sampling ist eine der wichtigsten Methoden zur Erzeu-gung von Stichproben in der Bayes-Statistik.

Die MCMC-Simulation hat wie man sieht eine große Bedeutung in der Mathematik und derPhysik.

20

Page 22: MCMC - Simulationsgerhold/pub_files/sem15/s_mueller.pdf · MCMC - Simulation Seminar f¨ur Finanz- und Versicherungsmathematik Wintersemester 2015/16 Betreuer Autor Dr. Stefan Gerhold

8 Literaturverzeichnis

[1] Diaconis, Persi; Lebeau, Gilles; Michel, Laurent; Geometric Analysis for the MetropolisAlgorithm on Lipschitz Domains.[2] Waldmann, K.-H., 2013; Stochastische Modelle Eine anwendungsorientierte Einfuhrung.

[3] Diaconis, Persi, 2009; The Markov Chain Monte Carlo Revolution.

[4] Neal, Radford M., 2011; MCMC using Hamiltonian dynamics.

[5] Zeuner, Michelle, 2011; Die Monte-Carlo-Methode.

[6] Robert, C. P., Casella, G., 2004; Monte Carlo Statistical Methods.

[7] Beichl, Isabel, 2000; The Metropolis Algorithm.

[8] Schomaker, Jens; Einfuhrung in die Theorie der Markov-Ketten.

[9] Spellecchia, Claudia, 2008; MCMC.

[10] Gurker, Werner, Version 2015; Angewandte Mathematische Statistik.

21

Page 23: MCMC - Simulationsgerhold/pub_files/sem15/s_mueller.pdf · MCMC - Simulation Seminar f¨ur Finanz- und Versicherungsmathematik Wintersemester 2015/16 Betreuer Autor Dr. Stefan Gerhold

9 Abbildungsverzeichnis

Abbildung 1.1: Verschlusselter Text;Diaconis, Persi; The Markov Chain Monte Carlo Revolution.

Abbildung 1.2: Entschlusselter Text;Diaconis, Persi; The Markov Chain Monte Carlo Revolution.

Abbildung 2.1: Markovkettengraph;URL: https://de.wikipedia.org/wiki/Markow-Kette, 24.02.2016

Abbildung 3.1: Boltzmann-Verteilung;

URL: http://www.chemieunterricht.de/dc2/fragen/kf-ka-279.htm, 24.02.2016Abbildung 4.1: Phasendiagramm;

URL: http://www.lernort-mint.de/Anorganische%20Chemie/Anwendungen/phasendiagramm.html, 24.02.2016Abbildung 4.2: Scheibenverteilung;

Diaconis, Persi; The Markov Chain Monte Carlo Revolution.

22