#Computer practical on Gibbs sampling
#Sampling from a bivariate normal distribution.
	
#Implement the Gibbs sampling algorithm to sample from the bivariate normal distribution with mean c(0,0) and covariance matrix matrix(c(4,1,1,4), nrow=2)

#Tasks 1 and 3
#set the mean parameter
mu <- c(0,0)
#set the covarince parameter
sigma <- matrix(c(4,1,1,4), nrow=2)

#Gibbs sampling
#set the desired sample size
n <- 3000
x <- matrix(nrow=n, ncol=2)

#set the starting value
cur.x <- c(0,0)

#run the Gibbs sampler; estimate P(X_1 >=0, X_2 >=0)
prop <- rep(0, times=n)

sd1 <- sqrt(sigma[1,1] - sigma[1,2]^2/sigma[2,2])
sd2 <- sqrt(sigma[2,2] - sigma[1,2]^2/sigma[1,1])  

for(i in 1:n){
  mean1 <- mu[1] + sigma[1,2]/sigma[2,2] * (cur.x[2] - mu[2])
  cur.x[1] <- rnorm(1, mean = mean1, sd = sd1)
  
  mean2 <- mu[2] + sigma[1,2]/sigma[1,1] * (cur.x[1] - mu[1])
  cur.x[2] <- rnorm(1, mean = mean2, sd = sd2)
  
  if(i == 1){
  	prop[i] <- as.numeric((cur.x[1] >= 0) && (cur.x[2] >= 0))
  	}
  	else{
  		prop[i] <- as.numeric((cur.x[1] >= 0) && (cur.x[2] >= 0)) + prop[i-1]
  		}

  x[i,] <- cur.x
}

prop <- prop/1:n

#Task 2
#plot the last 100 values, with subsequent values linked by a line
plot(x[901:1000,], type='b', xlab="X_1", ylab="X_2")

#look at sample plots of both variables
par(mfrow=c(2,1))
plot(x[,1], type="l", xlab="t", ylab="X_1")
plot(x[,2], type="l", xlab="t", ylab="X_2")

#look at the estimate of P(X_1 >=0, X_2 >=0)
plot(prop, ylim=c(0.1,0.4), type='l', xlab="t", ylab="Estimate of P(X_1>=0, X_2>=0)")

#Task 4
#Now set the covarince parameter as follows, and repeat.
sigma <- matrix(c(4,2.8,2.8,2), nrow=2)