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