Gibbs-Sampler

Der in Abschnitt 3.3.1 betrachtete MCMC-Algorithmus, um ,,zufällig ausgewählte'' zulässige Konfigurationen des Hard-Core-Modells zu erzeugen, ist ein Spezialfall des sogenannnten Gibbs-Samplers zur MCMC-Simulation von diskreten (hochdimensionalen) Zufallsvektoren.

• Sei $V$ eine endliche (nichtleere) Index-Menge, und sei ${\mathbf{X}}=(X(v),\,v\in V)$ ein diskreter Zufallsvektor,
  • der mit Wahrscheinlichkeit $1$ seine Werte in der endlichen (Zustands-) Menge $E\subset\mathbb{R}^{\vert V\vert}$ annimmt, wobei wir voraussetzen,
  • dass es für jedes Paar ${\mathbf{x}},{\mathbf{x}}^\prime\in E$ eine (endliche) Folge von Zuständen ${\mathbf{y}}_0,{\mathbf{y}}_1,\ldots,{\mathbf{y}}_n\in E$ gibt, so dass

    $${\mathbf{y}}_0={\mathbf{x}}\,,\quad {\mathbf{y}}_n={\mathbf{x}}^\prime \qquad \char93 \{v\in V:\,y_i(v)\not=y_{i+1}(v)\}=1\,,\quad\forall\, i=0,\ldots,n-1\,.$$ (34)

    
    

• Sei ${\boldsymbol{\pi}}=(\pi_{\mathbf{x}},\,{\mathbf{x}}\in E)$ die Wahrscheinlichkeitsfunktion des Zufallsvektors ${\mathbf{X}}$ mit $\pi_{\mathbf{x}}>0$ für jedes ${\mathbf{x}}\in E$, und für jedes $v\in V$ sei

  $$\pi_{x(v)\mid\, {\mathbf{x}}(-v)}=P\bigl(X(v)=x(v)\mid\, {\mathbf{X}}(-v)={\mathbf{x}}(-v)\bigr)$$ (35)

  
  
  • die bedingte Wahrscheinlichkeit, dass die Komponente $X(v)$ von ${\mathbf{X}}$ den Wert $x(v)$ annimmt,
  • unter der Bedingung, dass der Vektor ${\mathbf{X}}(-v)=(X(w),\,w\in V\setminus\{v\})$ der übrigen Komponenten gleich ${\mathbf{x}}(-v)$ ist, wobei vorausgesetzt wird, dass $(x(v),{\mathbf{x}}(-v))\in E$.

MCMC-Simulationsalgorithmus

 

• Ähnlich wie in Abschnitt 3.3.1 konstruieren wir eine Markov-Kette ${\mathbf{X}}_0,{\mathbf{X}}_1,\ldots$
  • mit dem Zustandsraum $E$ und mit einer (geeignet gewählten) irreduziblen und aperiodischen Übergangsmatrix ${\mathbf{P}}$,
  • so dass ${\boldsymbol{\pi}}$ die ergodische (Grenz-) Verteilung der zu konstruierenden Markov-Kette ${\mathbf{X}}_0,{\mathbf{X}}_1,\ldots$ ist.
• Dabei erzeugen wir erneut mit Hilfe der rekursiven Konstruktion, die bereits in Abschnitt 2.1.3 diskutiert wurde, einen ,,Pfad'' ${\mathbf{x}}_0,{\mathbf{x}}_1,\ldots$ der Markov-Kette:

  
  

  1. Wähle einen Anfangszustand ${\mathbf{x}}_0\in E$.
  2. Wähle eine Komponente $v\in V$ gemäß einer (vorgegebenen) Wahrscheinlichkeitsfunktion ${\mathbf{q}}=(q_v,\,v\in V)$, so dass $q_v>0$ für jedes $v\in V$.
  3. Generiere den (aktualisierten) Wert $x_{n+1}(v)$ der $v$-ten Komponente gemäß der (bedingten) Wahrscheinlichkeitsfunktion
     $${\boldsymbol{\pi}}_{\cdot\,\vert{\mathbf{x}}_n(-v)}=\bigl(\pi_{x(v)\mid\, {\mathbf{x}}_n(-v)},\, \forall\, x(v)\;$$$$\mbox{mit $(x(v),{\mathbf{x}}(-v))\in E$}$$$$\bigr)\,.$$

  4. Die Werte sämtlicher Komponenten $w\not= v$ bleiben unverändert, d.h., es gilt $x_{n+1}(w)=x_n(w)$ für jedes $w\not= v$.

Theorem 3.11   $\;$ Die Übergangsmatrix ${\mathbf{P}}=(p_{{\mathbf{x}}{\mathbf{x}}^\prime})$ sei gegeben durch

$$p_{{\mathbf{x}}{\mathbf{x}}^\prime}=\sum\limits_{v\in V} q_v \pi_... ...x}}^\prime(-v)\bigr)\,,\qquad\forall\, {\mathbf{x}},{\mathbf{x}}^\prime\in E\,,$$ (36)

wobei die bedingten Wahrscheinlichkeiten $\pi_{x^\prime(v)\mid\, {\mathbf{x}}(-v)}$ in $(\ref{def.deb.gib})$ gegeben sind. Dann ist ${\mathbf{P}}$ irreduzibel und aperiodisch, und das Paar $({\mathbf{P}},{\boldsymbol{\pi}})$ ist reversibel.

Beweis

$\;$ Die Behauptung ergibt sich auf ähnliche Weise wie im Beweis von Theorem 3.10.

• Um die Aperiodizität von ${\mathbf{P}}=(p_{{\mathbf{x}}{\mathbf{x}}^\prime})$ zu zeigen, genügt es zu beachten,
  • dass für jedes ${\mathbf{x}}\in E$
    $$p_{{\mathbf{x}},{\mathbf{x}}} = \sum\limits_{v\in V} q_v\pi_{x(v)... ..._{{\mathbf{z}}\in E:\,{\mathbf{z}}(-v)={\mathbf{x}}(-v)} \pi_{\mathbf{z}}}\;>0$$

  • und dass somit sämtliche Diagonalelemente $p_{{\mathbf{x}},{\mathbf{x}}}$ von ${\mathbf{P}}$ positiv sind.
• Die Irreduzibilität von ${\mathbf{P}}$ ergibt sich aus den folgenden Überlegungen.
  • Für beliebige, jedoch fest vorgegebene ${\mathbf{x}},{\mathbf{x}}^\prime\in E$ sei $k\le\vert V\vert$ die Anzahl derjenigen Komponenten $v\in V$ mit $x(v)\not=x^\prime(v)$.
  • Für $k=0$, d.h. ${\mathbf{x}}={\mathbf{x}}^\prime$, hatten wir bereits beim Nachweis der Aperiodizität gezeigt, dass $p_{{\mathbf{x}}{\mathbf{x}}}>0$.
  • Sei nun $k>0$, wobei wir o.B.d.A. annehmen, dass die Komponenten linear durchnumeriert sind und dass die ersten $k$ Komponenten von ${\mathbf{x}}$ und ${\mathbf{x}}^\prime$ voneinander verschieden sind.
  • Die Voraussetzung (34) über die Struktur des Zustandsraumes $E$ impliziert, dass es eine Folge ${\mathbf{y}}_0,\ldots,{\mathbf{y}}_k\in E$ gibt, so dass ${\mathbf{y}}_0={\mathbf{x}}$ und
    $${\mathbf{y}}_1=\bigl(x^\prime(v_1),x(v_2),\ldots,x(v_{\vert V\ver... ...\prime(v_k),x(v_{k+1}),\ldots,x(v_{\vert V\vert})\bigr)={\mathbf{x}}^\prime\,.$$

  • Für jedes $i=0,\ldots,k-1$ gilt dann

    $$p_{{\mathbf{y}}_i{\mathbf{y}}_{i+1}}= q_{v_i} \pi_{y_{i+1}(v_i)\m... ...athbf{z}}\in E:\,{\mathbf{z}}(-v_i)={\mathbf{y}}_i(-v_i)} \pi_{\mathbf{z}}}\;>0$$ (37)

    
    
    und somit auch $p^{(k)}_{{\mathbf{x}}{\mathbf{x}}^\prime}\ge\prod_{i=0}^{k-1} p_{{\mathbf{y}}_i{\mathbf{y}}_{i+1}}>0$.
• Es ist nun noch die Gültigkeit der Detailed-Balance-Gleichung (2.85) zu zeigen, d.h., dass

  $$\pi_{\mathbf{x}}\,p_{{\mathbf{x}}{\mathbf{x}}^\prime}=\pi_{{\math... ...\prime{\mathbf{x}}}\,, \qquad\forall\, {\mathbf{x}},{\mathbf{x}}^\prime\in E\,.$$ (38)

  
  
  • Wenn ${\mathbf{x}}={\mathbf{x}}^\prime$ gilt, dann ist die Gültigkeit von (38) offensichtlich.
  • Wenn $x(v)\not=x^\prime(v)$ für mehr als eine Komponente $v\in V$, dann gilt $p_{{\mathbf{x}}{\mathbf{x}}^\prime}=p_{{\mathbf{x}}^\prime{\mathbf{x}}}=0$, und (38) gilt somit auch in diesem Fall.
  • Es gelte nun $x(v)\not=x^\prime(v)$ für ein $v\in V$ und $x(w)=x^\prime(w)$ für alle $w\not= v$. Dann ergibt sich aus (37), dass
    $$\pi_{\mathbf{x}}\,p_{{\mathbf{x}}{\mathbf{x}}^\prime}=\pi_{\mathb... ...hbf{z}}}\; = \pi_{{\mathbf{x}}^\prime}\,p_{{\mathbf{x}}^\prime{\mathbf{x}}}\,.$$
    
     
      $\Box$
    
    
    

Sei ${\mathbf{X}}_0,{\mathbf{X}}_1,\ldots$ eine Markov-Kette mit dem Zustandsraum $E$ und mit der in (36) gegebenen Übergangsmatrix ${\mathbf{P}}=(p_{{\mathbf{x}}{\mathbf{x}}^\prime})$. Aus Theorem 3.11 ergibt sich dann, dass

$$\lim\limits_{n\to\infty} d_{\rm TV}({\boldsymbol{\alpha}}_n,{\boldsymbol{\pi}})=0$$ (39)

für jede Anfangsverteilung ${\boldsymbol{\alpha}}_0$, wobei ${\boldsymbol{\alpha}}_n$ die Verteilung von ${\mathbf{X}}_n$ bezeichnet. Außerdem gilt für den Gibbs-Sampler die folgende Monotonie-Eigenschaft.

Theorem 3.12   $\;$ Für jedes $n=0,1,\ldots$ gilt

$$d_{\rm TV}({\boldsymbol{\alpha}}_n,{\boldsymbol{\pi}})\ge d_{\rm TV}({\boldsymbol{\alpha}}_{n+1},{\boldsymbol{\pi}})\,.$$ (40)

Beweis

 

• Für beliebige $v\in V$ und ${\mathbf{x}}^\prime\in E$ ergibt sich aus (35), dass
  

  $$\sum\limits_{{\mathbf{x}}\in E:\,{\mathbf{x}}(-v)={\mathbf{x}}^\prime(-v)} \pi_{x^\prime(v)\mid\, {\mathbf{x}}(-v)}\;\pi_{\mathbf{x}} \stackrel{(\ref{def.deb.gib})}{=} \sum\limits_{{\mathbf{x}}\in E :\,{\mathbf{x}}(-v)={\mathbf{x}}^\prime(-v) } \pi_{{\mathbf{x}}^\prime} \pi_{x(v)\mid\, {\mathbf{x}}^\prime(-v)}$$

  $$= \pi_{{\mathbf{x}}^\prime} \underbrace{\sum\limits_{{\mathbf{x}}\i... ...)} \pi_{x(v)\mid\, {\mathbf{x}}^\prime(-v)}}_{=1} =\pi_{{\mathbf{x}}^\prime}\,.$$ (41)

  
  

• Hieraus und aus der Definitionsgleichung (36) der Übergangsmatrix ${\mathbf{P}}=(p_{{\mathbf{x}}{\mathbf{x}}^\prime})$ ergibt sich, dass
  

  $$2\;d_{\rm TV}({\boldsymbol{\alpha}}_{n+1},{\boldsymbol{\pi}}) = \sum\limits_{{\mathbf{x}}^\prime\in E} \bigl\vert\,\alpha_{n+1,\,{\mathbf{x}}^\prime}-\pi_{{\mathbf{x}}^\prime}\bigr\vert$$

  $$= \sum\limits_{{\mathbf{x}}^\prime\in E}\Bigl\vert\,\sum\limits_{{\... ...f{x}}}\,p_{{\mathbf{x}}{\mathbf{x}}^\prime}-\pi_{{\mathbf{x}}^\prime}\Bigr\vert$$

  $$\stackrel{(\ref{ein.ube.gib})}{=} \sum\limits_{{\mathbf{x}}^\prime\in E}\Bigl\vert\,\sum\limits_{{\... ...athbf{x}}(-v)={\mathbf{x}}^\prime(-v)\bigr)-\pi_{{\mathbf{x}}^\prime}\Bigr\vert$$

  $$= \sum\limits_{{\mathbf{x}}^\prime\in E}\Bigl\vert\sum\limits_{v\in... ...{\mathbf{x}}(-v)}\alpha_{n,\,{\mathbf{x}}} -\pi_{{\mathbf{x}}^\prime}\Bigr\vert$$

  $$\stackrel{(\ref{hil.equ.for})}{=} \sum\limits_{{\mathbf{x}}^\prime\in E}\Bigl\vert\sum\limits_{v\in... ...bf{x}}(-v)} \bigl(\alpha_{n,\,{\mathbf{x}}} -\pi_{{\mathbf{x}}}\bigr)\Bigr\vert$$

  $$\le \sum\limits_{{\mathbf{x}}^\prime\in E}\sum\limits_{v\in V} q_v\,\... ...bf{x}}(-v)} \bigl\vert\,\alpha_{n,\,{\mathbf{x}}} -\pi_{{\mathbf{x}}}\bigr\vert$$

  $$= \sum\limits_{{\mathbf{x}}\in E}\unde... ...^\prime}=1} \bigl\vert\,\alpha_{n,\,{\mathbf{x}}} -\pi_{{\mathbf{x}}}\bigr\vert$$

  $$= \sum\limits_{{\mathbf{x}}\in E} \bigl\vert\,\alpha_{n,\,{\mathbf{... ...x}}}\bigr\vert\;=\;2\;d_{\rm TV}({\boldsymbol{\alpha}}_n,{\boldsymbol{\pi}})\,.$$

  
  
  
   
    $\Box$
  
  
  

Beachte

 

• Eine modifizierte Version des bisher in diesem Abschnitt betrachteten Gibbs-Samplers ist der zyklische Gibbs-Sampler, bei dem die jeweils zu aktualisierende Komponente $v\in V$
  • nicht gemäß einer (vorgegebenen) Wahrscheinlichkeitsfunktion ${\mathbf{q}}=(q_v,\,v\in V)$, so dass $q_v>0$ für jedes $v\in V$, ausgewählt wird,
  • sondern die Komponenten $v\in V$ werden durchnumeriert und der Reihe nach (auf deterministische Weise) ausgewählt.
• Falls $k=n\vert V\vert+i$ für gewisse Zahlen $n=0,1,\ldots$ und $i=1,\ldots,\vert V\vert$, dann ist die Matrix ${\mathbf{P}}(k)=\bigl(p_{{\mathbf{x}}{\mathbf{x}}^\prime}(k)\bigr)$ der Übergangswahrscheinlichkeiten $p_{{\mathbf{x}}{\mathbf{x}}^\prime}(k)$ im $k$-ten Schritt gegeben durch

  $$p_{{\mathbf{x}}{\mathbf{x}}^\prime}(k)=\pi_{x^\prime(v_i)\mid\, {... ...}^\prime(-v_i)\bigr)\,,\qquad\forall\, {\mathbf{x}},{\mathbf{x}}^\prime\in E\,.$$ (42)

  
  

• Für einen gesamten (Scan-) Zyklus, bei dem jede Komponente genau einmal aktualisiert wird, ergibt sich die Übergangsmatrix

  $${\mathbf{P}}={\mathbf{P}}(1)\cdot\ldots\cdot{\mathbf{P}}(\vert V\vert)\,.$$ (43)

  
  

• Man kann sich leicht überlegen, dass die in (42) und (43) gegebene Matrix ${\mathbf{P}}=(p_{{\mathbf{x}}{\mathbf{x}}^\prime})$
  • irreduzibel und aperiodisch ist
  • und dass ${\boldsymbol{\pi}}$ die stationäre (Grenz-) Verteilung von ${\mathbf{P}}$ ist,
  • denn für jedes $i=1,\ldots,\vert V\vert$ und für jedes ${\mathbf{x}}^\prime\in E$ ergibt sich aus (41) und (42), dass
    $$\sum\limits_{{\mathbf{x}}\in E} \pi_... ...\prime(-v_i)\bigr) \stackrel{(\ref{hil.equ.for})}{=} \pi_{{\mathbf{x}}^\prime}$$

  • und dass somit auch
    $$\sum\limits_{{\mathbf{x}}\in E} \pi_{\mathbf{x}}p_{{\mathbf{x}}{\mathbf{x}}^\prime}= \pi_{{\mathbf{x}}^\prime}\,.$$

• Das Paar $({\mathbf{P}},{\boldsymbol{\pi}})$ ist im allgemeinen nicht reversibel. In Abschnitt 2.3.4 hatten wir jedoch gezeigt, dass das Paar $({\mathbf{M}},{\boldsymbol{\pi}})$ reversibel ist, wobei

  $${\mathbf{M}}={\mathbf{P}}\widetilde{\mathbf{P}} \qquad \widetilde{\mathbf{P}}={\,{\rm diag}}(\pi_{\mathbf{x}}^{-1}){\mathbf{P}}^\top{\,{\rm diag}}(\pi_{\mathbf{x}})$$ (44)

  
  
  die multiplikativ reversible Version von ${\mathbf{P}}$ bezeichnet.

Theorem 3.13   $\;$ Es gilt

$${\mathbf{M}}={\mathbf{P}}(1)\cdot\ldots\cdot{\mathbf{P}}(\vert V\vert)\cdot{\mathbf{P}}(\vert V\vert)\cdot\ldots\cdot{\mathbf{P}}(1)\,,$$ (45)

d.h., die multiplikativ reversible Version ${\mathbf{M}}$ der ,,Vorwärts-Scan-Matrix'' ${\mathbf{P}}$ stimmt mit der ,,Vorwärts-Rückwärts-Scan-Matrix'' überein.

Beweis

 

• Es genügt zu zeigen, dass $\widetilde{\mathbf{P}}={\mathbf{P}}(\vert V\vert)\cdot\ldots\cdot{\mathbf{P}}(1)$ für die in (44) definierte Matirx $\widetilde{\mathbf{P}}=(\widetilde p_{{\mathbf{x}}{\mathbf{x}}^\prime})$ gilt.
• Aus (42)-(44) ergibt sich, dass für beliebige ${\mathbf{x}},{\mathbf{x}}^\prime\in E$
  

  $${\widetilde{\mathbf{p}}_{{\mathbf{x}}{\mathbf{x}}^\prime} = \Bigg... ...^\top{\,{\rm diag}}(\pi_{\mathbf{x}})\Biggr)_{{\mathbf{x}}{\mathbf{x}}^\prime}}$$

  $$= \Biggl( {\,{\rm diag}}(\pi_{\mathbf{x}}^{-1}){\mathbf{P}}^\top(\v... ...op(1) {\,{\rm diag}}(\pi_{\mathbf{x}})\Biggr)_{{\mathbf{x}}{\mathbf{x}}^\prime}$$

  $$= \sum\limits_{{\mathbf{y}}_1,\ldots,{\mathbf{y}}_{\vert V\vert-1}\... ...\bigr)\pi_{y_1(v_{\vert V\vert-1})\mid\, {\mathbf{y}}_2(-v_{\vert V\vert-1})}\,$$

  $$\hspace{1cm}\times\,{1\hspace{-1mm}{\rm I}}\bigl({\mathbf{y}}_1(-... ...rt V\vert-1}(-v_1)={\mathbf{x}}^\prime(-v_1)\bigr) \pi_{{\mathbf{x}}^\prime}\,.$$

  
  

• Hieraus und aus (35) ergibt sich nun (ähnlich wie beim Beweis von (38)), dass
  

  $${\widetilde{\mathbf{p}}_{{\mathbf{x}}{\mathbf{x}}^\prime} = \sum\... ...)\bigr)\pi_{y_2(v_{\vert V\vert-1})\mid\, {\mathbf{y}}_1(-v_{\vert V\vert-1})}}$$

  $$\hspace{1cm}\times\,{1\hspace{-1mm}{\rm I}}\bigl({\mathbf{y}}_1(-... ...rm I}}\bigl({\mathbf{y}}_{\vert V\vert-1}(-v_1)={\mathbf{x}}^\prime(-v_1)\bigr)$$

  $$= \Bigl({\mathbf{P}}(\vert V\vert)\cdot\ldots\cdot{\mathbf{P}}(1)\Bigr)_{{\mathbf{x}}{\mathbf{x}}^\prime}\,.$$

  
  
  
   
    $\Box$
  
  
  

Beachte

 

• Bei der praktischen Anwendung von Gibbs-Samplern wird stets vorausgesetzt,
  • dass die in (36) bzw. (42) betrachteten bedingten Wahrscheinlichkeiten
    $$\pi_{x(v)\mid\, {\mathbf{x}}(-v)}=P\bigl(X(v)=x(v)\mid\, {\mathbf{X}}(-v)={\mathbf{x}}(-v)\bigr)$$
    für jedes $v\in V$ lediglich von dem Vektor $\bigl(x(w),\,w\in\mathcal{N}(v)\bigr)$ derjenigen Werte abhängen,
  • die der Zufallsvektor ${\mathbf{X}}=(X(w),\,w\in V)$ in einer gewissen (möglichst kleinen) Nachbarschaft $\mathcal{N}(v)\subset V$ von $v\in V$ annimmt.
• Dabei wird eine Familie $\mathcal{N}=\{\mathcal{N}(v),\,v\in V\}$ von Teilmengen von $V$ ein System von Nachbarschaften genannt, falls für beliebige $v,w\in V$

  (a)

  $v\not\in \mathcal{N}(v)$,

  (b)

  aus $w\in\mathcal{N}(v)$ stets $v\in\mathcal{N}(w)$ folgt.

• Für das in Abschnitt 3.3.1 diskutierte Beispiel des Hard-Core-Modells ist $\mathcal{N}(v)$ für jeden Eckpunkt $v\in V$ die Menge derjenigen Eckpunkte $w\not= v$, die (direkt) mit $v$ durch eine Kante verbunden sind.