#------------------------------------------------
#-----------------------------------------------------
# coppie di valori secondo una normale bivariata
# usando apposite librerie R o semplicemente rnorm()
# generando le due variabili a cascata (prima x,
# quindi y condizionata dal valore di x)
#
# GdA 25/1/11
#-----------------------------------------------------

mx = 2.5
my = 1.5
sx = 0.5
sy = 0.5
rho= -0.8

nevts = 10000

#---------------------------------------------------------------
# 1.1) usa package hacks
# ( se non esiste va installato con 'install.packages("hacks")' )

library(hacks)  # si carica la libreria

xy <- rbnorm(nevts, sx, sy, rho, mx, my)

summary(xy)

mean(xy[,1])
mean(xy[,2])
sd(xy[,1])
sd(xy[,2])
var(xy)
cor(xy)
plot(xy)
plot(xy[,1])
hist(xy[,1],nc=50)
hist(xy[,2],nc=50)

# 1.2) usa package mnormt (multidimensionale!)
library(mnormt)
mu <- c(mx, my)
V  <- rbind(  c( sx^2,      rho*sx*sy)
             ,c( rho*sx*sy, sy^2     )
           )       
xy <- rmnorm(nevts, mu, V)
plot(xy)

#--------------------------------------------------------------
# 2) medoto basato sulla condizionale
#
# a) estraiamo le x
x <- rnorm(nevts, mx, sx) 
#  b) estraiamo le y condizionate
y <- numeric()
for (i in c(1:length(x))) {
  y[i] <- rnorm(1, my + rho*sy/sx*(x[i]-mx), sy*sqrt(1-rho^2) )
}

mean(x)
mean(y)
sd(x)
sd(y)
cor(x,y)

plot(x,y)
hist(x, nc=50)
hist(y, nc=50)

# se vogliamo mettere i dati in una matrice per usare
# le stesse istruzioni usate precedentemente:
xy <- cbind(x,y)
# quindi: 
cor(xy)
plot(xy)