Gibbs sampling is a Markov chain Monte Carlo method for approximating high-dimensional posterior distributions. It updates one variable at a time, sampling from that variable’s conditional distribution given the latest values of all other variables. section 1.4.1 section 2.9
A Gibbs iteration proceeds as follows:
- Choose a starting state for the multidimensional variable. section 2.5
- Cycle through the dimensions, conditioning each update on the latest values of every other dimension. section 2.5
- Draw from each fully conditional distribution and record the resulting state. section 2.5
- Repeat the cycle to obtain samples that approximate the target distribution. section 2.9
The method is practical when the fully conditional distributions are easy to compute and sample. Naive Gibbs sampling is sequential, requiring a number of update steps proportional to the distribution’s dimensionality. section 2.5
Sampling a multivariate Gaussian
For a multivariate Gaussian, partition the variables into sets, compute the relevant conditional Gaussian distribution, and repeatedly sample each dimension while conditioning on the others. A from-scratch Python implementation separates conditional-distribution calculation into gauss_conditional and iteration into sample, using NumPy, SciPy, and a visualized sample cloud. section 2.5
Figure 2.11 Gibbs sampler pseudo-code
Use the Gibbs sampler pseudocode to inspect the update sequence and repeated conditional draws. section 2.5
Listing 2.4 Gibbs sampling
import numpy as np
import matplotlib.pyplot as plt
import itertools
from numpy.linalg import inv
from scipy.stats import multivariate_normal
np.random.seed(42)
class gibbs_gauss:
def gauss_conditional(self, mu, Sigma, setA, x): #1
dim = len(mu)
setU = set(range(dim))
setB = setU.difference(setA)
muA = np.array([mu[item] for item in setA]).reshape(-1,1)
muB = np.array([mu[item] for item in setB]).reshape(-1,1)
xB = np.array([x[item] for item in setB]).reshape(-1,1)
Sigma_AA = []
for (idx1, idx2) in itertools.product(setA, setA):
Sigma_AA.append(Sigma[idx1][idx2])
Sigma_AA = np.array(Sigma_AA).reshape(len(setA),len(setA))
Sigma_AB = []
for (idx1, idx2) in itertools.product(setA, setB):
Sigma_AB.append(Sigma[idx1][idx2])
Sigma_AB = np.array(Sigma_AB).reshape(len(setA),len(setB))
Sigma_BB = []
for (idx1, idx2) in itertools.product(setB, setB):
Sigma_BB.append(Sigma[idx1][idx2])
Sigma_BB = np.array(Sigma_BB).reshape(len(setB),len(setB))
Sigma_BB_inv = inv(Sigma_BB)
mu_AgivenB = muA + np.matmul(np.matmul(Sigma_AB, Sigma_BB_inv),
xB - muB)
Sigma_AgivenB = Sigma_AA - np.matmul(np.matmul(Sigma_AB,
Sigma_BB_inv), np.transpose(Sigma_AB))
return mu_AgivenB, Sigma_AgivenB
def sample(self, mu, Sigma, xinit, num_samples):
dim = len(mu)
samples = np.zeros((num_samples, dim))
x = xinit
for s in range(num_samples):
for d in range(dim):
mu_AgivenB, Sigma_AgivenB = self.gauss_conditional(mu, Sigma,
set([d]), x)
x[d] = np.random.normal(mu_AgivenB, np.sqrt(Sigma_AgivenB))
#end for
samples[s,:] = np.transpose(x)
#end for
return samples
if __name__ == "__main__":
num_samples = 2000
mu = [1, 1]
Sigma = [[2,1], [1,1]]
xinit = np.random.rand(len(mu),1)
num_burnin = 1000
gg = gibbs_gauss()
gibbs_samples = gg.sample(mu, Sigma, xinit, num_samples)
scipy_samples =
multivariate_normal.rvs(mean=mu,cov=Sigma,size=num_samples,random_state=42)
plt.figure()
plt.scatter(gibbs_samples[num_burnin:,0], gibbs_samples[num_burnin:,1],
label='Gibbs Samples')
plt.grid(True); plt.legend(); plt.xlim([-4,5])
plt.title("Gibbs Sampling of Multivariate Gaussian"); plt.xlabel("X1");
plt.ylabel("X2")
plt.show()
plt.figure()
plt.scatter(scipy_samples[num_burnin:,0], scipy_samples[num_burnin:,1],
label='Ground Truth Samples')
plt.grid(True); plt.legend(); plt.xlim([-4,5])
plt.title("Ground Truth Samples of Multivariate Gaussian");
plt.xlabel("X1");
plt.ylabel("X2")
plt.show()
Use the implementation listing to trace conditional Gaussian construction, dimension cycling, and sample generation. section 2.5
