#
# Normale bivariata campionata con un Gibbs Sampler
# (codice 'pedante' per far capire bene i passaggi)
#
#  GdA 23/4/09
#----------------------------------------------------

n = 10000

mu1 = 1
s1  = 1
mu2 = 2
s2  = 3
rho = 0.6

y1 <- numeric()
y2 <- numeric()

x1.in = -10
x2.in = -10

y1[1] = x1.in
y2[1] = x2.in

x1 = y1[1]   # although not used ..
x2 = y2[1]
for (i in 2:n) {
  y1[i] <- rnorm(1, mu1 + s1/s2*rho*(x2-mu2), s1*sqrt(1-rho^2))
  y2[i] <- rnorm(1, mu2 + s2/s1*rho*(y1[i]-mu1), s2*sqrt(1-rho^2))
  x1 = y1[i]  # although not used
  x2 = y2[i]
}