## ## ## Multivariate normal and ellipses ## ## source("/R/fun.R") source("/R/Plot3d.R") library(mvtnorm) mu = c(0,0) Sig = matrix( c(1,.9,.9,1.5),2) Sig dens <- function( X1 , X2 , ...) dmvnorm( c(X1,X2), mu , Sig) dens.rad <- function( rad , Sig) 1/ ( 2*pi* sqrt(det(Sig ))) * exp(- rad^2/2) maxdens <- dmvnorm( mu, mu, Sig) Plot3d( density ~ X1 + X2, list(density=0,X1=0,X2=0), xlim = c(-4,4), zlim =c(-4,4), ylim = c(0,maxdens), size=0) Lines3d( y = 0, x = c(-4,4), z =0, col = 'blue3', lwd = 1.5) Lines3d( y = 0, z = c(-4,4), x =0, col = 'blue3', lwd = 1.5) Lines3d( y = c(0,maxdens), z = 0, x =0, col = 'blue3', lwd = 1.5) Text3d( y = 0, x = 0, z = 0, text = " 0", col = 'black', cex = 2.5 ) Lines3d( y = 0, z = c(-4,4), x = 1, col = 'blue3') Lines3d( y = 0, z = c(-4,4), x =-1, col = 'blue3') Lines3d( y = 0, z = 1, x =c(-4,4), col = 'blue3') Lines3d( y = 0, z = -1, x =c(-4,4), col = 'blue3') Fit3d( dens, col = 'grey1', fill = F) for ( rad in seq(0,3,1)) { Lines3d( y = dens.rad( rad, Sig), xz = ell(mu,Sig, rad) , col = 'red', lwd = 2) } zform <- function( x ) paste( round(100*x,2),"%", sep = '') Text3d( y = dens.rad( 1, Sig), xz = ell(mu,Sig, 1)[1,,drop=F] , text = zform( pchisq( 1^2, df = 2)),col = 'red2', cex = 2) Text3d( y = dens.rad( 2, Sig), xz = ell(mu,Sig, 2)[1,,drop=F] , text = zform( pchisq( 2^2, df = 2)),col = 'red2', cex = 2) Text3d( y = dens.rad( 3, Sig), xz = ell(mu,Sig, 3)[1,,drop=F] , text = zform( pchisq( 3^2, df = 2)),col = 'red2', cex = 2)