MCMC-Methoden - Markov Chain Monte Carlo

Inhaltsverzeichnis

1  Grundlagen zu Markov-Ketten  4

1.1  Definition  4

1.2  Irreduzibilit at und Aperiodizit at  9

1.3  Station are Verteilungen und Reversibilit at  11

1.3.1  Existenz der station aren Verteilung  11

1.3.2  Reversibilit at einer Verteilung  18

1.4  Konvergenzsatz  19

1.4.1  Konvergenz gegen die station are Verteilung  19

1.4.2  Eindeutigkeit der station aren Verteilung  21

2  Metropolis-Hastings Algorithmus  23

2.1  Allgemeine Beschreibung des Metropolis-Hastings Algorithmus  23

2.2  Implementierung des Metropolis-Algorithmus im Beispiel der Exponential  

  verteilung  26

2.3  Fehler-Absch atzung im Beispiel der Exponetialverteilung  28

3  Gibbs-Sampler  30

3.1  Allgemeine Beschreibung des Gibbs-Samplers  30

3.2  Implementierung des Gibbs-Sampler Beispiels  33

3.3  Verallgemeinerung auf q-F arbungen  34

4  Approximate counting  36

4.1  Problemstellung  36

4.2  Existenz-Theorem  37

4.3  Beweis: erster Teil  38

4.4  Beweis: zweiter Teil  41

4.5  Implementierung  45

5  Literatur  49

6  Quelltexte  50

6.1  Metropolis-Hastings Algorithmus  50

  2  

6.2 Gibbs-Sampler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

6.3 Approximate-Counting Algorithmus . . . . . . . . . . . . . . . . . . . . . . 52

3

1 Grundlagen zu Markov-Ketten

1.1 Definition

Wir beginnen mit einem sehr einfachen Beispiel: Denken wir an einen zuf¨ alligen L¨ aufer in einer sehr kleinen Stadt, die nur aus vier Straßen besteht. Dabei werden die vier Stra- ßenecken wie in der untenstehenden Abbildung mit v 1 , v 2 , v 3 und v 4 bezeichnet. Zum Zeitpunkt 0 steht der zuf¨ allige L¨ aufer in der Ecke v 1 . Zum Zeitpunkt 1 wirft er eine faire M¨ unze und entscheidet je nach Ausfall, ob er weiter nach v 2 oder v 4 geht. Zum Zeitpunkt

2 wirft er wieder eine faire M¨ unze, um zu entscheiden, zu welcher benachbarten Ecke er

gehen soll. Dabei verwendet er die Entscheidungsregel, wenn die M¨ unze Kopf zeigt, einen Schritt im Uhrzeigersinn zu gehen und andernfalls, wenn die M¨ unze Zahl zeigt, einen Schritt gegen den Uhrzeigersinn zu gehen. Diese Prozedur wird fortgef¨ uhrt f¨ ur die Zeiten

3, 4, usw.

F¨ ur jedes n bezeichnet X n den Index der Straßenecke an der sich der L¨ aufer zur Zeit n befindet. Deswegen ist (X 0 , X 1 , ...) ein Zufallsprozess der Werte aus {v 1 , v 2 , v 3 , v 4 } an- nimmt. Weil der L¨ aufer zur Zeit 0 in v 1 startet, ergibt sich:

P (X 0 = v 1 ) = 1

Danach bewegt er sich nach v 2 oder v 4 mit jeweils der Wahrscheinlichkeit 1

2

1

4

und

1

2

Die Verteilung von X n f¨ ur n ≥ 2 zu berechnen, setzt ein wenig mehr Nachdenken voraus. An dieser Stelle ist es sinnvoll, sich auf bedingte Wahrscheinlichkeiten zu beziehen. Nehmen wir an, dass der L¨ aufer zur Zeit n an der Ecke v 2 steht. Dann erhalten wir die bedingten Wahrscheinlichkeiten

1

P (X n+1 = v 1 |X n = v 2 ) =

2

und

1

wegen des M¨ unzwurf-Mechanismus zur Entscheidung ¨ uber den n¨ achsten Schritt. Tats¨ achlich erhalten wir dieselben bedingten Wahrscheinlichkeiten, wenn wir in den Bedingungen die vollst¨ andige Vergangenheit des Prozesses bis zur Zeit n ber¨ ucksichtigen:

1

P (X n+1 = v 1 |X 0 = i 0 , X 1 = i 1 , ..., X n−1 = i n−1 , X n = v 2 ) =

2

und

1

P (X n+1 = v 3 |X 0 = i 0 , X 1 = i 1 , ..., X n−1 = i n−1 , X n = v 2 ) =

2

Dies gilt f¨ ur jede Wahl von i 0 , ..., i n−1 , vorausgesetzt, dass der Pfad i 0 , i 1 , ..., i n−1 ei- ne positive Wahrscheinlichkeit besitzt. Dieses Ph¨ anomen wird die Ged¨ achnislosigkeits- Eigenschaft genannt, die auch als Markov-Eigenschaft bekannt ist: Die bedingte Verteilung von X n+1 gegeben (X 0 , ..., X n ) h¨ angt nur von X n ab. Eine andere interessante Eigenschaft dieses Zufallsprozesses besteht darin, dass die bedingte Verteilung von X n+1 , gegeben dass beispielsweise X n = v 2 , f¨ ur alle n dieselbe ist. Diese Eigenschaft heißt Zeit-Homogenit¨ at, oder einfach Homogenit¨ at. Derartige Prozesse lassen sich mit Hilfe sogenannter ¨ Uber- gangsmatrizen beschreiben:

Definition Sei P eine k × k Matrix mit den Elementen {P i,j : i, j = 1, ..., k}. Ein Zu- fallsprozess (X 0 , X 1 , ...) mit endlichem Zustandsraum S = {s 1 , ..., s k } heißt dann (homo- gene) Markovkette mit ¨

Ubergangsmatrix P , wenn f¨ ur alle n, alle i, j ∈ {1, ..., k} und alle i 0 , ..., i n−1 ∈ {1, ..., k} gilt:

5

P (X n+1 = s j |X 0 = s i 0 , X 1 = s i 1 , ..., X n−1 = s i n−1 , X n = s i )

= P (X n+1 = s j |X n = s i ) = P i,j

Zum Beispiel ist das Beispiel mit dem zuf¨ alligen L¨ aufer von oben eine Markovkette, mit

dem Ereignisraum {1, ..., 4} und der ¨ Ubergangsmatrix

 

 

Jede ¨

Ubergangsmatrix erf¨ ullt dabei:

P i,j ≥ 0 ∀i, j ∈ {1, ..., k}

und

k

Die erste Eigenschaft bedeutet nur, dass alle bedingten Wahrscheinlichkeiten stets nicht-

negativ sind und die zweite Eigenschaft besagt, dass sie sich zu 1 aufsummieren.

Wir wenden uns nun einer anderen wichtigen Charakteristik von Markovketten zu,

n¨ amlich der sogenannten Anfangsverteilung, die beschreibt, wie die Markovkette startet.

Die Anfangsverteilung ist gegeben durch einen Zeilenvektor µ (0) , der sich folgendermaßen

zusammensetzt:

(0) (0) µ (0) = (µ 1 , µ 2 , ..., µ

= (P (X 0 = s 1 ), P (X 0 = s 2 ), ..., P (X 0 = s k )).

Da µ (0) eine Wahrscheinlichkeitsverteilung repr¨ asentiert, ergibt sich:

k

In unserem Beispiel mit dem zuf¨ alligen L¨ aufer gilt:

6

µ (0) = (1, 0, 0, 0)

Auf ¨ ahnliche Weise bezeichen die Zeilenvektoren µ (1) , µ (2) , ... die Verteilungen der Mar-

kovkette zu den Zeiten 1, 2, ..., so dass:

(n) (n) µ (n) = (µ 1 , µ 2 , ..., µ

= (P (X n = s 1 ), P (X n = s 2 ), ..., P (X n = s k ))

Im Beispiel mit dem zuf¨ alligen L¨ aufer ergibt sich beispielsweise:

1 1

Wenn wir die Anfangsverteilung µ (0) und die ¨

kennen, k¨ onnen wir alle Verteilungen µ (1) , µ (2) , ... der Markovkette berechnen: Das folgende

Resultat zeigt uns, dass dies lediglich eine Anwendung der Matrix-Multiplikation darstellt.

Wir schreiben P n f¨ ur die n-te Potenz der Matrix P .

Theorem F¨ ur eine Markovkette (X 0 , X 1 ...) mit dem Zustandsraum {s 1 , ..., s k }, Anfangs-

verteilung µ (0) und ¨

Ubergangsmatrix P , gilt f¨ ur jedes n, dass die Verteilung µ (n) zur Zeit

n Folgendes erf¨ ullt:

µ (n) = µ (0) P n

Beweis Wenden wir uns dem ersten Fall, n = 1 zu. F¨ ur j = 1, ..., k erhalten wir: Aus

dem Satz ¨ uber die totale Wahrscheinlichkeit folgt, dass

k

(1)

µ j = P (X 1 = s j ) =

denn die s i liefern eine disjunkte Zerlegung des Ereignisraumes S = {s 1 , ..., s k }. Nach

Definition der bedingten Wahrscheinlichkeit folgt, dass

k

(1)

µ j =

7

Nach Definition von µ und P gilt somit k

oder anders geschrieben

(1)

Hierzu ist ein wenig Wissen ¨ uber Matrix-Arithmetik n¨ otig. Schreibt man die Multiplika- tion nach dem Falk-Schema 1 aus, ergibt sich:  

Da wir die G¨ ultigkeit der Gleichung f¨ ur jede Komponente j einzeln nachweisen konnten, gilt also insgesamt:

µ (1) = µ (0) P Um diesen Sachverhalt f¨ ur den generellen Fall zu beweisen, benutzen wir vollst¨ andige Induktion. Wir fixieren ein m und nehmen an, dass die Behauptung f¨ ur n = m gilt. F¨ ur n = m + 1 erhalten wir dann analog zum Induktionsanfang: k

(m+1)

µ = P (X m+1 = s j ) = j

k

P (X m = s i )P (X m+1 = s j |X m = s i )

= i=1

k

=

Da wir auch f¨ ur den allgemeinen Fall die G¨ ultigkeit der Gleichung f¨ ur jede Komponente j einzeln nachweisen konnten, folgt insgesamt:

µ (m+1) = µ (m) P = µ (0) P m P = µ (0) P m+1 1 Siehe hierzu beispielsweise Seite 144 in: Gerd Fischer: Lineare Algebra: Eine Einf¨ uhrung f¨ ur Studi- enanf¨ anger, Vieweg Verlag

8

Das erste Gleichheitszeichen wurde zuvor bewiesen. Das zweite Gleichheitszeichen gilt nach Induktionsvoraussetzung.

1.2 Irreduzibilit¨ at und Aperiodizit¨ at

Wir beginnen mit der Irreduzibilit¨ at, was salopp gesprochen bedeutet, dass jedes Element s i des Ereignisraums S der Markovkette von jedem anderen Element s j aus erreicht wer- den kann. Um diese Beschreibung pr¨ aziser zu machen, stellen wir uns eine Markovkette (X 0 , X 1 , ...) mit Ereignisraum S = {s 1 , ..., s k } und ¨

Ubergangsmatrix P vor. Wir sagen, dass das Ereignis s i mit s j kommuniziert, geschrieben s i → s j , wenn die Kette eine positi- ve Wahrscheinlichkeit hat, irgendwann s j zu erreichen, wenn sie von s i startet. In anderen Worten: s i kommuniziert mit s j , wenn ein n existiert, so dass gilt:

P (X m+n = s j |X m = s i ) > 0

Diese Wahrscheinlichkeit l¨ asst sich auch schreiben als (P n ) i,j . Wenn sogar s i → s j und s j → s i gilt, sagen wir, dass s i und s j interkommunizieren und schreiben daf¨ ur s i ↔ s j . Dies leitet uns zu folgender Definition der Irreduzibilit¨ at:

Definition Eine Markovkette (X 0 , X 1 , ...) mit Ereignisraum S = {s 1 , ..., s k } und ¨

Uber- gangsmatrix P heißt irreduzibel, wenn f¨ ur alle s i , s j ∈ S gilt, dass s i ↔ s j . Ansonsten heißt die Kette reduzibel.

Wir kommen nun zum Konzept der Aperiodizit¨ at. F¨ ur eine endliche oder unendli- che Menge {a 1 , a 2 , ...} von positiven ganzen Zahlen, schreiben wir gcd{a 1 , a 2 , ...} f¨ ur den gr¨ oßten gemeinsamen Teiler von a 1 , a 2 , .... Die Periode d(s i ) von einem Zustand s i ∈ S ist nun definiert als:

d(s i ) = gcd{n ≥ 1 : (P n ) i,i > 0}, d.h. die Periode von s i ist der gr¨ oßte gemeinsame Teiler der Menge von Zeiten, zu denen die Kette nach s i zur¨ uckkehren kann, gegeben dass wir bei s i gestartet sind 2 . Wenn d(s i ) = 1 gilt, dann sagen wir, dass der Zustand s i aperiodisch ist.

2 Diese Menge k¨ onnte f¨ ur eine ung¨ unstige Wahl der Markovkette auch leer sein. Diesen Sonderfall be- trachten wir nicht als aperiodisch.

9

Definition Eine Markovkette heißt aperiodisch, wenn alle Elemente ihres Zustandsraums aperiodisch sind. Andernfalls heißt die Kette periodisch.

Theorem Angenommen, wir haben eine aperiodische Markovkette (X 0 , X 1 , ...) mit dem Ereignisraum S = {s 1 , ..., s k } und der ¨

Ubergangsmatrix P . Dann existiert stets ein N ∈ N, so dass (P n ) i,i > 0 ∀n ≥ N ∀i ∈ {1, ..., k}

Lemma Sei A = {a 1 , a 2 , ...} eine Menge von positiven ganzen Zahlen, f¨ ur die gilt:

1. Die Elemente von A sind teilerfremd, d.h. gcd{a 1 , a 2 , ...} = 1, und

2. A ist abgeschlossen unter der Addition, d.h. wenn a ∈ A und a ′ ∈ A, gilt a + a ′ ∈ A.

Dann existiert stets eine ganze Zahl N ∈ N, so dass n ∈ A ∀n ≥ N 3 .

Beweis des Theorems F¨ ur s i ∈ S, sei A i = {n ≥ 1 : (P n ) i,i > 0}, so dass A i in anderen Worten die Menge der m¨ oglichen R¨ uckkehr-Zeitpunkte zu s i darstellt, wenn von s i gestar- tet wird. Wir haben angenommen, das die Markovkette aperiodisch ist und deswegen ist auch das Ereignis s i aperiodisch, so dass die Elemente von A i teilerfremd sind. Außerdem ist A i abgeschlossen unter der Addition, aus dem folgenden Grund: Wenn a, a ′ ∈ A i , dann gilt P (X a = s i |X 0 = s i ) > 0 und P (X a+a ′ = s i |X a = s i ) > 0. Dies bedeutet aber dass

P (X a+a ′ = s i |X 0 = s i ) ≥ P (X a = s i , X a+a ′ = s i |X 0 = s i ),

P (X a+a ′ = s i , X a = s i , X 0 = s i )

= ,

P (X a+a ′ = s i |X a = s i )P (X a = s i |X 0 = s i )P (X 0 = s i )

= , = P (X a+a ′ = s i |X a = s i )P (X a = s i |X 0 = s i ) > 0,

da beide Faktoren nach Voraussetzung gr¨ oßer als 0 sind. Insgesamt folgt nun, dass a + a ′ ∈ A i . Zusammenfassend erf¨ ullt nun A i die Voraussetzungen 1 und 2 des Lemmas, was bedeutet, dass eine ganze Zahl N i < ∞ existiert, so dass n ∈ A i ∀n ≥ N i . Dies bedeutet

3 Ein Beweis dieses Lemmas befindet sich in: Br´ emaud, P. (1998) Markov Chains: Gibbs fields, Monte Carlo Simulation, and Queues, Springer, New York.

10

aber, dass (P n ) i,i > 0 ∀n ≥ N i . Das Theorem folgt nun, indem wir N folgendermaßen w¨ ahlen:

N = max{N 1 , ..., N k }

Korollar Sei (X 0 , X 1 , ...) eine irreduzible und aperiodische Markovkette mit Ereignis-

raum S = {s 1 , ..., s k } und ¨ Ubergangsmatrix P . Dann existiert stets ein M ∈ N, so dass

(P n ) i,j > 0 ∀n ≥ M ∀i, j ∈ {1, ..., k}.

Beweis Wegen der vorausgesetzten Aperiodizit¨ at existiert eine ganze Zahl N ∈ N, so

dass (P n ) i,i > 0 ∀n ≥ N ∀i ∈ {1, ..., k}. Fixiere nun zwei Zust¨ ande s i , s j ∈ S. Wegen der

vorausgesetzten Irreduzibilit¨ at k¨ onnen wir nun ein n i,j finden, so dass (P n i,j ) i,j > 0 gilt. Sei nun M i,j = N + n i,j . F¨ ur jedes m ≥ M i,j ergibt sich nun:

P (X m = s j |X 0 = s i ) ≥ P (X m−n i,j = s i , X m = s j |X 0 = s i ),

P (X m = s j , X m−n i,j = s i , X 0 = s i )

= ,

P (X m = s j |X m−n i,j = s i )P (X m−n i,j = s i |X 0 = s i )P (X 0 = s i )

= , = P (X m = s j |X m−n i,j = s i )P (X m−n i,j = s i |X 0 = s i ) > 0, wobei der erste Faktor wegen der geeigneten Wahl von n i,j gr¨ oßer 0 ist und der zweite

Faktor positiv ist, weil m − n i,j ≥ N . Gezeigt wurde also, dass (P m ) i,j > 0 ∀m ≥ M i,j . Das Korollar folgt nun, wenn wir M folgendermaßen festsetzen:

M = max{M 1,1 , M 1,2 , ..., M 1,k , M 2,1 , M 2,2 , ..., M 2,k , ..., M k,k }

1.3 Station¨ are Verteilungen und Reversibilit¨ at

1.3.1 Existenz der station¨ aren Verteilung

In diesem Abschnitt geht es um einige wichtige Aspekte der Markov-Ketten: Asymptoti-

sches Verhalten bei der Langzeit-Betrachtung von Markov-Ketten. Was k¨ onnen wir also

uber eine Markov-Kette sagen, die schon lange Zeit gelaufen ist? Wenn (X 0 , X 1 , ...) eine ¨

nichttriviale 4 Markov-Kette ist, dann wird der Wert von X n unendlich oft fluktuieren

4 Wir gehen hier von den ¨ ublicherweise auftretenden Markov-Ketten aus und betrachten nicht Beson-

derheiten wie absorbierende Zust¨ ande.

11