for (it in 1:10000) { xhat <- rnorm(1, sd=.2) + x0 yhat <- rnorm(1, sd=.2) + y0 if( hat(xhat,yhat)/hat(x0,y0) > runif(1)) { points(xhat, yhat, pch=".", cex=5, col="red") x0 <- xhat y0 <- yhat } else { points(x0 + .01*runif(1), y0 + .01*runif(1), pch=".",cex=5, col="red") } xm[it,] <- c(x0,y0) }