#ST407 Monte Carlo methods
#Practical 5 2012
#Using Gibbs samplers to sample from the Bivariate Normal distribution.
#The main difference between Gibbs and Metropolis Hastings as samplers is that we use the full conditional distributions as proposal distributions.
#Such proposals mean that we always accept the proposal when sampling from the proposal distribution.
#The full conditionals in this in this question are based on the conditional distribiutions of the Bivariate normal.
#the "full" part comes into the fact that WE CONDITIONAL ON THE MOST RECENT VALUES available
#Part 1
gibbs<-function(N,p){
#Gibbs for bivariate Normal
#N = number of iterations
#p = rho in the question, the covariance between X^(1) and X^(2)
x<-matrix(NA,N,2) #Storage to store our samples from the target distribution
x[1,]<-c(0,0) #Initialisation
for(i in 2:N){
x[i,1] <-rnorm(1, mean=p*x[i-1,2], sd=sqrt(1-p^2))#Full conditional for x^(1)
#We use the most recent value of x^(2) which is from the previous iteration i-1
x[i,2] <-rnorm(1, mean=p*x[i,1], sd=sqrt(1-p^2))#Full conditional for x^(2)
#We use the most recent value of x^(1) which is from the CURRENT iteration i
}
return(x)
}
Gibbs1.p0.01<-gibbs(N=1000, p=0.01)
plot(Gibbs1.p0.01, type="s") #Explores the statespace well.
#Comparable plot if we
acf(Gibbs1.p0.01) #No significant AC within components and between components
ts.plot(ts(Gibbs1.p0.01))
Gibbs2.p0.99<-gibbs(N=1000, p=0.99)
plot(Gibbs2.p0.99, type="s") #Very strong correlation, reflecting the large correlation in the covariance matrix
acf(Gibbs2.0.99) #Strong AC within components and between components.
#This is not desired if we want to produce independent samples from the Bivariate Normal.
#Thinning is possible to achieve independent samples, but this is computationally wasteful
ts.plot(ts(Gibbs2.p0.99)) #Time series indicates the high autocorrelation plot of the samples.
#Part 2
#Gibbs sampler Y^(1) Y^(2), a reparameterised version on X^(1) X^(2)
#The mean, variance and covariance of Y's can be computed by standard manipulation of expectations, variances and covariances
gibbs1<-function(N,p){
#Gibbs for bivariate Normal eigenvectors
#N=number of iterations
#p=rho, the correlation between the two components
#x in the Elke's code, is y in the notes and in my code
y<-matrix(NA,N,2) #Storage for our samples
y[1,]<-c(0,0) #Initialisation
for(i in 2:N){
y[i,1]<-rnorm(1, mean=0, sd=sqrt(1+p)) #Full conditional for component 1
y[i,2]<-rnorm(1, mean=0, sd=sqrt(1-p)) #Full conditional for component 2
#Note that these full conditionals are independent of the other component they are conditioned on+
}
return(y)
}
gibbs1.part2<-gibbs1(N=1000, p=0.99)
plot(gibbs1.part2, type="s")
acf(gibbs1.part2) #The AC is not significant both within and amongst components
ts.plot(gibbs1.part2) #Good set of independent samples
#Part 3
#The distribution of X^(1), X^(2) turns out to be initial joint distribution that we were interested in part 1
#So the idea is to to sample Y^(1), Y^(2) via Gibbs and then transform these samples to obtain X^(1), X^(2) samples
part3.xsamples<-gibbs1.part2 #Just some storage
part3.xsamples[,1]<-(gibbs1.part2[,1]+gibbs1.part2[,2])/sqrt(2)
part3.xsamples[,1]<-(gibbs1.part2[,1]-gibbs1.part2[,2])/sqrt(2)
plot(part3.xsamples, type="s")
acf(part3.xsamples) #The AC of the samples is now not significant between and amongst components
ts.plot(part3.xsamples) #And independent samples appear to be produced.
#So in general, independent samples are more frequent in the chain, and thus thinning is not required.
#As thinning does not seem necessary as much, it is more computationally efficient to produce indepent samples compared to part 1.
#Block is thus one of the ways in which we can improve MCMC samplers.