# Adapted from:
# http://www.mas.ncl.ac.uk/~ndjw1/teaching/sim/gibbs/gibbs.html
# and Computational Statistics with MATLAB

# Bivariate normal example

# Direct generation of samples
rbvn<-function (n,mu,sigma,rho) 
{
  z1 <- rnorm(n) # get a standard normal sample
  z2 <- rnorm(n)
  x <- mu[1] + sigma[1]*z1
  y <- mu[2] + sigma[2]*(rho*z1 + sqrt(1-rho^2)*z2)
  cbind(x, y)
}
mu <- c(1,4) # means of normal RVs
sigma <- c(3,2) # standard deviations
rho <- 0.5
bvn<-rbvn(10000,mu,sigma,rho)
plot(bvn[,1],bvn[,2],pch=".")

# Samples means and correlation
mean(bvn[,1]) 
mean(bvn[,2])
cor(bvn[,1],bvn[,2])

# Generate samples using Gibbs sampler
n <- 100000
xgibbs <- matrix(0,n,2)
# Initial point
xgibbs[1,] <- c(10,10)
# Run the chain
for (i in 2:n) {
  mean1 <- mu[1] + rho*sigma[1]/sigma[2]*(xgibbs[i-1,2]-mu[2])
  xgibbs[i,1] <- mean1 + sigma[1]*sqrt(1-rho^2)*rnorm(1) # only needs standard normal univariate sampling
  mean2 <- mu[2] + rho*sigma[2]/sigma[1]*(xgibbs[i,1] - mu[1])
  xgibbs[i,2] <- mean2 + sigma[2]*sqrt(1-rho^2)*rnorm(1)
}
xgibbs_adj <- xgibbs[seq(from=5000,to=n,by=10),] # Discard initial portion and thin
plot(xgibbs_adj[,1],xgibbs_adj[,2],pch=".")
mean(xgibbs_adj[,1]) 
mean(xgibbs_adj[,2])
cor(xgibbs_adj[,1],xgibbs_adj[,2])