#metroplis example: bivariate normal density with bivariate normal jumping kernel
# based on Francesca Dominici's code
metropolis <- function(theta,mu,Sigma){
  thetastar <- simulate.multnorm(theta,Sigma)
  Prec <- solve(Sigma)
  ratio <-  exp(-1/2*(t(thetastar-mu )%*%Prec%*%(thetastar-mu) -
			 t(theta-mu)%*%Prec%*%(theta-mu)))
  u <-  runif(1)
  if(u <= ratio) {theta <- thetastar}
  test <-  (u <= ratio)
  return(theta)
}

iterate <- function(theta0,NN){
  theta <- matrix(NA,2,NN)
  test <- NULL
  theta[,1] <- theta0
  for(n in 2:(NN-1)){
    theta[,n] <- metropolis(theta=theta[,n-1],mu=c(0,0),Sigma=diag(1,2))}
  return(theta)  
}

simulate.multnorm <- function(mu0,Q,M=1)
{
  s <- length(mu0)
  Z <-  matrix(rnorm(s*M),nrow=s)
  V <- t(chol(Q) )
  x <- (V%*%Z)+mu0
  return(x)
}



plotmcmc <- function(NN){
  par(oma=c(0,0,2,0))
  I <- iterate(theta0=c(0,0),NN) 
  I1 <- iterate(theta0=c(3,3),NN)
  I2 <- iterate(theta0=c(3,-3),NN) 
  I3 <- iterate(theta0=c(-3,3),NN)
  I4 <- iterate(theta0=c(-3,-3),NN)
  plot(I[1,],I[2,],type="l",xlim=c(-4,4), ylim=c(-4,4),xlab="theta1",ylab="theta2") 
  points(0,0,pch=0,cex=4)
  lines(I1[1,],I1[2,],lty=2) 
  points(3,3,pch=0,cex=4)
 lines(I2[1,],I2[2,],lty=2) 
  points(3,-3,pch=0,cex=4)
  lines(I3[1,],I3[2,],lty=2) 
  points(-3,3,pch=0,cex=4)
  lines(I4[1,],I4[2,],lty=2) 
  points(-3,-3,pch=0,cex=4)
  par(oma=c(0,0,0,0)) 

}