#-------------------------------------------------------
#
#  first example of MCMC chain
#
#  GdA 1/4/2016
#-------------------------------------------------------

pause <- function() { cat ("\n >> Press Enter to continue\n"); scan() }

#--------------------------------------------------------
# proposal on a square of side d ------------------------
prop.sq <- function(xy.old, d)  xy.old + runif(2,-d, d)

n  <- 1000    
dt <- 0.002     # time between moves: 0 -> manual
d  <- 0.5

xy <- c(0,0)
lim <- c(-1,1) # acceptance bounds
pl <- c(-3,3)  # limits for plotting

old.mar = par("mar")
par(mar=c(4,3,0.1,0.1))
plot(xy[1], xy[2], xlim=pl, ylim=pl, pch=19, cex=0.5,
     xlab='x', ylab='y', asp=0.5, col='red')
for(i in 1:n) {
  xy.old <- xy  
  xy <- prop.sq(xy.old, d)
  lines(c(xy.old[1], xy[1]), c(xy.old[2], xy[2]), col='blue')
  points(xy[1], xy[2],  pch=19, col='blue', cex=0.5)
  if(dt > 0) {
    Sys.sleep (dt) 
  } else {
    pause()  
  }
}

par(mar=old.mar)