#
# Bivariate normal explored by a Gibbs Sampler
#
#  GdA, May 2020
#----------------------------------------------------

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

n = 10000

mu1 = 1           # parameters
s1  = 2
mu2 = 1
s2  = 2
rho = 0.5

x1 = x2 = numeric(n)  # empty vectors
x1[1] = 10
x2[1] = 10

for (i in 2:n) {  
  x1[i] <- rnorm(1, mu1 + s1/s2*rho*(x2[i-1]-mu2), s1*sqrt(1-rho^2))    # x1 given x2
  x2[i] <- rnorm(1, mu2 + s2/s1*rho*(x1[i]-mu1), s2*sqrt(1-rho^2))      # x2 given x1
}

plot(NULL, NULL, xlab='', ylab='', xlim=c(-10,10), ylim=c(-10,10))
n.max = min(500, n)
points(x1[1], x2[1], pch=19, col='blue')
for (i in 2:n.max) {
  lines(c(x1[i-1],x1[i]), c(x2[i-1],x2[i]), col='cyan')
  points(x1[i], x2[i], pch=19, col='blue')
  pause()
}