#-------------------------------------------------------
#
#  MCMC on a square
#
#  GdA 1/4/2016
#-------------------------------------------------------

# functions -------------------------------------------------------------
pause <- function() { cat ("\n >> Press Enter to continue\n"); scan() }
delay <- function(dt)  if(dt > 0) {  Sys.sleep (dt) } else { pause() }
in.range <- function(x, r) x >= r[1] &  x <= r[2]

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

n  <- 10000
dt <- 0.0001     # time between moves: 0 -> manual
d  <- 2

multi = TRUE   # if TRUE, we record several times a point if acceptance fails
draw.lines = TRUE

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

x <- rep(0, n)  # vectors to store the history
x[1] <- xy[1]    
y <- rep(0, n) 
y[1] <- xy[2]

plot(xy[1], xy[2], xlim=pl, ylim=pl, pch=19, cex=0.5,
     xlab='x', ylab='y', asp=1, col='red')
if(draw.lines) lines(c(lim, rev(lim), lim[1]), c(lim[1], lim, rev(lim)), col='green')

n.yes = 1
for(i in 2:n) {
  xy.old <- xy  
  xy <- prop.sq(xy.old, d)
  # accept ?
  if(in.range(xy[1],lim) & in.range(xy[2],lim)) {
    n.yes <- n.yes + 1  
   if(draw.lines) lines(c(xy.old[1], xy[1]), c(xy.old[2], xy[2]), col='cyan')
    points(xy[1], xy[2],  pch=19, col='blue', cex=0.5)
    x[n.yes] <- xy[1]
    y[n.yes] <- xy[2] 
  } else {
    lines(c(xy.old[1], xy[1]), c(xy.old[2], xy[2]), lty=2, col='red')
    xy <- xy.old  # reset old point
    if (multi) { # record many times the same point if rejected
      n.yes <- n.yes + 1
      x[n.yes] <- xy[1]
      y[n.yes] <- xy[2]
    } 
  }    
  delay(dt)
}
x <- x[1:n.yes]
y <- y[1:n.yes]

# other plots and histograms
# plot(x,y, pch=19, cex=0.5, col='blue')
# hist(x, nc=50,col='cyan')
# hist(y, nc=50,col='cyan')

par(mfrow = c(1,2))
 hist(x, nc=50,col='cyan', main=NULL)
 hist(y, nc=50,col='cyan', main=NULL)
par(mfrow = c(1,1))