gibbs.update=function(start=c(1,1), n=10, rho=0.5, pause=T) {
# use pause=T for small n, to display step by step pressing return
# in between.
# use pause-F to run stragiht through and plot the pictures only
# at the end

# set up vectors for the samples
  x=rep(NA,n)
  y=rep(NA,n)
  x[1]=start[1]
  y[1]=start[2]

  #Set up plots
  xx=seq(-3,3,len=20)
  yy=seq(-3,3,len=20)
  z = outer(xx,yy,f,rho)
  if (pause) {
    layout(matrix(c(0,1,1,0,0,1,1,0,2,2,3,3,2,2,3,3),4,4,byrow=T))
    #plot contours 
    contour(xx,yy,z,col="blue",levels=c(0.01,0.02,0.05,0.1,0.14),
            drawlabels=F)
    title(xlab="x",ylab="y")
    #mark starting point:
    points(x[1],y[1],pch=4,col="red")
    #Time series plot so far:
    # X:
    plot(c(0,n),c(min(start[1],-3),max(start[1],3)),
         xlab="iterations",ylab="x",type="n")
    lines(c(0,n),c(1.96,1.96),lty=2,col="blue")
    lines(c(0,n),-c(1.96,1.96),lty=2,col="blue")
    # Y:
    plot(c(0,n),c(min(start[2],-3),max(start[2],3)),
         xlab="iterations",ylab="y",type="n")
    lines(c(0,n),c(1.96,1.96),lty=2,col="blue")
    lines(c(0,n),-c(1.96,1.96),lty=2,col="blue")
    # Wait for key press
    readLines(n=1)
  }

    
  #Main loop
  for (t in 2:n) {

    # sample x from p(x|y=y[t-1]):
    x[t]=rnorm(1, rho* y[t-1], sqrt(1-rho^2))
    # sample y from p(y|=x[t]):
    y[t]=rnorm(1, rho* x[t], sqrt(1-rho^2))

    if (pause || (t==n)) {
    
      #plot contours (use xx,yy and z from before)
      contour(xx,yy,z ,col="blue",levels=c(0.01,0.02,0.05,0.1,0.14),
              drawlabels=F)
      title(xlab="x",ylab="y")
      #mark starting point:
      points(x[1],y[1],pch=4,col="red")
      #draw path of the chain
      lines(c(x[1],rep(x[2:t],rep(2,t-1))),
            c(rep(y[1:(t-1)],rep(2,t-1)),y[t]),col="red")

      #Time series plot so far:
      # X:
      plot(c(0,n),c(min(start[1],-3),max(start[1],3)),
           xlab="iterations",ylab="x",type="n")
      lines(1:t,x[1:t],col="red")
      lines(c(0,n),c(1.96,1.96),lty=2,col="blue")
      lines(c(0,n),-c(1.96,1.96),lty=2,col="blue")
      # Y:
      plot(c(0,n),c(min(start[2],-3),max(start[2],3)),
           xlab="iterations",ylab="y",type="n")
      lines(1:t,y[1:t],col="red")
      lines(c(0,n),c(1.96,1.96),lty=2,col="blue")
      lines(c(0,n),-c(1.96,1.96),lty=2,col="blue")
      
    }

    if ((pause)  && (t<n)) {
      # Wait for key press
      readLines(n=1)
    }


  #end main loop
  }
}