#------------------------------------------------
# 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) medoto 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)
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)