# attenzione: in caso di errore nel caricamento delle
# librerie vanno installati i pacchetti relativi
# -> ed es. install.packages("mnormt")

library(mnormt)
library(rgl)

# parametri della gausiana bivariata
mx=0.4
my=2
sx=1
sy=0.5
rho=0.6

# vettore dei valori medi
mu <- c(mx, my)

# matrice di covarianza
V  <- rbind(  c( sx^2,      rho*sx*sy)
             ,c( rho*sx*sy, sy^2     )
           )       

# griglia nella quale calcolare la funzione
x <- seq(mx - 3*sx, mx + 3*sx, len=51)
y <- seq(my - 3*sy, my + 3*sy, len=51)

z <- outer(x,y) # inizializziamo con una matrice che contiene
                # la somma (tanto per metterci qualcosa: i valori
                # verranno riscritti nel loop)

# matrice dei valori della pdf bidimensionale
for(i in 1:length(x)) {
  for(j in 1:length(y)) {
    z[i,j] <- dmnorm(c(x[i],y[j]), mu, V)   
  }
}

# plot 3-D
persp3d(x, y, z, col='gray', xlab="x", ylab="y", zlab="f(x,y)")

# -> con mouse il plot si può ingrandire, quindi ruotare a piacere