#
# Bivariate normal explored by a Gibbs Sampler
# (pure academic exercise, even a bit pedantic...)
#
#  GdA April 2016 (rev. May 2020)
#----------------------------------------------------

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

n = 10000

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

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

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

n.min = 1
take <- n.min:n   # select the range to analyse

print(summary(x1[take]))
print(summary(x2[take]))

print(cor(x1[take],x2[take]))

#  insert 
plot(x1[take], col='blue', cex=0.8)
pause()
plot(x2[take],  col='blue',  cex=0.8)
pause()
plot(x2[take],x1[take],  col='blue', cex=0.8, asp=1)