Gibbs Sampler

The MCMC algorithm for the generation of ,,randomly picked'' admissible configurations of the hard core model (see Section 3.3.1) is a special case of a so-called Gibbs sampler for the simulation of discrete (high dimensional) random vectors.

• Let $V$ be a finite (nonempty) index set and let ${\mathbf{X}}=(X(v),\,v\in V)$ be a discrete random vector
  • taking values in the finite state space $E\subset\mathbb{R}^{\vert V\vert}$ with probability $1$ where we assume
  • that for every pair ${\mathbf{x}},{\mathbf{x}}^\prime\in E$ there is a finite sequence of states ${\mathbf{y}}_0,{\mathbf{y}}_1,\ldots,{\mathbf{y}}_n\in E$ such that

    $${\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)

    
    

• Let ${\boldsymbol{\pi}}=(\pi_{\mathbf{x}},\,{\mathbf{x}}\in E)$ be the probability function of the random vector ${\mathbf{X}}$ with $\pi_{\mathbf{x}}>0$ for all ${\mathbf{x}}\in E$, and for all $v\in V$ let

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

  
  
  • denote the conditional probability that the component $X(v)$ of ${\mathbf{X}}$ has the value $x(v)$
  • given that the vector ${\mathbf{X}}(-v)=(X(w),\,w\in V\setminus\{v\})$ of the other components equals ${\mathbf{x}}(-v)$ where we assume $(x(v),{\mathbf{x}}(-v))\in E$.

MCMC Simulation Algorithm

 

• Similar to Section 3.3.1 we construct a Markov chain ${\mathbf{X}}_0,{\mathbf{X}}_1,\ldots$
  • with state space $E$ and an (appropriately chosen) irreducible and aperiodic transition matrix ${\mathbf{P}}$,
  • such that ${\boldsymbol{\pi}}$ is the ergodic limit distribution of ${\mathbf{X}}_0,{\mathbf{X}}_1,\ldots$.
• Then we generate a ,,path'' ${\mathbf{x}}_0,{\mathbf{x}}_1,\ldots$ of the Markov chain by the recursive construction discussed in Section 2.1.3:

  
  

  1. Pick an initial state ${\mathbf{x}}_0\in E$.
  2. Pick a component $v\in V$ according to a given probability function ${\mathbf{q}}=(q_v,\,v\in V)$ such that $q_v>0$ for all $v\in V$.
  3. Generate the update $x_{n+1}(v)$ of the $v$th component according to the (conditional) probability function
     $${\boldsymbol{\pi}}_{\cdot\,\vert{\mathbf{x}}_n(-v)}=\bigl(\pi_{x(v)\mid\, {\mathbf{x}}_n(-v)},\, \forall\, x(v)\;$$$$\mbox{such that $(x(v),{\mathbf{x}}_n(-v))\in E$}$$$$\bigr)\,.$$

  4. The values of all components $w\not= v$ are not changed, i.e. $x_{n+1}(w)=x_n(w)$ for all $w\not= v$.

Theorem 3.11   $\;$ Let the transition matrix ${\mathbf{P}}=(p_{{\mathbf{x}}{\mathbf{x}}^\prime})$ be given as

$$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)

where the conditional probabilities $\pi_{x^\prime(v)\mid\, {\mathbf{x}}(-v)}$ are defined in $(\ref{def.deb.gib})$. Then ${\mathbf{P}}$ is irreducible and aperiodic and the pair $({\mathbf{P}},{\boldsymbol{\pi}})$ is reversible.

Proof

$\;$ The assertion can be proved similarly to the proof of Theorem 3.10.

• In order to see that ${\mathbf{P}}=(p_{{\mathbf{x}}{\mathbf{x}}^\prime})$ is aperiodic it suffices to notice
  • that for all ${\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$$

  • and hence all diagonal elements $p_{{\mathbf{x}},{\mathbf{x}}}$ of ${\mathbf{P}}$ are positive.
• The following considerations show that ${\mathbf{P}}$ is irreducible.
  • For arbitrary but fixed ${\mathbf{x}},{\mathbf{x}}^\prime\in E$ let $k\le\vert V\vert$ be the number of components $v\in V$ such that $x(v)\not=x^\prime(v)$.
  • For $k=0$, i.e. ${\mathbf{x}}={\mathbf{x}}^\prime$, we already showed while proving the aperiodicity that $p_{{\mathbf{x}}{\mathbf{x}}}>0$.
  • Let now $k>0$. Without loss of generality we may assume that the components are linearly ordered and that the first $k$ components of ${\mathbf{x}}$ and ${\mathbf{x}}^\prime$ differ.
  • By hypothesis (34) the state space $E$ contains a sequence ${\mathbf{y}}_0,\ldots,{\mathbf{y}}_k\in E$ such that ${\mathbf{y}}_0={\mathbf{x}}$ and
    $${\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\,.$$

  • Moreover, for each $i=0,\ldots,k-1$

    $$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)

    
    
    and thus $p^{(k)}_{{\mathbf{x}}{\mathbf{x}}^\prime}\ge\prod_{i=0}^{k-1} p_{{\mathbf{y}}_i{\mathbf{y}}_{i+1}}>0$.
• It is left to show that the detailed balance equation (2.85) holds, i.e.

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

  
  
  • If ${\mathbf{x}}={\mathbf{x}}^\prime$, then (38) obviously holds.
  • If $x(v)\not=x^\prime(v)$ for more than one component $v\in V$, then $p_{{\mathbf{x}}{\mathbf{x}}^\prime}=p_{{\mathbf{x}}^\prime{\mathbf{x}}}=0$ and hence (38) holds.
  • Let now $x(v)\not=x^\prime(v)$ for exactly one $v\in V$ (and hence $x(w)=x^\prime(w)$ for all $w\not= v$). Then (37) implies
    $$\pi_{\mathbf{x}}\,p_{{\mathbf{x}}{\mathbf{x}}^\prime}=\pi_{\mathb... ...hbf{z}}}\; = \pi_{{\mathbf{x}}^\prime}\,p_{{\mathbf{x}}^\prime{\mathbf{x}}}\,.$$
    
     
      $\Box$
    
    
    

Let ${\mathbf{X}}_0,{\mathbf{X}}_1,\ldots$ be a Markov chain with state space $E$ and the transition matrix ${\mathbf{P}}=(p_{{\mathbf{x}}{\mathbf{x}}^\prime})$ given by (36). As a consequence of Theorem 3.11 we get that in this case

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

for any initial concentration ${\boldsymbol{\alpha}}_0$ where ${\boldsymbol{\alpha}}_n$ denotes the distribution of ${\mathbf{X}}_n$. Furthermore, the Gibbs sampler shows the following monotonic behavior.

Theorem 3.12   $\;$ For all $n=0,1,\ldots$,

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

Proof

 

• For arbitrary $v\in V$ and ${\mathbf{x}}^\prime\in E$, formula (35) implies
  

  $$\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)

  
  

• Using this and the definition (36) of the transition matrix ${\mathbf{P}}=(p_{{\mathbf{x}}{\mathbf{x}}^\prime})$ we obtain
  

  $$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$
  
  
  

Remarks

 

• A modified version of the Gibbs sampler that was considered in this section is the so-called cyclic Gibbs sampler, which uses a different procedure for picking the component $v\in V$ that will be updated.
  • Namely, it is not chosen according to a (given) probability function ${\mathbf{q}}=(q_v,\,v\in V)$, where $q_v>0$ for all $v\in V$,
  • but the components $v\in V$ are sorted linearly and chosen one after another according to this order. The selection of the update candidates thus becomes a deterministic procedure.
• If $k=n\vert V\vert+i$ for some numbers $n=0,1,\ldots$ and $i=1,\ldots,\vert V\vert$, then the matrix ${\mathbf{P}}(k)=\bigl(p_{{\mathbf{x}}{\mathbf{x}}^\prime}(k)\bigr)$ of the transition probabilities $p_{{\mathbf{x}}{\mathbf{x}}^\prime}(k)$ in step $k$ is given as

  $$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)

  
  

• For an entire (scan) cycle, updating each component exactly once, one obtains the following transition matrix

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

  
  

• It is easy to show that the matrix ${\mathbf{P}}=(p_{{\mathbf{x}}{\mathbf{x}}^\prime})$ given by (42) and (43)
  • is irreducible and aperiodic
  • and that ${\boldsymbol{\pi}}$ is the stationary (limit) distribution of ${\mathbf{P}}$ as
  • for all $i=1,\ldots,\vert V\vert$ and for all ${\mathbf{x}}^\prime\in E$ formulae (41) and (42) imply that
    $$\sum\limits_{{\mathbf{x}}\in E} \pi_... ...\prime(-v_i)\bigr) \stackrel{(\ref{hil.equ.for})}{=} \pi_{{\mathbf{x}}^\prime}$$

  • and hence also
    $$\sum\limits_{{\mathbf{x}}\in E} \pi_{\mathbf{x}}p_{{\mathbf{x}}{\mathbf{x}}^\prime}= \pi_{{\mathbf{x}}^\prime}\,.$$

• The pair $({\mathbf{P}},{\boldsymbol{\pi}})$ is in general not reversible. However, in Section 2.3.4 we showed that the pair $({\mathbf{M}},{\boldsymbol{\pi}})$ is reversible where

  $${\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)

  
  
  denotes the multiplicative reversible version of ${\mathbf{P}}$.

Theorem 3.13   $\;$ The matrix ${\mathbf{M}}$ has the following representation

$${\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)

i.e., the multiplicative reversible version ${\mathbf{M}}$ of the ,,forward-scan matrix'' ${\mathbf{P}}$ coincides with the ,,forward-backward scan matrix''.

Proof

 

• It suffices to show that $\widetilde{\mathbf{P}}={\mathbf{P}}(\vert V\vert)\cdot\ldots\cdot{\mathbf{P}}(1)$ for the matrix $\widetilde{\mathbf{P}}=(\widetilde p_{{\mathbf{x}}{\mathbf{x}}^\prime})$ defined by (44).
• Formulae (42)-(44) imply for arbitrary ${\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}\,.$$

  
  

• This and (35) yield (similar to the proof of (38))
  

  $${\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$
  
  
  

Remarks

 

• If Gibbs samplers are used in practice it is always assumed
  • that the conditional probabilities considered in (36) and (42)
    $$\pi_{x(v)\mid\, {\mathbf{x}}(-v)}=P\bigl(X(v)=x(v)\mid\, {\mathbf{X}}(-v)={\mathbf{x}}(-v)\bigr)$$
    only depend on the vector $\bigl(x(w),\,w\in\mathcal{N}(v)\bigr)$ of the values
  • obtained by the random vector ${\mathbf{X}}=(X(w),\,w\in V)$ in a certain small neighborhood $\mathcal{N}(v)\subset V$ of $v\in V$.
• The family $\mathcal{N}=\{\mathcal{N}(v),\,v\in V\}$ of subsets of $V$ is called a system of neighborhoods if for arbitrary $v,w\in V$

  (a)

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

  (b)

  $w\in\mathcal{N}(v)$ implies $v\in\mathcal{N}(w)$.

• For the hard-core model from Section 3.3.1, $\mathcal{N}(v)$ is the set of those vertices $w\not= v$ that are directly connected to $v$ by an edge.