Overrelaxation for Gibbs Sampler



In [2]:

    
%matplotlib inline
import numpy as np
import numpy.random as nr
import matplotlib.pyplot as plt



In [3]:

    
Mu = np.array([0,0])
Sigma = np.array([[1,0.998],[0.998,1]])
Lambda = np.linalg.inv(Sigma)

1. Gibbs Sampler without overrelaxation



In [7]:

    
N = 40 # number of samples
x0_samples = [0.5]
x1_samples = [0.5]
x1 = 0.5
for i in range(N):
    cond_mu0 = Mu[0] - Lambda[1,0]*(x1 - Mu[1])/Lambda[0,0]
    cond_sigma0 = 1/Lambda[0,0]
    x0 = nr.normal(cond_mu0, cond_sigma0)
    x0_samples.append(x0)
    
    cond_mu1 = Mu[1] - Lambda[0,1]*(x0 - Mu[0])/Lambda[1,1]
    cond_sigma1 = 1/Lambda[1,1]
    x1 = nr.normal(cond_mu1, cond_sigma1)
    x1_samples.append(x1)



In [32]:

    
plt.plot(x0_samples, x1_samples, 'r')
lim = 0.6
plt.xlim((-lim,lim))
plt.ylim((-lim,lim))









    Out[32]:





(-0.6, 0.6)

2. Gibbs Sampler With Overrelaxation



In [35]:

    
N = 40 # number of samples
x0_samples_or = [0.5]
x1_samples_or = [0.5]
alpha = -0.89
x1 = 0.5
for i in range(N):  
    cond_mu0 = Mu[0] - Lambda[1,0]*(x1 - Mu[1])/Lambda[0,0]
    cond_sigma0 = 1/Lambda[0,0]
    x0 = cond_mu0 + alpha * (x0_samples_or[-1] - cond_mu0) + np.sqrt(1-alpha**2)*nr.normal(0, cond_sigma0)
    x0_samples_or.append(x0)
    
    cond_mu1 = Mu[1] - Lambda[0,1]*(x0 - Mu[0])/Lambda[1,1]
    cond_sigma1 = 1/Lambda[1,1]
    x1 = cond_mu1 + alpha * (x1_samples_or[-1] - cond_mu1) + np.sqrt(1-alpha**2)*nr.normal(0, cond_sigma1)
    x1_samples_or.append(x1)



In [36]:

    
plt.plot(x0_samples_or, x1_samples_or)
lim = 0.6
plt.xlim((-lim,lim))
plt.ylim((-lim,lim))









    Out[36]:





(-0.6, 0.6)

3. Plot



In [37]:

    
plt.plot(x0_samples, x1_samples,'r')
plt.plot(x0_samples_or, x1_samples_or,'b')
lim = 0.6
plt.xlim((-lim,lim))
plt.ylim((-lim,lim))









    Out[37]:





(-0.6, 0.6)



In [ ]: