# Rcode for the class # Chapter 2 ################################## # Section 2.1 Bivariate normal plots library("fMultivar") library("MASS") # First make a scatter plot of bivariate Gaussian # the iid standard normal case Sigma = matrix(c(1,0,0,1),ncol = 2) # The vector mean in this example is (1,2) # Below we draw 200 times a bivariate normal # distribution with mean (1,2) and variance Sigma xtemp = mvrnorm(500,c(0,0),Sigma) # Plotting the uncorrelated bivariate data plot(xtemp, main="",xlab="Z1",ylab="Z2", ylim=c(-5,5),xlim =c(-5,5), col="blue",asp=1,pch=16) # Next make a transformation lambda1 = sqrt(0.9) lambda2 = sqrt(0.1) lambda1 = sqrt(0.1) lambda2 = sqrt(0.9) AT = matrix(c(lambda1,lambda2,lambda1,-lambda2),ncol=2,byrow=T) Txtemp = t(AT%*%t(xtemp)) plot(Txtemp, main="",xlab="X1",ylab="X2", ylim=c(-5,5),xlim =c(-5,5), col="red",asp=1,pch=16) # Drawing # Correlated normal data # Here we define the variance matrix. Sigma = matrix(c(1,0.5,0.5,1),ncol = 2) # The vector mean in this example is (1,2) # Below we draw 200 times a bivariate normal # distribution with mean (1,2) and variance Sigma xtemp = mvrnorm(500,c(0,0),Sigma) # Plotting the bivariate data plot(xtemp, main="",xlab="X",ylab="Y", ylim=c(-5,5),xlim =c(-5,5), col="blue",asp=1,pch=16) # Note the option asp gives the same scale of the x and y-axis. ## dnorm2d - # Bivariate Normal Density: x <- (-40:40)/10 X <- grid2d(x) # Below we make a plot when the correlation is 0.5. z <- dnorm2d(X$x, X$y, rho = 0.8) ZD <- list(x = x, y = x, z = matrix(z, ncol = length(x))) # Perspective Density Plot: #persp(ZD, theta = -40, phi = 30, col = "orange") par(mfrow=c(1,2)) persp(ZD, theta = -20, phi = 20, col = "orange", main="Bivariate Normal: cor = 0.8") persp(ZD, theta = 20, phi = 20, col = "orange", main="Bivariate Normal: cor = 0.8") par(mfrow=c(1,1)) contour(ZD, main="Bivariate Normal Density cor = 0.8",asp=1) # Below we make a plot when the correlation is 0. # Perspective Density Plot: #persp(ZD, theta = -40, phi = 30, col = "orange") #persp(ZD, theta = 40, phi = 30, col = "orange") z <- dnorm2d(X$x, X$y, rho = 0.0) ZD <- list(x = x, y = x, z = matrix(z, ncol = length(x))) # Perspective Density Plot: #persp(ZD, theta = -40, phi = 30, col = "orange") #persp(ZD, theta = 40, phi = 30, col = "orange") par(mfrow=c(1,2)) persp(ZD, theta = -20, phi = 20, col = "orange", main="Bivariate Normal: cor = 0") #persp(ZD, theta = 0, phi = 20, col = "orange", main="Bivariate Normal: cor = 0") persp(ZD, theta = 20, phi = 20, col = "orange", main="Bivariate Normal: cor = 0") par(mfrow=c(1,1)) contour(ZD, main="Bivariate Normal Density cor = 0",asp=1) # Below we make a plot when the correlation is 0. # Perspective Density Plot: #persp(ZD, theta = -40, phi = 30, col = "orange") #persp(ZD, theta = 40, phi = 30, col = "orange") z <- dnorm2d(X$x, X$y, rho = -0.8) ZD <- list(x = x, y = x, z = matrix(z, ncol = length(x))) # Perspective Density Plot: #persp(ZD, theta = -40, phi = 30, col = "orange") #persp(ZD, theta = 40, phi = 30, col = "orange") par(mfrow=c(1,2)) persp(ZD, theta = -20, phi = 20, col = "orange", main="Bivariate Normal: cor = -0.8") #persp(ZD, theta = 0, phi = 20, col = "orange", main="Bivariate Normal: cor = -0.8") persp(ZD, theta = 20, phi = 20, col = "orange", main="Bivariate Normal: cor = -0.8") par(mfrow=c(1,1)) contour(ZD, main="Bivariate Normal Density cor = -0.8",asp=1) ################################# # Section 2.4.2 Making orthogonal basis # n = 2 y1 = c(0,1/sqrt(2)) x1 = c(0,1/sqrt(2)) x2 = c(0,-1/sqrt(2)) y2 = c(0,1/sqrt(2)) plot(x1,y1, type='l', lwd=4, main="orthonormal basis (e1 and e2)", xlab="x-axis",ylab="y-axis", ylim=c(-0.5,1),xlim =c(-1,1), col="blue") points(c(1/sqrt(2)),c(1/sqrt(2)),col="blue") lines(x2,y2,lwd=3,col="red") points(c(-1/sqrt(2)),c(1/sqrt(2)),col="red") grid(lwd=3) # Transforming iid normal data A = matrix(c(1/2,1/2,1/sqrt(2),-1/sqrt(2)),ncol=2,byrow=T) X = matrix(rnorm(400,mean = 3,sd=0.5),nrow=2) Y = A%*%X #par(mfrow=c(1,1)) plot(X[1,],X[2,],xlab="X1",ylab="X2",col="red",xlim = c(-2,5),ylim = c(-2,5)) lines(c(3,3),c(-2,5),col="green",lwd=2) lines(c(-2,5),c(3,3),col="green",lwd=2) plot(Y[1,],Y[2,],xlab="Y1",ylab="Y2",col="blue",xlim = c(-2,5),ylim = c(-2,5)) lines(c(3,3),c(-2,5),col="green",lwd=2) lines(c(-2,5),c(0,0),col="green",lwd=2) ################################### # normal quantiles (for n = 1000 replications) nq = rep(0,1000) for(i in c(1:1000)){ nq[i] = qnorm((i-0.5)/1000) } # t-quantiles (df=2) nqt = rep(0,1000) for(i in c(1:1000)){ nqt[i] = qt((i-0.5)/1000,df=2) } # remove the last observation # Now make QQplot againt a standard normal # Now construct a t-statisic library(MASS) ttestvec = rep(0,1000) Sigma = matrix(c(1,0,0,0,1,0,0,0,1),nrow=3) for(i in c(1:1000)){ test = mvrnorm(1,c(0,0,0),Sigma) est.mean = mean(test) est.sd = sd(test) ttestvec[i] = sqrt(3)*est.mean/est.sd } # sort estimates from smallest to largest xsort = sort(ttestvec) ## QQplot for standard normal plot(nq[-c(c(1:5),c(996:1000))],xsort[-c(c(1:5),c(996:1000))], type='p', main="",xlab="normal quantiles",ylab="quantiles of t-test") lines(nq[-c(c(1:5),c(996:1000))],nq[-c(c(1:5),c(996:1000))]) # QQplot for t-dist plot(nqt[-c(c(1:5),c(996:1000))],xsort[-c(c(1:5),c(996:1000))], type='p', main="",xlab="t-quantiles quantiles",ylab="quantiles of t-test",lwd=2) lines(nqt[-c(c(1:5),c(996:1000))],nqt[-c(c(1:5),c(996:1000))]) ################################################################################################### # Section 2.4.2 # Projections by demeaning Dimension 2 # Suppose that (X_{1},X_{2}) are iid with mean 5 and variance one. library("MASS") n=10 Sigma2 = matrix(c(1,0,0,1),ncol = 2) xtemp2 = t(mvrnorm(n,c(5,5),Sigma2)) ytemp2 = xtemp2 n = length(xtemp2[1,]) for(i in c(1:n)){ ytemp2[,i] = xtemp2[,i] - mean(xtemp2[,i]) } s=2 plot(x=xtemp2[1,s],y=xtemp2[2,s],lwd=4,type="o",xlab="X1",ylab="X2",xlim=c(-10,10), ylim=c(-10,10),asp=1) #plot(1,1,pch=16,ylim=c(-3,10),asp=1) lines(c(-10,10),c(-10,10),lwd=2,col="red",lty=2) u=mean(xtemp2[,s]) points(u,u, lwd=2,col="red",lty=2) lines(c(-10,10),c(10,-10),lwd=2,col="blue",lty=2) points(ytemp2[1,s],ytemp2[2,s],col="blue",lwd=3,pch=16) ########################################################################## # Section 2.4.3 # Chisq2 chi2 = function(x){ n = length(x) y = x for(i in c(1:n)){ y[i] = (1/2)*exp(-x[i]/2) } return(y) } x = seq(0,10,length=200) y = chi2(x) plot(x,y,col="blue",ylab = "chi square",lwd=3,type="l") x1 = seq(0,40,length=200) sigma2=1/2 ysigma = (2/sigma2)*chi2(x1) sigmax = (1/(2/sigma2))*x1 plot(sigmax,ysigma,col="red",ylab = "sigma chi square",lwd=3,type="l",xlab = "sigma*x",main = "sigma squared = 1/2,n=3") lines(sigmax,y,col="blue",lwd=3) ###################################### # Chisq1 chi1 = function(x){ n = length(x) y = x for(i in c(1:n)){ y[i] = (1/(sqrt(2)*gamma(1/2)))*(x[i]**(-1/2))*exp(-x[i]/2) } return(y) } x = seq(0.01,6,length=100) y = chi1(x) plot(x,y,col="blue",ylab = "chi square",lwd=3,type="l",xlab = "x",main = "Density of Chi-square with 1 df") ########################################3 # chisquares K chiK = function(x,k){ n = length(x) y = x for(i in c(1:n)){ y[i] = (1/((2**(k/2))*gamma(k/2)))*(x[i]**((k/2)-1))*exp(-x[i]/2) } return(y) } x = seq(0.01,25,length=100) x2 = seq(0.01,70,length=100) x3 = seq(0.01,140,length=100) k1 = 10 k2=25 k3 = 50 y = chiK(x,k=k1) y25 = chiK(x2,k=k2) sigma2=1/2 ysigma = (k1/sigma2)*chiK(x,k=k1) ysigma2 = (k2/sigma2)*chiK(x2,k=k2) ysigma3 = (k3/sigma2)*chiK(x3,k=k3) sigmax = (1/(k1/sigma2))*x sigmax2 = (1/(k2/sigma2))*x2 sigmax3 = (1/(k3/sigma2))*x3 plot(sigmax3,ysigma3,col="blue",ylab = "sigma chi square",lwd=3,type="l",xlab = "sigma*x",main = "sigma squared = 1/2,n=11,26 and 51") lines(sigmax2,ysigma2,col="red",lwd=3) lines(sigmax,ysigma,col="green",lwd=3)