############################################################################ # NAME: Chris Bilder # # DATE: 3-22-16 # # PURPOSE: Show the basics of matrix algebra in R # # Program is based on one used for STAT 873 by the same name # # # ############################################################################ ############################################################################ # Basics A <- matrix(data = c(1, 2, 3, 4, 5, 6), nrow = 2, ncol = 3, byrow = TRUE) A class(A) B <- matrix(data = c(-1, 10, -1, 5, 5, 8), nrow = 2, ncol = 3, byrow = TRUE) A[1,1] A[1,] A[,1] A + B A - B #Show what happens with the byrow = TRUE option matrix(data = c(1, 2, 3, 4, 5, 6), nrow = 2, ncol = 3) #Vector example y <- matrix(data = c(1,2,3), nrow = 3, ncol = 1, byrow = TRUE) y class(y) x <- c(1,2,3) x class(x) is.vector(x) ############################################################################# # Matrix multiplication A <- matrix(data = c(1, 2, 3, 4, 5, 6), nrow = 2, ncol = 3, byrow = TRUE) B <- matrix(data = c(3, 0, 1, 2, 0, 1), nrow = 3, ncol = 2, byrow = TRUE) A %*% B B %*% A #Of course, you can save the results in a matrix too C <- A %*% B C #What is A*B? A * B #This shows that element-by-element multication take place A * A 3 * A #Work with vectors x <- c(1,2,3) A %*% x # Inner product x %*% x # Outer product x %o% x ############################################################################# # Inverse A <- matrix(data = c(1, 2, 3, 4), nrow = 2, ncol = 2, byrow = TRUE) solve(A) A%*%solve(A) solve(A)%*%A round(solve(A)%*%A, 2) ############################################################################# # More matrix operations # Find diagonal elements of a matrix A <- matrix(data = c(1, 2, 3, 4), nrow = 2, ncol = 2, byrow = TRUE) A diag(A) # Create diagonal matrix diag(c(1,2,3,4)) # Determinant det(A) # Transpose t(A) ############################################################################# # Eigenvalues A <- matrix(data = c(1, 0.5, 0.5, 1.25), nrow = 2, ncol = 2, byrow = TRUE) A eigen(A) save.eig <- eigen(A) names(save.eig) save.eig$values save.eig$vectors # Verify A%*%save.eig$vectors[,1] save.eig$values[1]*save.eig$vectors[,1] # Verify A%*%save.eig$vectors[,2] save.eig$values[2]*save.eig$vectors[,2] # Length of a vector sqrt(sum(save.eig$vectors[,1]^2)) sqrt(sum(save.eig$vectors[,2]^2)) ############################################################ # Eigenvector plot # Square plot (default is "m" which means maximal region) par(pty = "s") # Set up some dummy values for plot b1 <- c(-1,1) b2 <- c(-1,1) plot(x = b1, y = b2, type = "n", main = "Eigenvectors of A", xlab = expression(b[1]), ylab = expression(b[2]) , panel.first = grid()) # Draw line on plot - h specifies a horizontal line abline(h = 0, lty = "solid", lwd = 2) # v specifies a vertical line abline(v = 0, lty = "solid", lwd = 2) # Draw eigenvectors arrows(x0 = 0, y0 = 0, x1 = save.eig$vectors[1,1], y1 = save.eig$vectors[2,1], col = "red", lty = "solid") arrows(x0 = 0, y0 = 0, x1 = save.eig$vectors[1,2], y1 = save.eig$vectors[2,2], col = "red", lty = "solid") # Other eigenvectors: new.vec <- -10*save.eig$vectors[,1] A %*% new.vec save.eig$values[1] * new.vec sqrt(sum(new.vec^2)) ############################################################################# # Additional items # Add a leading column of 1's mat1 <- matrix(data = 1:20, nrow = 10, ncol = 2, byrow = TRUE) mat1 X <- cbind(1, mat1) X class(X) #Coerce data frame to matrix set1 <- data.frame(x1 = 1:10, x2 = 11:20) class(set1) mat2 <- as.matrix(set1) mat2 class(mat2) # Include row and column names A <- matrix(data = c(1, 2, 3, 4, 5, 6), nrow = 2, ncol = 3, byrow = TRUE, dimnames = list(my.row.names = c("row1", "row2"), my.col.names = c("col1", "col2", "col2"))) A #