#------------------------------------------------
# coppie di valori secondo una normale bivariata
#
# GdA 3/3/09
#-------------------------------------------------

mx=0.4
my=2
sx=1
sy=0.5
rho=0.6

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

library(hacks)  # si carica la libreria

xy <- rbnorm(10000, 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(10000, mu, V)

# 2) metodo basato sulla condizionale
#
# a) estraiamo le x
x <- rnorm(10000, mx, sx) 

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)
plot(x, y, pch=19, cex=0.1)
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)
# etc.

# 3) usiamo trasformazione di variabili
# z1 ~ N(0,1)
# z2 ~ N(0,1)
# v = cos(theta) * z1 + sin(theta) * z2  [ con cos(theta)=rho ]
z1 <- rnorm(10000)
z2 <- rnorm(10000)
v <- rho*z1 + sqrt(1-rho^2)*z2
x <- z1*sx + mx
y <-  v*sy + my
mean(x)
mean(y)
sd(x)
sd(y)
plot(x,y)