##################################################
#
# Sampling from a multivariate normal via
# univariate draws from the full conditional
# normal distributions.
#
##################################################

p   = 6
mm  = 0
vv  = 1
rho = 0.99
mu  = rep(mm,p)
Sig = diag(vv,p)
for (i in 1:p)
  for (j in 1:p)
    if (i!=j){
      Sig[i,j]=rho
    }

# Sampling jointly from p(x1,...xp)    
library(MASS)
set.seed(12345)
M  = 10000
xx = mvrnorm(M,mu,Sig)

par(mfrow=c(2,p))
for (i in 1:p)
  ts.plot(xx[,i],ylab="",xlab="Iteration",main=paste("Component ",i))
for (i in 1:p)
  acf(xx[,i],main="")


# Sampling from p(x1,...,xp) via univariate Gibbs steps
set.seed(12345)
M0= 1000
M = 5000
niter = M0+M
x = rnorm(p)
xs = matrix(0,niter,p)
xs[1,]=x
for (iter in 2:niter){
for (i in 1:p){
  mean = mu[i] + Sig[i,-i]%*%solve(Sig[-i,-i])%*%(x[-i]-mu[-i])
  var  = Sig[i,i] - Sig[i,-i]%*%solve(Sig[-i,-i])%*%Sig[-i,i] 
  x[i] = rnorm(1,mean,sqrt(var))
}
xs[iter,] = x
}

par(mfrow=c(2,p))
for (i in 1:p)
  ts.plot(xs[,i],ylab="",xlab="Iteration",main=paste("Component ",i))
for (i in 1:p)
  acf(xs[(M0+1):M,i],main="")

par(mfrow=c(2,3))
for (i in 1:p){
  plot(density(xx[,i]),main=paste("Component ",i),xlab="")
  lines(density(xs[(M0+1):M,i]),col=2)
}