Metropolis-Hastings-Algoritmus

28.09.2005, 18:08

Michael78

Metropolis-Hastings-Algoritmus

Hallo!
Ich kämpfe mich gerade mit dem Metropolis-Hastings-Algoritmus ab - kann mir jemand dazu vielleicht ein numerisches Beispiel geben?

Was mir nicht klar ist, ist diese Akzeptanzwahrscheinlichkeit:
Da ist ja erst ein Wert X(t) und dann wird ein darauf folgender Wert X(t+1) gezogen, der überprüft wird.
Zur Überprüfung wird dann die Akzeptanzwahrscheinlichkeit berechnet, wobei zwei Dichtefunktionen zur Anwendung kommen:
So weit ich das verstanden habe, die Prioriwahrscheinlichkeit und eine so genannte "Vorschlagsdichte" (proposal density).
Ich verstehe das nicht: Ist es bei diesen Bayes-Sachen nicht so, dass die Prioridichte bereits eine vorgeschlagene Dichte ist, der sich die Posterioridichte dann annähert?
Worin unterscheiden sich denn in dem Algoritmus die zwei Dichten??

Ich wäre jedem herzlich dankbar, der mir hier ein wirklich idiotensicheres Beispiel Schritt für Schritt zeigen würde.
Wikipedia hat mir leider nicht weitergeholfen bzw. brauche ich wohl ein Beispiel...

Vielen Dank und schönen Abend den hilfreichen Experten! :-)

28.09.2005, 18:20

Auf diesen Beitrag antworten »

Von welchen Dichten sprichst du? Von denen der Markov-Kette $\begin{eqnarray*} X(t) \end{eqnarray*}$ , oder von der der Vorschlagsverteilung? Das sind bei MH ganz verschiedene Dinge. Also pack mal paar Symbole und Bezeichnungen hier rein, auf die man sich dann beziehen kann - ohne das geht es bei dieser komplexen Materie nicht.

29.09.2005, 09:06

Michael78

Auf diesen Beitrag antworten »

a=a1 * a2

a1 = P(x')/P(xt)
a2 = Q(xt;x')/Q(x'/xt)

Siehe http://en.wikipedia.org/wiki/Metropolis-Hastings_algorithm

Mir ist da einfach nicht klar:
P(x) ist ja die Funktion, die approximiert werden soll.
Wird dann für die Vorschlagsdichte Q(x';xt) eine andere Verteilung herangezogen?

Ich wäre für ein Zahlenbeispiel wirklich dankbar.

29.09.2005, 10:13

Auf diesen Beitrag antworten »

Das Problem ist, dass ein einigermaßen sinnvolles Beispiel schon enorm komplex ist. Ich wähle daher ein eher "unsinniges" Beispiel einer Simulation, d.h., eine die man in der Praxis nie und nimmer so durchführen würde:

Wir wollen mit MH die Standardnormalverteilung $\begin{eqnarray*} \mathcal{N}(0,1) \end{eqnarray*}$ simulieren, d.h., reelle Zufallsgrößen mit der Dichte

$\begin{eqnarray*} f(x)=\frac{1}{\sqrt{2\pi}}\cdot \exp\left\{ -\frac{x^2}{2} \right\} \end{eqnarray*}$

Dann ist gemäß MH z.B. folgendes Vorgehen möglich: Von einem Startwert $\begin{eqnarray*} x_0 \end{eqnarray*}$ aus konstruiert man eine Folge $\begin{eqnarray*} (x_n)_{n=1,2,\ldots} \end{eqnarray*}$ in der Weise, dass man am Punkt $\begin{eqnarray*} x=x_n \end{eqnarray*}$ angelangt einen Vorschlagswert $\begin{eqnarray*} x'=x+u \end{eqnarray*}$ bestimmt, wobei $\begin{eqnarray*} u \end{eqnarray*}$ eine im Intervall $\begin{eqnarray*} [-h,h] \end{eqnarray*}$ liegende gleichverteilte Zufallsgröße sein möge. In Symbolen ausgedrückt entspricht das der Vorschlagsverteilungsdichte

$\begin{eqnarray*} q(x,x') = \begin{cases} \frac{1}{2h} & \;\text{für}\; x-h<x'<x+h\\ 0 & \;\text{sonst} \end{cases} \end{eqnarray*}$

Der MH-Algorithmus geht nun so, dass dieses $\begin{eqnarray*} x' \end{eqnarray*}$ mit Wahrscheinlichkeit

$\begin{eqnarray*} p = \min\left\{ 1, \frac{f(x')q(x,x')}{f(x)q(x',x)} \right\} \end{eqnarray*}$

als neuer Wert der Folge, also $\begin{eqnarray*} x_{n+1}:=x' \end{eqnarray*}$ akzeptiert werden soll. Andernfalls bleibt alles beim alten, d.h., $\begin{eqnarray*} x_{n+1}:=x (=x_n) \end{eqnarray*}$ .

Berechnen wir also mal diese Akzeptanzwahrscheinlichkeit für unser konkretes Beispiel: Hier liegt nun gerade Symmetrie $\begin{eqnarray*} q(x,x')=q(x',x) \end{eqnarray*}$ vor, das ist der sogenannte Metropolis-Algorithmus (ein Spezialfall, der historisch eher entwickelt wurde). Also erhalten wir als Akzeptanzwahrscheinlichkeit

$\begin{eqnarray*} p = \min\left\{ 1, \frac{f(x')}{f(x)} \right\} = \min\left\{ 1, \exp\left\{ \frac{x^2-x'^2}{2} \right\} \right\} = \begin{cases} 1 & \;\text{für}\; |x'|\leq |x|\\ \exp\left\{ \frac{x^2-x'^2}{2} \right\} & \;\text{sonst} \end{cases} \end{eqnarray*}$

Die entstehende Folge $\begin{eqnarray*} (x_n) \end{eqnarray*}$ liefert nun im stationären Fall (also nach einer Einschwingzeit, englisch "burn-in", d.h. für $\begin{eqnarray*} n\geq n_0 \end{eqnarray*}$ ) tatsächlich standardnormalverteilte Zufallsgrößen, wobei allerdings aufeinander folgende Glieder dieser Folge nicht unabhängig sind.

Spannend ist die Frage, wie diese Vorschlagsverteilung zu wählen ist. Das ist in sehr weiten Grenzen variabel - es muss nur garantiert sein, dass in endlich vielen Schritten der Wertebereich der Zufallsgröße "ausgelotet" werden kann, d.h., alle "Ecken" müssen erreichbar sein. Aber natürlich ist nicht jede Vorschlagsverteilung gleich gut:

Nehmen wir unser Beispiel, da ist ja noch das $\begin{eqnarray*} h \end{eqnarray*}$ variabel. Wählen wir ein sehr kleines $\begin{eqnarray*} h \end{eqnarray*}$ , z.B. $\begin{eqnarray*} h=0.1 \end{eqnarray*}$ , dann erreichen wir sehr hohe Akzeptanzraten, bewegen uns aber von Schritt zu Schritt kaum vom Fleck. Wählen wir andererseits ein sehr großes $\begin{eqnarray*} h \end{eqnarray*}$ , z.B. $\begin{eqnarray*} h=10 \end{eqnarray*}$ , dann springen wir mit den Vorschlägen x' sehr oft aus dem "wahrscheinlichen" Bereich (3Sigma-Regel!) der Normalverteilung heraus, was sich in einer geringen Akzeptanzwahrscheinlichkeit äußert - d.h. die Folge "steht" regelrecht eine ganze Weile, bis sie von Zeit zu Zeit einen größeren "Satz" macht.

Richtig verstehen kannst du das vermutlich erst, wenn du das mal konkret simulierst, und dir eine durche Simulatuon entstandene Folge $\begin{eqnarray*} (x_n) \end{eqnarray*}$ anschaust, am besten im Plot.

29.09.2005, 10:33

Michael78

Auf diesen Beitrag antworten »

Herzlichen Dank für die ganze Mühe und die wirklich detaillierte Erklärung!!

Ich glaube, ich bin jetzt wirklich einen sehr großen Schritt weiter.

Doch noch ein paar Fragen:
- Auf das x folgt ja in der Kette immer ein neu gezogenes x' oder wieder das x. Was bedeutet hier die Wahrscheinlichkeit für "sonst":
Wird die neue Zahl dann summiert aus x*(1-a) + x' * a?

- Woran merkt man, dass die Einschwingphase vorüber ist?
- Wenn die nachfolgenden Zahlen voneinander nicht unabhängig sind, liegt das daran, dass sie weiterhin dem "Akzeptanz-Test" unterliegen?
Oder braucht es diesen dann gar nicht mehr, weil er nur die Konvergenz sichern soll?

29.09.2005, 11:10

Auf diesen Beitrag antworten »

Zitat:

Original von Michael78
- Auf das x folgt ja in der Kette immer ein neu gezogenes x' oder wieder das x. Was bedeutet hier die Wahrscheinlichkeit für "sonst":
Wird die neue Zahl dann summiert aus x*(1-a) + x' * a?

Nein - das ist eine direkte Wahrscheinlichkeitsdichte, besser wäre vielleicht die Schreibweise $\begin{eqnarray*} q_x(x') \end{eqnarray*}$ , um deutlich zu kennzeichnen, was das Argument ist. Das "sonst" kennzeichnet den Fall, dass $\begin{eqnarray*} x' \end{eqnarray*}$ außerhalb des Intervalls $\begin{eqnarray*} (x-h,x+h) \end{eqnarray*}$ liegen möge - der Dichtewert 0 in diesem Fall zeigt allerdings, dass das nicht passieren soll.

D.h., dieses $\begin{eqnarray*} q(x,x')=q_x(x') \end{eqnarray*}$ ist nur eine anderere Schreibweise für das vorher geschriebene $\begin{eqnarray*} x'=x+u \end{eqnarray*}$ . Simulieren würde man die Berechnung von $\begin{eqnarray*} x' \end{eqnarray*}$ so, dass man ein $\begin{eqnarray*} u\sim \mathrm{Glm}[-h,h] \end{eqnarray*}$ auswürfelt und zu dem "alten" $\begin{eqnarray*} x \end{eqnarray*}$ addiert. Ich habe das $\begin{eqnarray*} q(x,x') \end{eqnarray*}$ nur deswegen aufgeschrieben, um diese Simulation richtig in die MH-Schreibweise einordnen zu können.

Tatsächlich verwendet man MH eher bei viel komplizierteren bzw. höherdimensionalen Verteilungen, und nicht wie hier bei so einer schnöden eindimensionalen Geschichte. Erstmalig aufgetaucht ist der Metropolis-Algorithmus meines Wissens nach in der Arbeit

N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. . Teller, and E. Teller. Equations of state calculations by fast computing machines. J. Chem. Phys., 21:1087--1091, 1953.

(Edward Teller dürfte ein Begriff sein - der Vater der Wasserstoffbombe), und zwar bei Teilchensystem-Simulationen. Hier besteht also ein Zustand aus mehreren Hundert oder Tausend Dimensionen (!), nämlich den Positionen (dreidimensional) und evtl. sogar Geschwindigkeiten aller Teilchen des Systems!!!

Ein Übergang $\begin{eqnarray*} x\to x' \end{eqnarray*}$ ändert dann nicht die Position aller Teilchen, sondern nur eines einzigen! Für solche Anwendungsfälle ist der MH bestens gerüstet.

Zitat:

Original von Michael78
- Woran merkt man, dass die Einschwingphase vorüber ist?

Das ist die Frage, die in den wunden Punkt tippt - das kann man adhoc nicht sagen. Dazu muss man sich Zeitreihen gewisser Charakteristiken (hier in diesem einfachen Beispiel die Folge selbst; gleitende Mittelwerte; Autokorrelation; ...) anschauen und beurteilen, wann die Einschwingphase vorbei ist. Da ist gerade bei komplexen Simulationen eine Menge Heuristik dabei.

Zitat:

Original von Michael78
- Wenn die nachfolgenden Zahlen voneinander nicht unabhängig sind, liegt das daran, dass sie weiterhin dem "Akzeptanz-Test" unterliegen?
Oder braucht es diesen dann gar nicht mehr, weil er nur die Konvergenz sichern soll?

Die Frage ist, was du mit den simulierten Werten anstellen willst. Für viele Sachen stört die Abhängigkeit gar nicht, wenn man z.B. den Erwartungswert $\begin{eqnarray*} E(g(X)) \end{eqnarray*}$ irgendeines Funktionals $\begin{eqnarray*} g \end{eqnarray*}$ über dieser Verteilung simulieren will. Braucht man wirklich unabhängige (oder besser gesagt "nahezu unabhängige") Werte, dann muss man sich die Autokorrelation dieser Folge genauer unter die Lupe nehmen. Ist die nach, sagen wir $\begin{eqnarray*} k \end{eqnarray*}$ Schritten hinreichend niedrig (<10^(-6) o.ä.), dann kann man die Folgenglieder $\begin{eqnarray*} x_{n_0},x_{n_0+k},x_{n_0+2k},\ldots \end{eqnarray*}$ als unabhängig ansehen, d.h., man nimmt nur jeden $\begin{eqnarray*} k \end{eqnarray*}$ -ten Wert der Folge und wirft die anderen weg.

29.09.2005, 11:19

Michael78

Auf diesen Beitrag antworten »

Entschuldigung, bei der ersten Frage war ich ungenau:

Mit "sonst" habe ich nicht die 0-Variante der Vorschlagsdichte gemeint, sondern weiter unten in Deinem Beispiel das "sonst" in der Formel der Akzeptanzwahrscheinlichkeit.
Ich habe gelesen, dass entweder x oder x' verwendet wird, was bedeutet dann zum Beispiel eine Akzeptanzwahrscheinlichkeit von 30%:
Dass der nachfolgende Wert in der Kette errechnet wird durch x*0,7+x' *o,3 ?

Und meine (hoffentlich, versprochen smile

) letzte Frage:
wenn also das Burn-in abgeschlossen ist, hört man dann auch mit den Akzeptant-Tests auf, oder?
Und falls ja: Wodurch entsteht dann die Abhängigkeit?

Danke!

29.09.2005, 11:25

Auf diesen Beitrag antworten »

Zitat:

Original von Arthur Dent
Der MH-Algorithmus geht nun so, dass dieses $\begin{eqnarray*} x' \end{eqnarray*}$ mit Wahrscheinlichkeit $\begin{eqnarray*} p \end{eqnarray*}$ als neuer Wert der Folge, also $\begin{eqnarray*} x_{n+1}:=x' \end{eqnarray*}$ akzeptiert werden soll. Andernfalls bleibt alles beim alten, d.h., $\begin{eqnarray*} x_{n+1}:=x (=x_n) \end{eqnarray*}$ .

Eindeutiger geht's nicht. Also keine Linearkombinationen, sondern du musst "würfeln":

* Du simulierst einen in [0,1] gleichverteilten Wert z.
* Ist z<p, dann setzt du $\begin{eqnarray*} x_{n+1}:=x' \end{eqnarray*}$ , andernfalls setzt du $\begin{eqnarray*} x_{n+1}:=x_n \end{eqnarray*}$ .

Alles klar?

29.09.2005, 11:36

Michael78

Auf diesen Beitrag antworten »

Klar, danke!
Irgendwie kam mir Würfeln so.. zufällig vor ;-)

Letzte Frage daher: Wenn also das Burn-in abgeschlossen ist, hört man dann auch mit den Akzeptanz-Tests auf, oder?
Falls ja: Wodurch entsteht dann die Abhängigkeit?

29.09.2005, 11:48

Auf diesen Beitrag antworten »

Nein, das Vorgehen bleibt bestehen!

Schau dir das doch nochmal genau an: In der Vorschlagsverteilung im Beispiel steckt doch keinerlei Information über die zu simulierende Normalverteilung, also kannst du nicht erwarten, dass bei "Weglassen" der Akzeptanz-Prozedur diese Normalverteilung simuliert wird!!!

Du musst mal den Grundgedanken von MH, oder allgemein vom übergeordneten "Markov Chain Monte Carlo" (MCMC - denn dazu gehört ja MH) verstehen: Zu einer Verteilung $\begin{eqnarray*} P \end{eqnarray*}$ , z.B. mit gegebener Dichte $\begin{eqnarray*} f \end{eqnarray*}$ wie hier, wird eine Markov-Kette $\begin{eqnarray*} (X_n) \end{eqnarray*}$ gerade so konstruiert, dass sie als stationäre Verteilung genau dieses $\begin{eqnarray*} P \end{eqnarray*}$ besitzt.

Und zur Abhängigkeit: Schau dir nochmal mein Beispiel mit dem Vielteilchensystem an:

Zitat:

Original von Arthur Dent
Ein Übergang $\begin{eqnarray*} x\to x' \end{eqnarray*}$ ändert dann nicht die Position aller Teilchen, sondern nur eines einzigen!

Da dürfte doch offensichtlich sein, wieso $\begin{eqnarray*} x \end{eqnarray*}$ und $\begin{eqnarray*} x' \end{eqnarray*}$ abhängig sind: Wenn du von 1000 Teilchen nur die Position eines einzigen Teilchens änderst, dann ist das System ja noch "fast" dasselbe!!!

05.10.2005, 11:10

Michael78

Auf diesen Beitrag antworten »

Danke!
Noch eine Frage zu den komplizierteren Fällen - ich bin leider kein Physiker, aber danke für den Literaturhinweis (Metropolis et al. 1953) - :

Hier habe ich gelesen, dass man blockweise vorgehen kann, also zB einmal 5 Variablen updatet und dann wieder 5 usw.
Mit ist dabei nicht ganz klar:
Wenn ich 5 Variablen update, heißt das also, dass ich 5 mal aus der gleichen Verteilung ziehe und für jede den Akzeptanz-Test mache.

1.) Stimmt das so?
2.) Ändert sich dann für die nächsten 5 irgendwas zB an der Zielverteilung/Posterioriverteilung?
3) Falls nein, wo liegt dann der Vorteil des blockweisen Vorgehens gegenüber dem einzelnen Vorgehen?

----------------------------------------- Z U S A M M E N G E F Ü G T -----------------------------------------

Vielleicht damit ihr seht, wovon ich als Input ausgehe:

http://www2.gsb.columbia.edu/divisions/m..._metroplois.pdf

(Wesentlich schwieriger und umfangreicher ist:
http://www.princeton.edu/~yacine/eco575/...ib-shephard.pdf
)

Nur zur Einordnung, ich studiere Wirtschaft.

\\edit by mercany: Doppelpost zusammengefügt. Bitte benutze die edit-Funktion!

Vielleicht damit ihr seht, wovon ich als Input ausgehe:

http://www2.gsb.columbia.edu/divisions/m..._metroplois.pdf

(Wesentlich schwieriger und umfangreicher ist:
http://www.princeton.edu/~yacine/eco575/...ib-shephard.pdf
)

Nur zur Einordnung, ich studiere Wirtschaft.

edit: Zweiten Doppelpost zusammengefügt, bitte benutze die edit-Funktion! (MSS)

05.10.2005, 11:49

Auf diesen Beitrag antworten »

Also den Metropolis-Artikel habe ich mehr aus historischen Gründen angeführt, ich habe ihn auch hier momentan nicht vorliegen, also kann ich zu Details des Inhalts wenig sagen. Zum Verständnis von MH (bzw. allgemeiner MCMC) sind andere Artikel oder Bücher sicher viel besser geeignet.

Außer Historikern bzw. historisch interessierten Mathematikern liest heute ja auch keiner mehr die Original-Arbeiten von Leonard Euler... Augenzwinkern

05.10.2005, 12:06

Michael78

Auf diesen Beitrag antworten »

Das erleichtert mich :-)

Trotzdem checke ich im Moment noch nicht, wie das im Falle mehrerer Variablen geht, die in einem einzigen Iterationsschritt gesamplet werden. (?)

UND - Wenn ich das Princeton-Dokument genannt habe:
Da werden außerdem noch unbeobachtbare Parameter gesamplet, ich blicke dabei gar nicht mehr durch, wo denn die beobachteten Daten eigentlich ins Spiel kommen...

(vielleicht sollte ich meinen Nick kurzerhand in "verwirrt" umtaufen, schlimm ist das)

Neue Frage »

Antworten »

Metropolis-Hastings-Algoritmus

Verwandte Themen