#bivariate normal density f<-function (x, y = x, a=0,b=0,s1=1,s2=2,rho = 0) { xoy = ((x-a)^2/s1^2 - 2 * rho * (x-a) * (y-b)/(s1*s2) + (y-b)^2/s2^2)/(2 * (1 - rho^2)) density = exp(-xoy)/(2 * pi * s1*s2*sqrt(1 - rho^2)) density } x<-y<-seq(-4,4,length=100) z<-outer(x,y,function(x,y)f(x,y)) #static plot persp(x,y,z,theta=45,phi=25,col='lightblue') #3d interactive plot require(rgl) persp3d(x,y,z,theta=45,phi=25,col='lightblue') M <- par3d("userMatrix") bg3d("white") play3d(par3dinterp(userMatrix=list(M,rotate3d(M, angle=pi,x=1,y=0,z=0))), duration=10) #different rho for(i in seq(-0.9, 0.9, by=.1) ) { z=outer(x,y,function(x,y)f(x,y,a=0,b=0,s1=1,s2=1,rho=i)) persp(x,y,z,theta=0,phi=25,col='lightblue',main=paste("rho=",as.character(i))) locator(1) } #different s1 for(s1 in seq(0.1, 5, by=.5) ) { z=outer(x,y,function(x,y)f(x,y,a=0,b=0,s1=s1,s2=1,rho=0)) persp(x,y,z,theta=45,phi=25,col='lightblue',main=paste("s1=",as.character(s1))) locator(1) } ##contour plot par(mfrow=c(2,2)) x<-y<-seq(-4,4,length=100) z<-outer(x,y,function(x,y)f(x,y)) contour(z,main="independent") z<-outer(x,y,function(x,y)f(x,y,rho=0.5)) contour(z,main="rho=0.5") z<-outer(x,y,function(x,y)f(x,y,s1=4,rho=0.5)) contour(z,main="s1=4 and rho=0.5") z<-outer(x,y,function(x,y)f(x,y,s1=4,rho=-0.9)) contour(z,main="s1=4 and rho=-0.9") ##simulate from a multivariate normal distribution mu<-c(0,1) Sigma<-matrix(c(1,0.5,0.5,1),2,2) n<-1000 library(MASS) biv<-mvrnorm(n,mu,Sigma) colnames(biv)<-c("X","Y") #alternative, use rnorm sig.eigen<-eigen(Sigma) Sig<-sig.eigen\$vectors%*%diag(sqrt(sig.eigen\$values))%*%t(sig.eigen\$vectors) biv2<-mu+t(matrix(rnorm(2*n),ncol=2)%*%Sig) biv2<-t(biv2) par(mfrow=c(1,2)) hist(biv2[,1],main="Histogram of X") hist(biv2[,2],main="Histogram of Y") install.packages("mixtools") library(mixtools) plot(biv) ellipse(mu<-colMeans(biv),sigma<-cov(biv),alpha=.05,col='red') points(t(mu),col='red',pch=19) plot(biv2) ellipse(mu<-colMeans(biv2),sigma<-cov(biv2),alpha=.05,col='red') points(t(mu),col='red',pch=19) #bivariate normal distribution install.packages("fMultivar") library(fMultivar) rnorm2d(n,rho=0) pnorm2d(x=0,y=0,rho=0) dnorm2d(x=0,y=0,rho=0) install.packages("mvtnorm") library(mvtnorm) m <- 3 sigma <- diag(3) sigma[2,1] <- 3/5 sigma[3,1] <- 1/3 sigma[3,2] <- 11/15 pmvnorm(lower=rep(-Inf, m), upper=c(1,4,2), mean=rep(0, m), corr=sigma) qmvnorm(0.95, sigma = diag(2), tail = "both") x <- rmvnorm(n=500, mean=c(1,2,1), sigma=sigma) colMeans(x) var(x) ### bivariate normal density with matlab http://personal.kenyon.edu/hartlaub/MellonProject/Bivariate2.html