gibbs1=function( n=1000 , rho=0.99 , start1=-1 , start2=-1 )
 {
  x1 <- c(1:n)
  x2 <- c(1:n)
  x1[1] <- start1
  x2[1] <- start2
  par( mfrow=c( 2,2 ) )
  for ( i in 2:n )
   {
    x1[i] <- rnorm( 1 , x2[i-1] * rho , 1-rho^2 )
    x2[i] <- rnorm( 1, x1[i] * rho , 1 - rho^2 )

  plot( x1[1:i] , type="l" , xlab="Iteration" ,xlim=c(1,n),ylab="x1")
  plot( x2[1:i] , type="l" , xlab="Iteration",xlim=c(1,n),ylab="x2" )
  plot( x1[1:i] , x2[1:i] , type="p" , pch="." ,xlab="x1",ylab="x2")
  plot( x1[1:i] , x2[1:i] , type="l" ,col="red",xlab="x1",ylab="x2")
 }
}