# # 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] }