#Gibbs example: bivariate normal distribution:example pag:328

KK _ data()
II _ initialize()
II _ gibbs(KK,II)
mc _ iterate(N=100)

data _ function(){
  y1 _ y2 _ 0
  rho_ .81
  return(y1,y2,rho)
}


initialize_function(){
theta1 _ theta2 _ 0
return(theta1,theta2)}


gibbs_function(KK,II){
  theta1 _ II$theta1
  theta2 _ II$theta2
  theta1 _   rnorm(1,KK$y1 + KK$rho*(theta2-KK$y2),sqrt(1-KK$rho^2))
  theta2 _ rnorm(1,KK$y2 + KK$rho*(theta1-KK$y1),sqrt(1-KK$rho^2))
  return(theta1,theta2)
}


iterate _ function(NN=1000){
  theta1 _ theta2 _ NULL
  for( n in 1:NN){
    II _ gibbs(KK,II)
    theta1[n] _ II$theta1
    theta2[n] _ II$theta2
  }
  return(theta1,theta2)
}


##################
simu1_ function(NN=1000) {
  theta _ matrix(NA,2,NN)
  for( nn in 1:NN){
    theta[,nn] _ simulate.multnorm(c(0,0),cbind(c(1,0.81),c(0.81,1)))
  }   
  return(theta)
}





gibbs.iterate_function(a=0,b=0,NN){
  theta1_NULL
  theta2_NULL
  theta1[1]_a
  theta2[1]_b
  for(n in 2:(NN-1)){
    theta1[n]_gibbs(theta1=theta1[n-1],theta2=theta2[n-1],y1=0,y2=0,rho=.81)$theta1
    theta2[n]_gibbs(theta1=theta1[n],theta2=theta2[n-1],y1=0,y2=0,rho=.81)$theta2
  }
  return(theta1,theta2)
}


I_gibbs.iterate(0,0,100)






plotgibbs_function(w=0,NN){
  if (w == 1) postscript("/home/biostats/fdominic/course/mcmc.gibbs.ps")
  par(oma=c(0,0,2,0))
  I_gibbs.iterate( a=0,b=0,NN) 
  I1_gibbs.iterate( a=3,b=3,NN)
  I2_gibbs.iterate( a=3,b=-3,NN) 
  I3_gibbs.iterate( a=-3,b=3,NN)
  I4_gibbs.iterate( a=-3,b=-3,NN)
 
plot(I$theta1,I$theta2,type="l",xlim=c(-4,4), ylim=c(-4,4),xlab="theta1",ylab="theta2") 
  points(0,0,pch=0,cex=4)
  lines(I1$theta1,I1$theta2,lty=2) 
  points(3,3,pch=0,cex=4)
 lines(I2$theta1,I2$theta2,lty=2) 
  points(3,-3,pch=0,cex=4)
  lines(I3$theta1,I3$theta2,lty=2) 
  points(-3,3,pch=0,cex=4)
  lines(I4$theta1,I4$theta2,lty=2) 
  points(-3,-3,pch=0,cex=4)
  par(oma=c(0,0,0,0)) 
  if (w == 1) dev.off()
}