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

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

    09-Aug-2020
  • Category

    Documents

  • view

    2
  • download

    0

Embed Size (px)

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

  • MCMC - Simulation

    Seminar für Finanz- und VersicherungsmathematikWintersemester 2015/16

    Betreuer AutorDr. Stefan Gerhold Christoph Müller

    e1226165

  • Inhaltsverzeichnis

    1 MCMC und die Kryptographie 3

    2 Einführung in die Theorie der Markovketten 5

    3 MCMC-Simulationen 10

    4 Hard disk in a box 13

    5 Mathematische Ausführung 15

    6 MCMC in der Statistik 17

    7 Zusammenfassung und Schlussfolgerung 20

    8 Literaturverzeichnis 21

    9 Abbildungsverzeichnis 22

    1

  • Einleitung

    Simulationen sind in der Mathematik nicht mehr wegzudenken. Viele komplexe Fragestellun-gen können nur numerisch berechnet werden. Die folgende Arbeit wird sich mit der MCMC-Simulation beschäftigen. Durch die Einteilung der mathematischen Fragestellung in Zustände,welche mit zufälligen 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. Zunächst werde ich ein Beispiel zur Anwendungdieser Methode geben und anschließend auf die Theorie und das Themengebiet genauere einge-hen.

    2

  • 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 verschlüsselte Texte zu entschlüsseln. Essoll versucht werden dieses Problem mathematisch zu lösen. Hier ein Beispiel.

    Abbildung 1.1: Verschhlüsselter Text

    Um folgenden Text zu entschlüsseln benötigt man eine Funktion f mit:

    f : {coded} → {decoded}

    welche für alle Zeichen im verschlüsselten Text Buchstaben oder Satz- und Sonderzeichenzurückgibt, wobei diese zwei Gruppen fortan als Symbole deklariert werden. Um dieses Pro-blem zu lösen ist zunächst eine Statistik von Nöten. Man nimmt sich ein Standardliteraturwerkund zählt alle Eins-zu-Eins Übergänge von Symbolen. Man startet mit dem ersten Buchstabendes Alphabets und zählt wie oft im Text ein A nachgestellt ist. Dies wird für alle Symbolkombi-nationen durchgeführt und anschließend eine Matrix erstellt, welche so angeordnet ist, dass dieAnzahl von dem Startsymbol in der Zeile aufgeschrieben sind. Es entsteht eine Übergangsmatrixvon Eins-zu-Eins Übergängen. Man erhält 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 Verschlüsselung Zeichen Groß- und Kleinschreibung vernacchlässigtwurden oder ob die Satz- und Sonderzeichen verschlüsselt wurden. Falls ersichtlich wird, dassSatzzeichen nicht verschlüsselt wurden, könnte 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 Übergangswahrscheinlichkeit

    3

  • der Eins-zu-Eins Übergänge angibt.Die Entschlüsselung erfolgt nach dem Trail-and-Error Prinzip. Deshalb wird eine Möglichkeitbenötigt, um festzustellen wie gut eine Entschlüsselung tatsächhlich ist. Da schon eine Übergangsmatrixvorhanden ist, kann diese zur Klassifizierung der Entschlüsselung herangzogen werden. Eine ge-eignete Funktion wäre

    W (f) =∑i

    Bf(si),f(si+1)

    wobei si alle Zeichen des verschlüsselten Textes durchläuft.Nun kann folgender Algorithmus gestartet werden.

    I Wähle ein fI Berechne W (f)I durch Mischen von Zuordnungen erstelle ein f∗I Berechne W (f∗)I wenn W (f∗) >W (f) wähle f = f∗ und starte erneutI falls W (f) > W (f∗) wirf eine Münze und bei Kopf gehe zu f∗ und bei Zahl bleibe bei f

    Diese Rückschrittsmöglichkeit 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 zufälligen Vertauschungen zu verbessern. Die Laufzeit die-ses Algorithmus kann sehr stark schwanken. In den meisten Fällen werden ein paar tausendDurchläufe benötigt, um das globale Maximum zu erreichen. Der in Abbildung 1.1 gezeigteverschlüsselte Text, wird nach Anwendung des Algorithmus entschlüsselt:

    Abbildung 1.2: Entschlüsselter Text

    In diesem Fall handelt es sich um ein interessantes Beispiel. Der Text wurde von einem Gefängnisinsassenverschlüsselt und sollte verschickt werden. Die Gefängniswärter konnten den Text abfangenund er wurde nach obriger Methode entschlüsselt. Die Sinnhaftigkeit der Übersetzung kannman leicht erkennen. Der Text ist eine Mischung aus Englisch, Spanisch und einem eigenenGefängnisjargon.Im folgenden Kapitel wird es eine allgemeine Einführung in die Theorie der Markovketten ge-ben und wichtige Sätze und Definitionen werden formuliert, um das MCMC-Verfahren bessernachvollziehen zu können.

    4

  • 2 Einführung 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 Darstellungsformfür 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 lässt 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 für die weitere Arbeit einzuführen, braucht man zunächst 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

  • Das bedeutet, dass jeder Zustand in endlichen Schritten erreicht werden kann. Ein wichtigerSatz für 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. Zunächst 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 Einträge nichtnegativ. Daraus folgt auch, dass dieEinträge 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 für allePotenzen von A zeigen. Nach dem Satz von Perron Frobenius weiß man, dass der Spektralra-dius ein positiver und einfacher Eigenwert ist und er erfüllt die Gleichung %(A) = lim

    n→∞‖An‖

    1n .

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

    An. Falls für dieMatrixnorm nun die Zeilensummennorm gewählt 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. Für diesen gilt Aπ = π. Zu beachten ist, dass falls man einmal diesen Vektor alsVerteilung angibt, verlässt man die Verteilung nicht mehr, da A ∗Aπ = Aπ = π. Dieser Vektorwird als stationäre Verteilung bezeichnet. Eine weitere Eigenschaft, welche die Existenz einerstationären Verteilung impliziert, ist die Detailed-Balance-Eigenschaft. Diese lautet Aijπi =Ajiπj . Diese Eigenschaft ist eine Voraussetzung für die Definiton von MCMC-Simuationen.Das wichtigste Resultat das man wir für spätere Konvergenzbetrachtungen benötigen, ist, dasslimn→∞

    (An)ij = πj . Dies bedeutet, dass alle Zustände gegen die stationäre Verteilung konvergie-ren. Das motiviert den folgenden Satz:

    Satz 2.2. Sei Xt eine Markovkette und A die zugehörige Übergangsmatrix. Weiters sei derZustandsraum χ endlich. Dann gilt:

    limn→∞

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

    Beweis. Wie zuvor bewiesen, existiert ein Eigenwert λ = 1 und der dazugehörige Eigenvektorπ. Der Beweis für 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

  • Als nächstes wird die Konvergenz der Markovkette betrachtet. Die Frage ist, wie schnell konver-giert die Markovkette. Dies wird im folgenden Satz behandelt. Für 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 Übergangsmatrix und es existiert eine stationäre Verteilung. Dann existierenEigenvektoren und Eigenwerte, sodass:

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

    β2ni ψ2i .

    Beweis. Wir führen 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

  • 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 ψ2i .

    Man kann nun sehen, dass die Konvergenz von den Eigenwerten und Eigenvektoren abhängt.Falls eine besser Konvergenzrate erreicht werden soll, dann müssen geeignete Manipulationenfür die EW und EV gemacht werden.Es ist einfach die Übergangswahrscheinlichkeiten, jedoch die stationäre Verteilung oftmals schwie-rig zu berechnen. Genau aus diesem Grund wird auf die MCMC-Simulation zurückgegriffen. DieMCMC-Simulation wird im nächsten Kapitel vorgestellt.Zunächst folgt noch eine kurze Einführung in die Theorie der Markovketten auf allgemeinenZustandsräumen. Bis hierher wurden nur endliche Zustandsräume 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 erfüllt:

    · P (x, ·) ist eine Wahrscheinlichkeit auf (Ψ,F) für alle x ∈ Ψ,· P (·, A) ist eine messbare Funktion für 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 über dieWahrscheinlichkeiten aller möglichen Zwischennzustände z integriert wird. Eine stationäre Ver-teilung erfüllt dann folgende Gleichung:

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

    Im späteren Verlauf der Arbeit wird diese Verallgemeinerung verwendet, um komplexere Bei-spiele zu lösen.

    8

  • Zur Veranschaulichung sieht man in der nächsten Grafik den Graphen einer Markovkette.

    Abbildung 2.1: Markovkettengraph

    Wie man in diesem Graphen sieht, sind die einzelnen Zustände mit einer Übergangswahrscheinlichkeitverbunden. Diese Markovkette ist auch irreduzibel, da es keine getrennte Kommunikationsklassegibt. Die Übergangsmatrix würde wie folgt lauten:

    A =

    0 15 4512 0 12110

    910 0

    .Man sieht, dass die Zeilensumme immer 1 beträgt. Da sie irreduzibel ist, existiert die stationäreVerteilung. Die stationäre Verteilung ist π = (0.22, 0.4, 0.38). Im nächsten Kapitel wird derMCMC-Algorithmus eingeführt.

    9

  • 3 MCMC-Simulationen

    In diesem Kapitel werden die MCMC-Simulationen eingeführt und Motivation und Beispielebehandelt. MCMC steht für Markov-Chains-Monte-Carlo. Die MCMC-Simulation wird verwen-det um eine Markovkette zu erzeugen, welche eine gewünschte Wahrscheinlichkeitsverteilungals stationäre Verteilung hat. Diese erzeugten Markovketten erfüllen alle die Detailed-Balance-Eigenschaft, welche im vorigen Kapitel vorgestellt wurde. Die Eigenschaft impliziert die Existenzeiner stationären Verteilung.Wir betrachten nun einen Metropolis-Hastings-Algorithmus:

    Satz 3.1. Sei χ ein endlicher Zustandsraum und π(x) eine Wahrscheinlichkeit. Sei Axy dieÜbergangsmatrix einer Markovkette auf χ mit Axy > 0⇔ Ayx > 0. Sei Bxy := min{1, πyAyxπxAxy }.Dann sei

    Kxy :={AxyBxy für x 6= yJxy +

    ∑zAxz(1−Bxz) für x = y.

    Dann erfülle Kxy folgende Gleichung: πxKxy = πyKyx

    Beweis. für x = y trivial erfüllt. Für 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 stationären Verteilunghatte, zu einer neuen Markovkette umgewandelt wurde, welche die gewünschte stationäre Ver-teilung aufweist.Zunächst wird auf die Entstehung und die Motivation des Metropolisalogrithmus und desMetropolis-Hastings-Algorithmus eingegangen. Der Metropolisalogrithmus geht auf NicholasMetropolis zurück. 1953 veröffentlichte 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 erfüllt 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

  • Dieser Algorihtmus generiert eine stationäre Verteilung welche der Boltzman-Verteilung ent-spricht. In der nächsten Grafik sieht man eine Boltzman-Verteilung:

    Abbildung 3.1: Boltzmann-Verteilung

    Wie man sieht ist die Boltzmann-Verteilung für niedrigere Temperaturen nach links gestaucht.W. Keith Hastings verallgemeinerte diesen Algorithmus, in dem er die Vorschlagsdichte P (x|y)generiert, welche Vorschläge mit Wahrscheinlichkeit p = min(1, W (y)P (x|y)W (x)P (y|x)) akzeptiert und auchvom nächsten Zustand abhängig 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 stationäre Verteilung konvergiert.Dies angewendet auf das Kryptographiebeispiel im 1. Kapitel ergibt folgenden Algorithmus.Zuvor benötigt man noch ein paar Eigenschaften. Sei χ die Menge aller Eins-zu-Eins Übergängevom Raum der verschlüsselten Zeichen zum Raum der Symbole. Die stationäre Verteilung hatfolgende Form:

    π(f) = z−1∏i

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

    wobei z der Normierungsfaktor ist. M ist die Matrix der Einschrittübergänge der Sprache diefür das Enschlüsseln verwendet wird. Der Normierungsfaktor lautet wie folgt:

    z =∑f

    ∏i

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

    11

  • Es wird über alle möglichen Enschlüsselungsfunktionen f summiert. Dies macht die Bestimmungvon z fast unmöglich. Ein weiteres Problem stellt die Größe des Zustandsraums χ dar. Diese zweiProblematiken versucht man mit der MCMC-Simulation zu umgehen. Sei nun die Mächtigkeitder Menge aller Eins-zu-Eins Übergänge im verschlüsselten Text gleich m und Mächtigkeit derMenge aller Eins-zu-Eins Übergänge 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. Für den Algorithmus generiert man J(f, f∗) was einen zufälligen 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 unterscheiden

    0 sonst.

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

    12

  • 4 Hard disk in a box

    In diesem Kapitel wird eine weitere Anwendung für 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 1n .I Nun wird ein Punkt im Radius h ausgewählt.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 möglichen Positionen der Scheiben ist. Das Problem ist dasWählen einer uniform-verteilten Stichporbe von x ∈ τ(n, r). Dieser Algorithmus verschiebtzufällig Koordinaten. Falls nun X1, X2, . . . die erfolgreichen Positionen sind so wird, Xk fürkleine k und r uniform-verteilt sein. Für 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 zunächst die Frage gestellt werden, wie diese aussehenkönnen. Die Motivation für dieses Beispiel kommt von der Betrachtung von Phasendiagrammen.Diese geben den Verlauf der Aggregatzustände im Verhältnis zu Druck und Temperatur an. Einsolches sieht man in folgender Abbildung:

    Abbildung 4.1: Phasendiagramm

    13

  • Diese Phasendiagramme sind bereits in der Forschung intensiv analysiert worden. Es gibt eineKurve endlicher Länge zwischen dem flüssigen und gasförmigen Zustand, welche im Tripletpunktstartet. An diesem coexistieren alle drei Zuständ. Der Endpunkt dieser Kurve ist der kritischePunkt. Man kann eine gegen Unendlich laufende fest-flüssig Kurve erkennen. Eine interessanteFolgerung dieses Phasendiagramms ist, dass egal welche Temperatur vorgegeben wird, man mitdem richtigen Druck eine fest-flüssige Transformation, oder umgekehrt erzeugen kann. Ab einemgewissen Druck sind die starken innermolekularen abstoßenden Kräfte nicht mehr relevant. Inder nächsten Grafik sieht man wie die Scheibenverteilung bei gewissen Druck aussieht.

    Abbildung 4.2: Scheibenverteilung

    η bezeichnet die dynamische Viskosität. Dies entspricht Druck in einer Zeiteinheit. Je höherdiese ist, desto höher ist der Druck. Man sieht nun, dass für höheren Druck die Scheiben eineGitterform annehmen. Eine sinnvolle Funktion f welche oben gesucht wurde, lässt sich nunaufgrund der Gitterstruktur darstellen. Diese lautet:

    f(x) =

    ∣∣∣∣∣∣ 1NN∑j=1

    1Nj

    ∑k

    e6iθjk

    ∣∣∣∣∣∣ ,wobei über 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 für diese Funktion könnenmit Simulationen erzeugt werden. Im nächsten Abschnitt wird dieses Problem mathematischgenauer ausgeführt.

    14

  • 5 Mathematische Ausführung

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

    Ω p∗(x)dx

  • Wähle X0 = x ∈ A(h) −→ wähle X1 ∈ P (X0, dy) −→ wähle 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)

    für 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−c2kh

    2,

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

    16

  • 6 MCMC in der Statistik

    Simulationen sind eine häufig verwendete Methode zur Erzeugung von Zufallszahlen. In derStatistik ist die Manipulation von Stichproben und Betrachtung ihrer Verteilungsfunktionenein zentrales Thema. Es gibt einige Möglichkeiten solche Zufallszahlen zu erzeugen. Methodendafür 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 zurück. Benannt wurde es nach dem amerikanischen Physiker J. W. Gibbs.Nachstehend folgt eine kurze Einführung 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 erhält 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. Zunächst wird ein kleines Motivationsbeispiel für 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

  • 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 erzeugtzunächst 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 ähnliches Prinzip. Für 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) für 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) abhängig. Für i→∞ folgtnun für die Gibbs-Sampler:

    XiD−−→ fX(x),

    YiD−−→ fY (y).

    Die Konvergenzgeschwindigkeit für die Kette kann stark varieren. Üblicherweise wird zur Schätzungzunächst die ersten m Elemente nicht verwendet und dann erst die nächsten n-m Paare. Dieswird als Einbrennen der Kette bezeichnet. Der Grund dafür ist, dass diese Paare sich anfangsnoch stark von der stationären Verteilung unterscheiden können. Je mehr Schritte der Algo-rithmus läuft, desto besser wird die stationäre Verteilung erreicht. Deshalb lässt man die erstenm Paare weg, da diese die größte Varianz zur stationären Verteilung besitzen. Demnach ist eingeeigneter Schätzer für den Erwartungswert von X:

    Ē(X) = 1n−m

    n∑i=m+1

    XiP−−→ E(X).

    18

  • Für die Randdicht würde sich dann folgender Schätzer anbieten:

    f̄X(x) =1

    n−m

    n∑i=m+1

    f(x|yi).

    Dieser bivariante Sampler kann einfach auf einen mehrdimensionalen Fall erweitert werden. DieHauptanwendung für 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).

    Für 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 schätzt man dann den Erwartungswert von u(θ̃) wiefolgt ab:

    1n−m

    n∑i=1

    u(θi), n > m.

    19

  • 7 Zusammenfassung und Schlussfolgerung

    Markovketten sind mathematische Modelle, welche für 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 zukönnen.

    MCMC-Simulationen nutzen die Theorie der Markovketten aus und versuchen diese zum Lösenvon komplexen Aufgabenstellungen zu verwenden. Die Anwendung der MCMC-Simulation istweit gefächert. In dieser Arbeit wurde bewiesen, dass sie zur Entschlüsselung von Texten ver-wendet werden können. MCMC kann jedoch auch verwendet werden, um Teichlenbewegungund deren resultierenden Zustände zu beschreiben. Deshalb findet dieses Verfahren auch inder Thermodynamik Verwendung. Man könnte auch andere physikalische Phänomene 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

  • 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 Einführung.

    [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; Einführung in die Theorie der Markov-Ketten.

    [9] Spellecchia, Claudia, 2008; MCMC.

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

    21

  • 9 Abbildungsverzeichnis

    Abbildung 1.1: Verschlüsselter Text;Diaconis, Persi; The Markov Chain Monte Carlo Revolution.

    Abbildung 1.2: Entschlüsselter 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