# come gaussiana_bivariata_3D.R, ma fa uso di # dmgauss() autocostruita al posto di dmnorm() # controlla vettore dei valori attesi e matrice di covarianza check.V <- function(mu, V) { n <- length(mu) # controlla dimensioni di V if ( sum(dim(V) != n) ) return(FALSE) # controlla che V sia definita positiva (nessun autovettore <= 0) if( sum( eigen(V)$values <= 0 ) > 0) return(FALSE) return(TRUE) } # pdf multidimensionale (x è un vettore!) dmgauss <- function(x, mu, V) { # controlla lunghezze dei vettori mu e x # [eventuali altri controlli possono essere fatti, una tantum # dalla funzione check.V() ] n <- length(mu) if ( length(x) != n ) { cat( sprintf("Errore: i vettori x e mu devono essere della stessa lunghezza!\n")) return(NULL) } # pdf (2*pi)^(-n/2) * det(V)^(-1/2) * exp( - t(x-mu) %*% solve(V) %*% (x-mu) / 2 ) } # libreria per la grafica 3D 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) # variante (rispetto a 'outer' usata in gaussiana_bivariata_3D.R) # per costruire una matrice vuota delle dimensioni giuste z <- matrix( rep(0, length(x)*length(y)), c(length(x), length(y)) ) # matrice dei valori della pdf bidimensionale if( check.V(mu, V) ) { for(i in 1:length(x)) { for(j in 1:length(y)) z[i,j] <- dmgauss(c(x[i],y[j]), mu, V) } } else stop ("Problemi con la matrice di covarianza!") # plot 3-D persp3d(x, y, z, col='gray', xlab="x", ylab="y", zlab="f(x,y)") # -> con il mouse il plot si può ingrandire, quindi ruotare a piacere #------------------------------------------------------------ # function pausa pausa <- function() { cat ("\n >> Guarda il plot e dai enter per continuare\n") scan() } pausa() # marginali f(x|y) per tre valori di y, definiti in yc yc <- c(my-sy, my, my+sy) # ridefiniamo x su un intervallo più ampio x <- seq(mx - 5*sx, mx + 5*sx, len=201) plot(x, dnorm(x, mx, sx), ty='l', col='blue', lwd=1.5, ylim=c(0, 1.5*max( dnorm(x, mx, sx) ) ), ylab='f(x)' ) # plottiamo le condizionali "fatte a mano" for (y in yc) { fxnn <- rep(0, length(x)) # fxnn: 'f(x)' non normalizzata for (i in 1:length(x)) fxnn[i] <- dmgauss(c(x[i], y), mu, V) norma <- sum(fxnn) * (x[2]-x[1]) fx <- fxnn / norma points(x, fx, ty='l', lty=2, col='red') } pausa() # plottiamo le condizionali usando la regola generale # di shift e strizzamento for (y in yc) { points(x, dnorm(x, mx + rho*sx/sy*(y-my), sx*sqrt(1-rho^2)), ty='l', lty=3, col='magenta') }