Gibbs Sampler

Example 1 (Glauber Dynamics for Ising Model)

Related: Ising Model

States:

$$S=\{\underline{\sigma }=({\sigma }_1,\dots ,{\sigma }_m)\mid {\sigma }_j=+1,-1\}$$

Target:

$$\pi ,\pi (\underline{\sigma })=\frac{b(\underline{\sigma })}{Z}>0\forall \underline{\sigma }\in S$$

Notation:

$$\begin{array}{rl}{\underline{\sigma }}^{j,+1}& =({\sigma }_1,\dots ,{\sigma }_{j-1},+1,{\sigma }_{j+1},\dots ,{\sigma }_m)\\ {\underline{\sigma }}^{j,-1}& =({\sigma }_1,\dots ,{\sigma }_{j-1},-1,{\sigma }_{j+1},\dots ,{\sigma }_m)\end{array}$$

Here we have a high dimensional state space. The goal is to sample from a distribution $\pi$ proportional to the weight function $b$. The constant $Z$, of course, is the partition function, which is often impossible to compute. We use the following notation to represent updates to the entire system.

Glauber Dynamics (Special Case of Gibbs Sampler)

The idea is to use conditional distributions to determine the probability of spin ${\sigma }_j$ or $V_j$ taking value $\tau$ where $\tau \in \{+1,-1\}$. For simplicity, let $\underline{\sigma }=\underline{V}$ and $\underline{V}\sim \pi$, so $\underline{V}=(V_1,\dots ,V_m)\in S$. Consider the condition distribution of $V_j$ given by $V_i$ for $i\ne j$:

$$\begin{array}{rl}\mathbb{P}(V_j=\tau \mid V_1,\dots V_{j-1},V_{j+1},\dots ,V_m)& =\frac{\pi (V_1,\dots ,V_{j-1},\tau ,V_{j+1},\dots ,V_m)}{\pi (V_1,\dots ,+1,\dots ,V_m)+\pi (V_1,\dots ,-1,\dots ,V_m)}\end{array}$$

where the $+1,-1$ in the denominator denotes the $j$th entry. This is only nonzero for $\tau =+1,-1$. Let

$$\begin{array}{rl}{\pi }_j(\tau \mid {\sigma }_1,\dots ,{\sigma }_{m-1})& =\mathbb{P}(V_j=\tau \mid (V_1,\dots ,V_{j-1},V_{j+1},\dots ,V_m)=({\sigma }_1,\dots ,{\sigma }_{m-1}))\end{array}$$

for $\underline{V}\sim \pi$. Caution: Some conventions use $\pi (\cdot \mid \cdot )$ instead of ${\pi }_j(\cdot \mid \cdot )$.

Glauber Dynamics (Random-scan Gibbs sampler)

Given $X_k={\underline{x}}_k=(x_{k,1},\dots ,x_{k,m})\in S$.

1. Choose $J$ uniformly at random from $\{1,2,\dots ,m\}$.
2. Sample $Y$ from the conditional distribution $${\pi }_J(\cdot \mid x_{k,1},\dots ,x_{k,J-1},x_{k,J+1},\dots ,x_{k,m})$$
3. Set $X_{k+1}\leftarrow (x_{k,1},\dots ,x_{k,J-1},Y,x_{k,J+1},\dots ,x_{k,m})$.

Algorithmic Implementation (Discrete State Space)

For a discrete system like the Ising Model where ${\sigma }_j\in \{-1,+1\}$, Step 2 is implemented by calculating the explicit scalar probabilities for the two possible states of the $J$-th coordinate.

Assume $X_k=\underline{\sigma }\in S$.

1. Set ${\alpha }_{+}$ as the probability of the $J$-th spin being $+1$:

$${\alpha }_{+}={\pi }_J(+1\mid {\sigma }_{-J})=\frac{\pi ({\underline{\sigma }}^{J,+1})}{\pi ({\underline{\sigma }}^{J,-1})+\pi ({\underline{\sigma }}^{J,+1})}$$

2. Set ${\alpha }_{-}$ as the probability of the $J$-th spin being $-1$:

$${\alpha }_{-}={\pi }_J(-1\mid {\sigma }_{-J})=\frac{\pi ({\underline{\sigma }}^{J,-1})}{\pi ({\underline{\sigma }}^{J,-1})+\pi ({\underline{\sigma }}^{J,+1})}$$

Note ${\alpha }_{+}+{\alpha }_{-}=1$. The partition function $Z$ cancels out in these fractions.

3. Branch execution:

• Set $X_{k+1}\leftarrow {\underline{\sigma }}^{J,+1}$ with probability ${\alpha }_{+}$.
• Set $X_{k+1}\leftarrow {\underline{\sigma }}^{J,-1}$ with probability ${\alpha }_{-}$.

Here, this requires a discrete state space, where ${\sigma }_j\in \{-1,+1\}$.

Example 2

Each “spin” is in $\{-1,+1\}$. Set $J$th spin conditional on the value of every other spin. For example, let $m=2$, such that

$$\pi (-1,-1)=\frac{1}{Z},\pi (+1,-1)=\frac{2}{Z},\pi (-1,+1)=\frac{3}{Z},\pi (+1,+1)=\frac{5}{Z}$$

Suppose $X_k=(-1,-1)$ and $J=2$. Then

$$\begin{array}{rl}{\underline{\sigma }}^{J,+1}& =(-1,+1)\\ {\underline{\sigma }}^{J,-1}& =(-1,-1)\end{array}$$

such that

$$\begin{array}{rl}{\alpha }_{+}& =\frac{3}{1+3}=\frac{3}{4}\\ {\alpha }_{-}& =\frac{1}{1+3}=\frac{1}{4}\end{array}$$

We can see that ${\alpha }_{+}+{\alpha }_{-}=1$. Then,

$$\mathbb{P}(X_{k+1}=(-1,+1)\mid X_k=(-1,-1))=\frac{1}{2}\cdot \frac{3}{4}=\frac{3}{8}$$

where the first product term is the probability of selecting particle $2$, or $1/m$, and the second is ${\alpha }_{+}$.

Example 3

Assume we have the same setup as Example 2. Then the conditional distribution ${\pi }_2$ is

$$\begin{array}{rl}{\pi }_2(-1\mid +1)& =\mathbb{P}(\text{entry 2}=-1\mid \text{entry 1}=+1)\\ & =\frac{\pi (+1,-1)}{\pi (+1,-1)+\pi (+1,+1)}\\ & =\frac{2}{7}\end{array}$$

Likewise,

$$\begin{array}{rl}{\pi }_2(+1\mid +1)& =\frac{\pi (+1,+1)}{\pi (+1,-1)+\pi (+1,+1)}\\ & =\frac{5}{7}\end{array}$$

Remarks:

• ${\pi }_2(\cdot \mid +1)$ is a PMF on $\{-1,+1\}$, just like ${\pi }_1(\cdot \mid -1),{\pi }_1(\cdot \mid +1),{\pi }_2(\cdot \mid -1)$
• Conditional distributions do not require knowledge on the normalization $Z$.

General Case of the Gibbs Sampler

The goal of this algorithm is to sample from complex, high-dimensional probability distributions that are difficult to compute directly.

Ising model	$⇝$	General case
$\underline{\sigma }=({\sigma }_1,\dots ,{\sigma }_m)$	$⇝$	$\underline{x}=(x_1,\dots ,x_m)$
${\sigma }_j\in \{-1,+1\}$	$⇝$	$x_j\in \mathbb{R}$
$\pi (\underline{\sigma })=\pi ({\sigma }_1,\dots ,{\sigma }_m)$	$⇝$	$f(\underline{x})=f(x_1,\dots ,x_m)$ (unknown normalization)
$\pi ({\sigma }_j\mid {\sigma }_1,\dots ,{\sigma }_{j-1},{\sigma }_{j+1},\dots ,{\sigma }_m)$	$⇝$	$f(x_j\mid x_1,\dots ,x_{j-1},x_{j+1},\dots ,x_m)$
PMF on ${\sigma }_j\in \{-1,+1\}$	$⇝$	PDF on $x_j\in \mathbb{R}$

The conditional distribution does not require the global normalization constant. The continuous conditional distribution is computed as:

$$f_j(y\mid x_1,\dots ,x_{m-1})=\frac{f(x_1,\dots ,x_{j-1},y,x_{j+1},\dots ,x_m)}{{\int }_{\mathbb{R}}f(x_1,\dots ,x_{j-1},y,x_{j+1},\dots ,x_m)dy}$$

So $f_j(\cdot \mid x_1,\dots ,x_{m-1})$ is a PDF on $\mathbb{R}$. Note that the denominator is an integral, in the same way that in the discrete space, we sum the possibilities for that flip. This denotes the space of all possibilities that $x_j$ can become. In our case, we care about when $\mathbb{P}(x_j=y)$ specifically.

Example 4

Let

$$f(x,y)=c\exp \left(-(x^2+1)(y^2+1)\right),(x,y)\in {\mathbb{R}}^2$$

$c$ is a normalization, so

$${\int }_{{\mathbb{R}}^2}c\exp \left(-(x^2+1)(y^2+1)\right)dxdy=1$$

What are the conditional distributions? We see that

$$\begin{array}{rl}{f}_1(x\mid y)& =\frac{c\exp \left(-(x^2+1)(y^2+1)\right)}{{\int }_{\mathbb{R}}f(x,y)dx}\\ & =c_1(y)\exp \left(-(y^2+1)x^2\right)\end{array}$$

by treating $y$ as a constant and expanding out $f(x,y)$. Here, we apply the continuous conditional distribution from General Case of the Gibbs Sampler. Then, $c_1(y)$ is such that

$${\int }_{-\infty }^{\infty }{c}_1(y)\exp \left(-(y^2+1)x^2\right)dx=1$$

for all $y\in \mathbb{R}$. In other words, $c_1(y)$ is the normalization of the PDF

$$f_1(x\mid y)\propto \exp \left(-(y^2+1)x^2\right)$$

Recall that $Z\sim N(0,{\sigma }^2)$ has a PDF of

$$f_Z(z)\propto \exp \left(\frac{-1}{2{\sigma }^2}{z}^2\right)$$

This implies that $f_1(x\mid y)$ is the PDF of $N\left(0,\frac{1}{2(y^2+1)}\right)$. Similarly, $f_2(y\mid x)$ is the PDF of $N\left(0,\frac{1}{2(x^2+1)}\right)$.

Algorithm (Random-Scan Gibbs Via Normal Distribution)

The previous examples gives the following insight: We can model the probabilities as part of the Normal distribution.

Given $(X_k,Y_k)\in {\mathbb{R}}^2$.

1. Choose $J$ uniformly at random from $\{1,2\}$. Note that we are in ${\mathbb{R}}^2$.
2. Generate $Z\sim N(0,1)$. Recall Box-Muller Method.
3. If $J=1$, set $$W=\frac{1}{\sqrt{2(Y_k^2+1)}}Z(\text{i.e. }W\sim f_1(\cdot \mid Y_k))$$ and set $(X_{k+1},Y_{k+1})\leftarrow (W,Y_k)$.
4. If $J=2$, set $$W=\frac{1}{\sqrt{2(X_k^2+1)}}Z(\text{i.e. }W\sim f_2(\cdot \mid X_k))$$ and set $(X_{k+1},Y_{k+1})\leftarrow (X_k,W)$.

Algorithm (Systematic-Scan Gibbs Sampler)

Given ${\underline{X}}_k\in {\mathbb{R}}^m$:

1. Sample $Y_1$ from $f_1(\cdot \mid X_{k,2},\dots ,X_{k,m})$
2. For $i=2$ to $m$: Sample $Y_i$ from $f_i(\cdot \mid \underset{\in {\mathbb{R}}^{m-1}}{\underbrace{Y_1,\dots ,Y_{i-1},X_{k,i+1},\dots ,X_{k,m}}})$
3. Set ${\underline{X}}_{k+1}=(Y_1,\dots ,Y_m)\in {\mathbb{R}}^m$

Output: ${\underline{X}}_0,{\underline{X}}_1,\cdots \in {\mathbb{R}}^m$

Remarks:

• The Systematic-Scan chooses coordinates $j=1,2,\dots ,m$ in strict order.
• The Random-Scan Gibbs Sampler is a special case of the Metropolis algorithm.

Metropolis Algorithm (Random-Walk Sampler)

We can extend the Metropolis Algorithm.

• Goal: Sample from target PDF $f(\underline{x})$ where $\underline{x}\in {\mathbb{R}}^m$.
• Given ${\underline{X}}_k\in {\mathbb{R}}^m$,

1. Generate $\underline{Z}=(Z_1,\dots ,Z_m)\in {\mathbb{R}}^m$ where $Z_1,\dots ,Z_m\sim N(0,1)$ iid.

   I.e., $\underline{Z}\sim N(\underline{0},\Sigma )$ is Multivariate Normal. This is also our “noise”, or a random directional perturbation.

2. Set proposal state $\underline{Y}={\underline{X}}_k+\underline{Z}$.

   The proposal transition matrix $q({\underline{X}}_k,\cdot )$ is a PDF of $N({\underline{X}}_k,I_m)$

3. Set $$\alpha =\min \left(1,\frac{f(\underline{Y})}{f({\underline{X}}_k)}\right)$$ and generate Uniform $U\sim \mathcal{U}[0,1]$.
   • Set ${\underline{X}}_{k+1}\leftarrow \underline{Y}$ if $U\le \alpha$.
   • Set ${\underline{X}}_{k+1}\leftarrow {\underline{X}}_k$ if $U<\alpha$.

Remarks:

• The output ${\underline{X}}_0,{\underline{X}}_1,\dots$ is a Markov Chain on a continuous state space ${\mathbb{R}}^m$.
• It is “sensitive” to the typical “jump” size $\underline{Z}$.
• What if happens if we take jumps $\ell Z$ and let $\ell \to 0$? We get SDEs!
• $\ell$ is the spatial exploration rate of th random walk.

Remarks

1. We first define the Ising Model and use Glauber Dynamics algorithm to figure out how to sample from these difficult distributions. This is the most restricted class of MCMC: a discrete-space, single-coordinate update mechanism.
2. We generalize this via the Metropolis Algorithm to generalize transitions by permitting global state proposals via a transition matrix $Q$. At first, we let the transitions be symmetric, to show detailed balance.
3. We further generalized this to show asymmetry via the Metropolis Hastings Algorithm by having a ratio to penalize directional bias in the proposal distribution, forcing detailed balance.
4. We expand the problem domain from discrete states $\{-1,+1\}$ to the infinite continuous domain of $\mathbb{R}$. The target PMF $\pi (\underline{\sigma })$ becomes PDF $f(\underline{x})$. We generalize the discrete scalar fractions with continuous integrals to compute the conditional $f_j(y\mid x_{-j})$. The core mechanism is to update one dimension at a time while holding all others constant to bypass the global normalization integral.
5. We can generalize to multidimensional continuous transitions using additive noise in Metropolis Algorithm (Random-Walk Sampler). The proposal mechanism becomes a continuous Gaussian transition kernel.
6. Finally, we can further generalize this to continuous time steps by introducing the scaling factor $\ell$ to modify the variance of the Gaussian proposal jump to $\ell \underline{Z}$. By taking the limit $\ell \to 0$ and the discrete time step interval $\Delta t\to 0$ proportionally, the discrete Markov Chain converges into a continuous-time mathematical model. This leads us to Stochastic Differential Equations (SDEs) over $t\in {\mathbb{R}}^{+}$.