#   6.1   Gibbs sampler example


# (x,y) ~ N(theta.x,theta.y, (1,.8,.8,1))

#  x = 0, y = 0

x <- 0
y <- 0

theta.x <- rep(-5,50)
theta.y <- rep(-5,50)

for(i in 2:50)
{
	theta.x[i] <- rnorm(1,mean = x + 0.8*(theta.y[i-1] - y),sd = sqrt(1 - 0.8^2))
	theta.y[i] <- rnorm(1,mean = y + 0.8*(theta.x[i]   - x),sd = sqrt(1 - 0.8^2))
}	

par(mfcol=c(1,2))
plot(theta.x[1],theta.y[1],col=3,xlim=c(-5,5),ylim=c(-5,5),pch="+",cex=2)
lines(theta.x[1:50],theta.y[1:50],col=3)

theta.x <- rep(-5,50)
theta.y <- rep(5,50)

for(i in 2:50)
{
	theta.x[i] <- rnorm(1,mean = x + 0.8*(theta.y[i-1] - y),sd = sqrt(1 - 0.8^2))
	theta.y[i] <- rnorm(1,mean = y + 0.8*(theta.x[i]   - x),sd = sqrt(1 - 0.8^2))
}	

points(theta.x[1],theta.y[1],col=4,pch="+",cex=2)
lines(theta.x[1:50],theta.y[1:50],col=4)


theta.x <- rep(5, 12000)
theta.y <- rep(-5,12000)

for(i in 2:12000)
{
	theta.x[i] <- rnorm(1,mean = x + 0.8*(theta.y[i-1] - y),sd = sqrt(1 - 0.8^2))
	theta.y[i] <- rnorm(1,mean = y + 0.8*(theta.x[i]   - x),sd = sqrt(1 - 0.8^2))
}	

points(theta.x[1],theta.y[1],col=2,pch="+",cex=2)
lines(theta.x[1:50],theta.y[1:50],col=2)


plot(theta.x[10000:12000],theta.y[10000:12000],col=3,xlim=c(-5,5),ylim=c(-5,5),cex=1)

summary(theta.x[10000:120000])

hist(theta.x[10000:120000], freq = FALSE)
curve(dnorm(x), col = 2, lty = 2, lwd = 2, add = TRUE)

qqnorm(theta.x[10000:120000])
abline(0,1,col=2,lty=2,lwd=2)

 

#   6.2  Metropolis sampler


# (x,y) ~ N(theta.x,theta.y, (1,.8,.8,1))
#  x = 0, y = 0

library(mvtnorm)

Mean  <- c(0,0)
Sigma <- matrix(c(1,.8,.8,1),nrow=2,byrow=T)

theta.x <- rep(-5,10000)
theta.y <- rep(-5,10000)

for(i in 2:10000)
{
	theta.star.x <- theta.x[i-1] + rnorm(1,mean=0,sd=.2)
	theta.star.y <- theta.y[i-1] + rnorm(1,mean=0,sd=.2)

	d.im1 <- dmvnorm(c(theta.x[i-1],theta.y[i-1]), mean=Mean, sigma=Sigma, log=T)
	d.star<- dmvnorm(c(theta.star.x,theta.star.y), mean=Mean, sigma=Sigma, log=T)

	r <- min(c(1,exp(d.star - d.im1)))
	if(runif(1) < r)
	{
		theta.x[i] <- theta.star.x 
		theta.y[i] <- theta.star.y 
	}
	else
	{
		theta.x[i] <- theta.x[i-1] 
		theta.y[i] <- theta.y[i-1]
	}
}

par(mfcol=c(1,2))
plot(theta.x[1],theta.y[1],col=1,xlim=c(-5,5),ylim=c(-5,5),pch="+",cex=2)
lines(theta.x[1:500],theta.y[1:500],col=1)

plot(theta.x[1],theta.y[1],col=1,xlim=c(-5,5),ylim=c(-5,5),pch="+",cex=2)
lines(theta.x[1:10000],theta.y[1:10000],col=1)

hist(theta.x[1000:10000], freq = FALSE)
curve(dnorm(x), col = 2, lty = 2, lwd = 2, add = TRUE)

qqnorm(theta.x[1000:10000])
abline(0,1,col=2,lty=2,lwd=2)




#  6.2.1 Application of Metropolis-Hastings algorithm on a discrete MC

#  a. The random walks 

P.RW <- array(0.1,dim=c(6,6)) + diag(rep(.4,6))

P.RW <- rbind(c(.5,.1,.1,.1,.1,.1),array(1/6,dim=c(4,6)),c(.1,.1,.1,.1,.1,.5))

P.RW <- matrix(c(0,1,0,0,0,0,.5,0,.5,0,0,0,0,.5,0,.5,0,0,0,0,.5,0,.5,0,0,0,0,.5,0,.5,0,0,0,0,1,0),nrow=6,ncol=6,byrow=TRUE)

P.RW

s <- rep(NA,50000)
s[1] <- 3
for (j in 2:50000)
   s[j] <- sample(1:6,size=1,prob=P.RW[s[j-1],])

m <- c(500,2000,8000,50000)
for (i in 1:4)
   print(table(s[1:m[i]])/m[i])


P.RW%*%P.RW
P.RW%*%P.RW%*%P.RW%*%P.RW
P.RW%*%P.RW%*%P.RW%*%P.RW%*%P.RW%*%P.RW%*%P.RW%*%P.RW

P.tmp	<- P.RW
for(i in 1:1000) P.tmp <- P.tmp %*% P.RW
P.tmp

P.eigen	<- P.tmp[1,]
P.eigen%*%P.RW



# b. Use Metropolis-Hastings MCMC algorithm to create MC with 
#    uniform stationary distribution from RW 

P.MH <- matrix(0,ncol=6,nrow=6)

for(i in 1:6)
{
	p.stay <- P.RW[i,i]
	for(j in (1:6)[-i])
		if(P.RW[i,j] > 0)
		{	
			p.ACC <- min(c(1,(1 / P.RW[i,j]) / (1 / P.RW[j,i])))
			P.MH[i,j] <- P.RW[i,j] * p.ACC
			p.stay    <- p.stay + P.RW[i,j] * (1 - p.ACC)
		}
		
	P.MH[i,i] <- p.stay
}

P.MH
rbind(rep(1/6,6)) %*% P.MH


# c. Use Metropolis-Hastings MCMC algorithm to create MC with Bin(5,.3) stationary distribution from RW 

P.MH <- matrix(0,ncol=6,nrow=6)
(w.bin <- dbinom(0:5,5,.3))

for(i in 1:6)
{
	p.stay <- P.RW[i,i]
	for(j in (1:6)[-i])
		if(P.RW[i,j] > 0)
		{	
			p.ACC <- min(c(1,(w.bin[j] / P.RW[i,j]) / (w.bin[i] / P.RW[j,i])))
			P.MH[i,j] <- P.RW[i,j] * p.ACC
			p.stay    <- p.stay + P.RW[i,j] * (1 - p.ACC)
		}
		
	P.MH[i,i] <- p.stay
}

P.MH
P.MH%*%P.MH


s <- rep(NA,50000)
s[1] <- 3
for (j in 2:50000)
	 s[j] <- sample(1:6,size=1,prob=P.MH[s[j-1],])

m <- c(500,2000,8000,50000)
for (i in 1:4)
   print(table(s[1:m[i]])/m[i])

w.bin %*% P.MH