Hinweise zu den Programmierprojekten

Abgabe bis 17.6. im CIP-Pool. Bitte vereinbaren Sie rechtzeitig einen Termin.

Schreiben Sie Ihre Lösungen bitte direkt in die Notebooks mit den Projektbeschreibungen (download auf der Vorlesungshomepage), oder kopieren Sie die Aufgabenstellungen zu den Lösungen.

Jedes der Projekte zählt vierfach. Wir empfehlen die Bearbeitung von einem der beiden Projekte - natürlich können Sie auch beide abgeben.

Projekt 1: Simulation des Isingmodells

Einleitung

Wir betrachten ein großes Rechteck

({{({{Λ={0,1,…,2d+1}×}}){0,1,...,2d+1}}})

im (ℤ^2), und den Konfigurationsraum

({{({{S = (({0,1})^Λ) = {x=(Subscript[((Subscript[x, i,j])), (i,j)∈Λ]): (Subscript[x, i,j])∈{0,1} für (i,j)∈Λ}}})}})

In Mathematica können wir ein Element aus S als zweidimensionales Feld von Nullen und Einsen darstellen, und mithilfe von ArrayPlot graphisch darstellen. Eine zufällige (gleichverteilte) Konfiguration aus S erhält man beispielsweise mit:

d=20; x=Table[RandomInteger[],{i,0,2d+1},{j,0,2d+1}]; ArrayPlot[x]

Wir wollen nun Wahrscheinlichkeitsverteilungen auf S betrachten, bei denen im Gegensatz zur Gleichverteilung benachbarte Pixel vorzugsweise im selben Zustand sind. Dazu definieren wir die "Energie" einer Konfiguration x als

({{({{H(x) =  ∑((|(Subscript[x, i,j])-(Subscript[x, k,l])|)^2)}}) ,}})

wobei die Summe über alle i,j,k,l=0,1,...,2d+1 mit ‖(i,j)-(k,l)‖=1 läuft. Jedes Paar von benachbarten Pixeln wird also zweimal berücksichtigt. Man kann H(x) auch als ein Maß für die Größe der Grenze zwischen der schwarzen und der weißen Region in der Konfiguration x interpretieren. Wir definieren nun Wahrscheinlichkeitsverteilungen, unter denen Konfigurationen mit kleiner Energie bevorzugt werden. Dazu beschränken wir uns auf Konfigurationen, für die alle Pixel am Rand des Rechtecks schwarz gefärbt sind. Sei

({{({{(SuperPlus[S])={x∈S}}) : (Subscript[x, i,j])=1  für alle (i,j) mit i=0 oder i=2d+1 oder j=0 oder j=2d+1} .}})

Für β≥0 betrachten wir die Wahrscheinlichkeitsverteilung (Subscript[μ, β]) auf (SuperPlus[S]) mit Massenfunktion

({{({{(Subscript[μ, β])(x) = (1/(Subscript[Z, β]))(e^(-β·H(x)))}})    ,       x∈(SuperPlus[S]),}})

wobei

({{({{(Subscript[Z, β]) = (∑_(x∈(SuperPlus[S])))(e^(-β·H(x)))}})}})

eine Konstante ist, die (Subscript[μ, β])  zur Massenfunktion einer Wahrscheinlichkeitsverteilung normiert. Je größer der Wert von β ist, umso stärker werden Konfigurationen mit vielen gleichfarbigen Nachbarpixeln unter der Verteilung (Subscript[μ, β])  bevorzugt. Insbesondere ist (Subscript[μ, 0]) die Gleichverteilung auf (SuperPlus[S]), während für β→∞  die Wahrscheinlichkeit, daß alle Pixel wie der Rand schwarz sind, gegen 1 konvergiert. Die folgenden Grafiken zeigen Stichproben von (Subscript[μ, β])  für d=20 und β=0.2, 0.44 bzw. 0.6 :

Hintergrund

Das beschriebene Modell ist eines der fundamentalen Modelle der statistischen Physik, und wurde von Ising 1925 eingeführt. Heute spielt das Ising-Modell auch eine wichtige Rolle in vielen anderen Bereichen, zum Beispiel als A-priori-Verteilung in Bayesschen Modellen der Bildverarbeitung. In der physikalischen Interpretation beschreibt das Modell (auf stark vereinfachte Weise) ein ferromagnetisches Material. An jeder Position des Gitters Λ befindet sich ein Atom, dessen Spin nach oben (+1) oder nach unten (0) zeigen kann. Der Parameter β ist die inverse Temperatur. Bei hoher Temperatur sind die Spins ungeordnet, bei niedriger Temperatur zeigen Sie vorzugsweise in eine Richtung. 
Sei (Subscript[m, β])(d) der Erwartungswert unter (Subscript[μ, β]) des Wertes (Subscript[x, d,d]) in der Mitte des Rechtecks Λ. Ein tiefes Resultat von Onsager besagt, daß

({{({{(lim_(d→∞))(Subscript[m, β])(d)=(1/2)                                                                                   falls  sinh(2β)≤1, während
(lim_(d→∞))(Subscript[m, β])(d)=(1/2)+(1/2)(((1-sinh(((2β))^(-4))))^(1/8)) >(1/2)            falls sinh(2β)>1.}})}})

Im Fall sinh(2β)>1 beeinflussen also die Randwerte den Wert in der Mitte sogar im Limes d→∞, während für kleinere β der Einfluss der Randwerte auf den Mittelpunkt für d→∞ verschwindet. Im physikalischen Modell bedeutet dies, daß bei der kritischen inversen Temperatur

({{(Subscript[({{β}}), crit]) = (1/2)ArcSinh(1) ≈ 0.44}})

ein Phasenübergang stattfindet.

Simulation des Ising-Modells

Ein einfaches Acceptance-Rejection-Verfahren ist zur Simulation des Isingmodells schon für relativ kleine Werte von d nicht mehr praktikabel. Man verwendet daher meist Markov Chain Monte Carlo Verfahren. Tatsächlich ist die beim Gibbs-Sampler simulierte Markovkette auf dem Raum (SuperPlus[S]) sehr ähnlich zu der physikalischen Glauberdynamik des Spinsystems, und wurde ursprünglich auch von dieser motiviert. Um einen Gibbs-Sampler effizient implementieren zu können, unterteilen wir das Gitter Λ (ohne Rand) schachbrettartig in gerade und ungerade Felder:

({{({{(Λ^g) = {(i,j): i,j=1,2,...,2d   mit  i+j  gerade},
(Λ^u)  = {(i,j): i,j=1,2,...,2d   mit  i+j  ungerade},}})}})

und setzen

({{(({{x}})^g)=((Subscript[x, i,j])| (i,j)∈(Λ^g)),          (x^u) = ((Subscript[x, i,j])| (i,j)∈(Λ^u)).}})

Da die Konfigurationen x∈(SuperPlus[S]) auf dem Rand von Λ überall den Wert 1 annehmen, können wir x identifizieren mit dem Paar ((x^g),(x^u)). Zudem haben die bedingten Verteilungen von (x^g) gegeben den Wert von (x^u) und umgekehrt eine einfache Form. Bezeichnet man mit

({{(Subscript[({{s}}), i,j]) = (Subscript[x, i-1,j])+(Subscript[x, i+1,j])+(Subscript[x, i,j-1])+(Subscript[x, i,j+1])}})

die Anzahl der schwarzen Pixel in der Umgebung von (i,j), dann gilt

(Subscript[μ, β])((x^g)|(x^u))=(∏_((i,j)∈(Λ^g)))(Subscript[ν, β])((Subscript[x, i,j])|(Subscript[s, i,j])),        und

(Subscript[μ, β])((x^u)|(x^g))=(∏_((i,j)∈(Λ^u)))(Subscript[ν, β])((Subscript[x, i,j])|(Subscript[s, i,j])),

wobei (Subscript[ν, β])(·|s) die durch

({{({{                         (Subscript[ν, β])(1|s) = (e^(4β·(s-2)))}})/(1+(e^(4β·(s-2)))) ,}})({{({{(Subscript[ν, β])(0|s) = 1}})/(1+(e^(4β·(s-2)))) ,}})

definierte Bernoulliverteilung auf {0,1} bezeichnet. Unter der bedingten Verteilung gegeben (x^u) sind die (Subscript[x, i,j]),  (i,j)∈(Λ^g), also unabhängig mit Verteilung (Subscript[ν, β])(·|(Subscript[s, i,j])), und analog sind unter der bedingten Verteilung gegeben (x^g) die (Subscript[x, i,j]),  (i,j)∈(Λ^u), unabhängig mit Verteilung (Subscript[ν, β])(·|(Subscript[s, i,j])). Beide bedingten Verteilungen können wir also leicht simulieren. Damit liegt folgender Algorithmus zur Simulation von (Subscript[μ, β]) nahe:

Algorithmus (Gibbs-Sampler für das zweidimensionale Isingmodell)

Initialisierung:

Erzeuge (Subscript[x, i,j]^((0)))  (1≤i,j≤2d) unabhängig, =0,1 jeweils mit W'keit 1/2.

Setze (Subscript[x, 0,j]^((0)))=(Subscript[x, 2d+1,j]^((0)))=(Subscript[x, i,0]^((0)))=(Subscript[x, i,2d+1]^((0)))=1  (0≤i,j≤2d+1).

For n:=1,2,... do

x:=(x^((n-1)))

Update (x^g)=((Subscript[x, i,j]): i,j∈(Λ^g))  ~ (Subscript[μ, β])((x^g)|(x^u))

Update (x^u)=((Subscript[x, i,j]): i,j∈(Λ^u))  ~ (Subscript[μ, β])((x^u)|(x^g))

(x^((n))) :=x

Die Wahl der Anfangsbedingung (x^((0))) ist dabei relativ willkürlich  - man kann auch andere Anfangsbedingungen verwenden.

Theoretische Resultate

Die in dem Algorithmus simulierte Markovkette (X^((n))) hat die Verteilung (Subscript[μ, β]) als Gleichgewicht. Aus dem Konvergenzsatz für Markovketten folgt, daß (X^((n))) für große n näherungsweise gemäß (Subscript[μ, β]) verteilt ist - also ist (x^((n))) eine approximative Stichprobe von (Subscript[μ, β]). 
Zudem folgt aus dem Ergodensatz (Gesetz der großen Zahlen) für Markovketten, daß für jede Funktion f:(SuperPlus[S])→ℝ und für jedes b∈ℕ die empirischen Mittelwerte

({{(Subscript[(OverHat[({{ϑ}})]), n]) := (1/n) (∑_(k=b)^(b+n-1))f((x^((k))))}})

für n→∞ mit Wahrscheinlichkeit 1 gegen den Erwartungswert

({{({{ϑ   := }}) ∑f(x)(Subscript[μ, β])(x)}})

von f unter (Subscript[μ, β]) konvergieren. Durch hinreichend große Wahl von b kann man die Konvergenz beschleunigen, s.u.

Aufgabenstellung

Zeigen Sie, daß die bedingte Verteilung

({{({{(Subscript[μ, β])((x^g)|(x^u)) := (((Subscript[μ, β])((x^g),(x^u)))/((∑_(x^g)) (Subscript[μ, β])((x^g),(x^u))))}})}})

durch (1) gegeben ist, und folgern Sie, daß (Subscript[μ, β]) ein Gleichgewicht für die durch obigen Algorithmus definierte Markovkette (X^((n))) ist.  (Hinweis: Da auf beiden Seiten von Gleichung (1) Wahrscheinlichkeitsverteilungen stehen, genügt es zu zeigen, daß

({{({{(Subscript[μ, β])((x^g)|(x^u))  ∝ (Subscript[μ, β])((x^g),(x^u))  ∝ }})(∏_((i,j)∈(Λ^g))) exp(4β·(Subscript[x, i,j])·((Subscript[s, i,j])-2))  ∝ (∏_((i,j)∈(Λ^g))) (Subscript[ν, β])((Subscript[x, i,j])|(Subscript[s, i,j]))}})

gilt, wobei "∝" bedeutet, daß beide Seiten als Funktionen in (x^g) proportional sind - sich also nur durch eine von (x^u) abhängende multiplikative Konstante unterscheiden).

Implementieren Sie den oben beschriebenen Gibbs-Sampler in Mathematica. Gehen Sie dazu folgendermaßen vor:

Als Vorübung: Erzeugen Sie ein Schachbrettmuster (d.h. die Konfiguration x mit (x^g)≡1 und (x^u)≡0), und erstellen Sie einen ArrayPlot.

Erstellen Sie Mathematica-Funktionen updateg[x_] und updateu[x_], die die Werte (Subscript[x, i,j]) für (i,j)∈(Λ^g) bzw.  (i,j)∈(Λ^u) gemäß den jeweiligen bedingten Verteilungen updaten. Hierbei können β und d als Konstanten verwendet werden, denen vor Aufruf der Funktionen ein fester Wert zugewiesen wird. Wenden Sie jeweils einen Update auf das Schachbrettmuster an, erstellen Sie einen ArrayPlot der erhaltenen Konfigurationen für verschiedene Werte von β und d, und überprüfen Sie, ob Sie die erwarteten Ergebnisse erhalten.

Erzeugen Sie nun (z.B. mit NestList[ ] ) eine Liste von ArrayPlots für die ersten 20 Schritte des Algorithmus mit d=20 und β=0.2, 0.44 bzw. 0.6. Führen Sie anschließend 300 Schritte durch, und stellen Sie den Verlauf der Konfigurationen in diesen 300 Schritten mithilfe von Manipulate[ ] graphisch dar. Die Berechnung kann nun etwas Zeit in Anspruch nehmen - immerhin müssen 300 mal 1600 Pixel erneuert werden. Was beobachten Sie ? Für welchen Wert von β dauert es am längsten, bis sich ein Gleichgewicht einstellt ?

Experimentieren Sie auch mit anderen Randbedingungen, z.B. (Subscript[x, 0,j])=(Subscript[x, i,0])=0, (Subscript[x, 2d+1,j])=(Subscript[x, i,2d+1])=1 für alle i,j.

Verwenden Sie den Gibbs-Sampler, um den Wert von (Subscript[m, β])(d) zu schätzen. Gehen Sie dabei wie folgt vor:

Schreiben Sie eine Mathematica-Funktion MCMCSchaetzer[b_,nmax_], die die Schätzer
({{({{(Subscript[(OverHat[m]), β]^((n)))}}) := (1/n) (∑_(k=b)^(b+n-1))(Subscript[x, d,d]^((k)))          (n=1,2,...,nmax)}})
berechnet. Wählen Sie b für die obigen Werte von d und β jeweils so groß, daß Sie den Eindruck haben, daß sich nach b Schritten in guter Näherung ein Gleichgewicht eingestellt hat ("burn-in-time"). Stellen Sie den Verlauf der Schätzer als Funktion von n graphisch dar, wobei Sie nmax möglichst groß wählen.

Erzeugen Sie eine Graphik, in der der Verlauf von (Subscript[m, β])(d) als Funktion von β approximiert wird. Schätzen Sie dazu (Subscript[m, β])(d) für beispielsweise 40 verschiedene Werte von β zwischen 0.2 und 0.6, und interpolieren Sie die erhaltenen Schätzwerte linear.