# come gaussiana_bivariata_3D.R
# con i valori della pdf calcolati direttamente
# (e più velocemente) in outer(), invece che nel
# loop bidimensionale

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)

f <- function(x,y) dmnorm(cbind(x,y), mu, V)
z <- outer(x,y, f)

# Nota: questo si può fare con dmnorm(),
# ma non con dmgauss(), in quanto quest'ultima, nella sua
# semplicità, non accetta input matriciale.

# 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