################################################################################ #### Demo of Bivariate Normal Distribution #### ################################################################################ #### By Jimin Ding, 02/8/2017 library(rgl) ## For interactive 3D plots library(mvtnorm) ## For multivariate normal distributions ?dmvnorm ## One may type the function name to see how the author define the fuction in details. dmvnorm ## Now let's draw the density for the bivariate normal distribution with ## E(X)=0, E(Y)=1, Var(X)=1, Var(Y)=4, Corr(X,Y)=0.5 CovXY=matrix(c(1,1,1,4),nrow=2) x=seq(from=-3,to=3,by=0.1) y=seq(from=-7,to=7,by=0.1) xy=cbind(x=rep(x,each=length(y)),y=rep(y,times=length(x))) fxy=dmvnorm(xy,mean=c(0,1),sigma=CovXY) z=matrix(fxy,byrow=TRUE,nrow=length(x),ncol=length(y)) library(rgl) persp3d(x,y,z,col='steelblue',alpha=0.7) grid3d(c("x", "y", "z")) planes3d(a=0,b=0,c=-1,d=0.08,alpha=0.5, col='green') planes3d(a=0,b=0,c=-1,d=0.06,alpha=0.5, col='green') planes3d(a=0,b=0,c=-1,d=0.04,alpha=0.5, col='green') clipplanes3d(a=0,b=0,c=-1,d=0.08) clipplanes3d(a=0,b=0,c=-1,d=0.06) clipplanes3d(a=0,b=0,c=-1,d=0.04) ## Generate a random sample from above Binormal distributions. XYi=data.frame(rmvnorm(1000,mean=c(0,1),sigma=CovXY)) names(XYi)=c('x','y') plot(XYi[,1],XYi[,2],pch=20,col=rgb(0.2,0.2,1,alpha=0.5), xlab='x',ylab='y',main='Bivariate Normal Random Variables') points(0,1,col='red',pch=19) contour(x,y,z,add=T,nlevel=4,labcex=1,lwd=2) source('http://www.math.wustl.edu/~jmding/math3200/hist3d.R') hist3d(XYi\$x,XYi\$y) ## Marginal Distributions ## We push the observations to X and Y-axes. ## (The points are jiggled a little to see the frequency.) library(ggplot2) scatter <- qplot(x,y, data=XYi) + scale_x_continuous(limits=c(min(x),max(x))) + scale_y_continuous(limits=c(min(y),max(y))) + geom_rug(col=rgb(.5,0,0,alpha=.2)) scatter ## We see the darkest part of x-axis is around 0 and that of y-axis is around 1. hist(XYi\$x) hist(XYi\$y) ## Now we project the all curves to x plot(x,z[,80],type='l',ylab='f(x,fixed y)') for (i in 1:dim(z)[2]){ lines(x,z[,i]) } ## If we aggregate all curves in above figure, we will have the density of x, ## which is can be also viewed from the histgoram of x. ##### Conditional Distributions. ## We consider f(x,1), so the density at y=1 persp3d(x,y,z,col='steelblue',alpha=0.7) planes3d(a=0,b=-1,c=0,d=1,alpha=0.5, col='green') clipplanes3d(a=0,b=-1,c=0,d=1) ## if we scale this curve by the probability that y at(around) y=1, ## then we have the conditional density f(x|y=1)