4 Working With Matrices in R

For the computational assignment concerning pmfs of compound distributions it is needed to be able to work with matrices in R. Below I talk about some of functions in R that you may find useful.

4.1 Creating Matrices

  1. Creating an Identity Matrix
diag(5)
##      [,1] [,2] [,3] [,4] [,5]
## [1,]    1    0    0    0    0
## [2,]    0    1    0    0    0
## [3,]    0    0    1    0    0
## [4,]    0    0    0    1    0
## [5,]    0    0    0    0    1
  1. Arbitrary Diagonal Matrix
A=diag(1:5)
A
##      [,1] [,2] [,3] [,4] [,5]
## [1,]    1    0    0    0    0
## [2,]    0    2    0    0    0
## [3,]    0    0    3    0    0
## [4,]    0    0    0    4    0
## [5,]    0    0    0    0    5
class(A)
## [1] "matrix" "array"
  1. Using the matrix Function
matrix(1:9,nrow=3)
##      [,1] [,2] [,3]
## [1,]    1    4    7
## [2,]    2    5    8
## [3,]    3    6    9
matrix(1:9,nrow=3,byrow=TRUE)
##      [,1] [,2] [,3]
## [1,]    1    2    3
## [2,]    4    5    6
## [3,]    7    8    9
  1. Brute Force
rbind(1:3,4:6,7:9)
##      [,1] [,2] [,3]
## [1,]    1    2    3
## [2,]    4    5    6
## [3,]    7    8    9
cbind(1:3,4:6,7:9)
##      [,1] [,2] [,3]
## [1,]    1    4    7
## [2,]    2    5    8
## [3,]    3    6    9

4.2 Matrix Algebraic Operations and Manipulations

  1. Column or Row Permutation and Transpose of a Matrix
A=diag(1:5);
# Row Permutation 
A[5:1,]
##      [,1] [,2] [,3] [,4] [,5]
## [1,]    0    0    0    0    5
## [2,]    0    0    0    4    0
## [3,]    0    0    3    0    0
## [4,]    0    2    0    0    0
## [5,]    1    0    0    0    0
# Column Permutation
A[,5:1]
##      [,1] [,2] [,3] [,4] [,5]
## [1,]    0    0    0    0    1
## [2,]    0    0    0    2    0
## [3,]    0    0    3    0    0
## [4,]    0    4    0    0    0
## [5,]    5    0    0    0    0
# Transpose 
t(A[,5:1])
##      [,1] [,2] [,3] [,4] [,5]
## [1,]    0    0    0    0    5
## [2,]    0    0    0    4    0
## [3,]    0    0    3    0    0
## [4,]    0    2    0    0    0
## [5,]    1    0    0    0    0
  1. Extracting the Main Diagonal or a Column or a Row or anyset of elements
A=diag(1:5)
A
##      [,1] [,2] [,3] [,4] [,5]
## [1,]    1    0    0    0    0
## [2,]    0    2    0    0    0
## [3,]    0    0    3    0    0
## [4,]    0    0    0    4    0
## [5,]    0    0    0    0    5
# Main Diagnoal 
diag(A)
## [1] 1 2 3 4 5
# 2nd Column
A[1:5,2]
## [1] 0 2 0 0 0
# 4th Row
A[4,1:5]
## [1] 0 0 0 4 0
# Cross Diagonal
A[cbind(1:5,5:1)]
## [1] 0 0 3 0 0
  1. Matrix Multiplication
A=matrix(1:9,nrow=3);
A
##      [,1] [,2] [,3]
## [1,]    1    4    7
## [2,]    2    5    8
## [3,]    3    6    9
B=diag(1:3);
B
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    2    0
## [3,]    0    0    3
# Scaling Columns
A%*%B
##      [,1] [,2] [,3]
## [1,]    1    8   21
## [2,]    2   10   24
## [3,]    3   12   27
# Scaling Rows
B%*%A
##      [,1] [,2] [,3]
## [1,]    1    4    7
## [2,]    4   10   16
## [3,]    9   18   27
  1. Matrix Addition
A+B
##      [,1] [,2] [,3]
## [1,]    2    4    7
## [2,]    2    7    8
## [3,]    3    6   12
A-B
##      [,1] [,2] [,3]
## [1,]    0    4    7
## [2,]    2    3    8
## [3,]    3    6    6

4.3 Matrix Functions

  1. Eigen Values, Eigen Vectors and Inverse
A<-matrix(1:9,nrow=3);
A<-A%*%t(A) # A Positive Semi Definite Matrix
# Eigen Values
eigen(A)$values
## [1] 2.838586e+02 1.141413e+00 3.739411e-14
# Eigen Vectors 
eigen(A)$vectors # Each Row is an Eigen Vector
##            [,1]        [,2]       [,3]
## [1,] -0.4796712  0.77669099  0.4082483
## [2,] -0.5723678  0.07568647 -0.8164966
## [3,] -0.6650644 -0.62531805  0.4082483
# Check
round(A%*%eigen(A)$vectors[,1]-eigen(A)$values[1]*eigen(A)$vectors[,1],5)
##      [,1]
## [1,]    0
## [2,]    0
## [3,]    0
# Another Check
round(t(eigen(A)$vectors)%*%eigen(A)$vectors,5)
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    1    0
## [3,]    0    0    1
# Cholesky Decomposition
chol(A+0.5*diag(3))
##          [,1]     [,2]      [,3]
## [1,] 8.154753 9.564974 11.036508
## [2,] 0.000000 1.418195  1.717739
## [3,] 0.000000 0.000000  1.320931
# Inverse
solve(A+0.5*diag(3))
##             [,1]       [,2]        [,3]
## [1,]  0.70166042 -0.6298875  0.03856454
## [2,] -0.62988752  1.3379754 -0.69416176
## [3,]  0.03856454 -0.6941618  0.57311194
  1. Determinant and Trace and Some Fun
# A Random Matrix 
A<-matrix(rnorm(9),nrow=3);
A<-A%*%t(A)
# Is it Positive Definite
eigen(A)$values
## [1] 5.435727942 3.928923632 0.000224315
# Determinant and Trace
c(det(A),sum(diag(A)))
## [1] 0.004790597 9.364875890
# Check
c(prod(eigen(A)$values),sum(eigen(A)$values))
## [1] 0.004790597 9.364875890

Page built: 2021-04-09