Introduction

It is not always easy to tune the proposal distribution. The Gibbs sampling is a procedure for multivariate distributioin in which all samples are accepted therefore we do not have to specify a proposal distribution, leaving some guessing work out of MCMC procedure. The Gibbs sampling can only be applied in situations where we know the full conditional distribution of each component in the multivariate distribution conditioned on all other components.

Sampling from conditionals

Given your posterior distribution has the following functional form: \begin{equation} f(x,y)= \frac{1}{C}e^{-\frac{(x^2y^2+x^2+y^2-8x-8y)}{2}} \end{equation}

After certain simplification, we find the conditional distribution have a relatively simple form. \begin{equation} p(x|y)= g(y)e^{-(x-\frac{4}{(1+y^2)})^2\frac{(1+y^2)}{2}} \end{equation} \begin{equation} p(y|x)= g(x)e^{-(y-\frac{4}{(1+x^2)})^2\frac{(1+x^2)}{2}} \end{equation}

Gibbs sampling in summary:

  • Set i = 1
  • Initialize $x^{(i)}=\textbf{u}$, where $\textbf{u}=(u_1,u_2)$
  • Repeat
    • Sample $x_1^{(i+1)}$ from the conditional $p(x_1|x_2=x_2^{(i)})$
    • Sample $x_2^{(i+1)}$ from the conditional $p(x_2|x_1=x_1^{(i+1)})$

Code Implementation

In [2]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import scipy as sp
import pandas as pd
pd.set_option('display.width',500)
pd.set_option('display.max_columns',100)

import seaborn as sns
In [3]:
f = lambda x,y: np.exp(-(x*x*y*y+x*x+y*y-8*x-8*y)/2)
In [23]:
xx = np.linspace(-1,8,100)
yy = np.linspace(-1,8,100)

xg, yg = np.meshgrid(xx,yy)

z = f(xg.ravel(), yg.ravel())
z2 = z.reshape(xg.shape)

plt.contourf(xx,yy,z2,cmap='BuGn')
plt.colorbar()
plt.xlabel('x')
plt.ylabel('y')
plt.show()

Gibbs sampler

In [40]:
unif = np.random.uniform
rnorm = np.random.normal

N = 50000

x = np.zeros(N+1)
y = np.zeros(N+1)

x[0] = unif()
y[0] = unif()

sig = lambda z,i: np.sqrt(1/(1+z[i]**2))
mu = lambda z,i: 4/(1+z[i]**2)

for i in np.arange(1,N+1):
    x[i] = rnorm(mu(y,i-1),sig(y,i-1))
    y[i] = rnorm(mu(x,i),sig(x,i))
In [55]:
# show histgram of x
plt.figure(figsize=(10,4))
plt.subplot(1,2,1)
plt.hist(x,bins=50)
plt.subplot(1,2,2)
plt.hist(y,bins=50)

plt.show()
In [72]:
# fitting effect through scatter plot
plt.contourf(xg,yg,z2,cmap='BuGn',alpha = 0.8)
plt.plot(x,y,'r.',alpha =0.1,markersize=1)

plt.show()
In [ ]: